Czy w GLM prawdopodobieństwo dziennika modelu nasyconego zawsze wynosi zero?

Jako część danych wyjściowych uogólnionego modelu liniowego do oceny modelu wykorzystywane są odchylenie zerowe i rezydualne. Często widzę formuły dla tych wielkości wyrażone jako prawdopodobieństwo dziennika modelu nasyconego, na przykład: /stats//a/113022/22199 , Regresja logistyczna: jak uzyskać model nasycony

Model nasycony, o ile rozumiem, jest modelem, który doskonale pasuje do obserwowanej odpowiedzi. Dlatego w większości miejsc, które widziałem, prawdopodobieństwo logarytmiczne nasyconego modelu jest zawsze podawane jako zero.

Jednak sposób, w jaki podawana jest formuła dewiacji, sugeruje, że czasami ta liczba nie jest równa zero. (Jakby zawsze było zero, po co zawracać sobie tym głowę?)

W jakich przypadkach może być różna od zera? Jeśli nigdy nie jest różna od zera, to po co uwzględniać ją w formule dewiacji?

— Alex
źródło

Odpowiedzi:

Jeśli naprawdę miałeś na myśli prawdopodobieństwo dziennika , odpowiedź brzmi: nie zawsze jest to zero.

Weźmy na przykład dane Poissona: . Prawdopodobieństwo dziennika dla jest określone przez: $y_i \sim \text{Poisson}(\mu_i), i = 1, \ldots, n$ $Y = (y_1, \ldots, y_n)$

\begin{matrix} (*) & ℓ (μ; Y) = - \sum_{i = 1}^{n} μ_{i} + \sum_{i = 1}^{n} y_{i} \log μ_{i} - \sum_{i = 1}^{n} \log (y_{i}!) . \end{matrix}

$\ell(\mu; Y) = -\sum_{i = 1}^n \mu_i + \sum_{i = 1}^n y_i \log \mu_i - \sum_{i = 1}^n \log(y_i!). \tag{$*$}$

Zróżnicuj w względem i ustaw na (w ten sposób otrzymujemy MLE dla modelu nasyconego): rozwiązać ten problem za uzyskać zastępując powrotem do dla daje że Log-Likelihood nasyconego modelu to: chyba że weźmiesz specjalnego wartości. $\ell(\mu; Y)$ $(*)$ $\mu_i$ $0$

- 1 + \frac{y_{i}}{μ_{i}} = 0.

$-1 + \frac{y_i}{\mu_i} = 0.$

μ_{i}

$\mu_i$

{\hat{μ}}_{i} = y_{i}

$\hat{\mu}_i = y_i$

{\hat{μ}}_{i}

$\hat{\mu}_i$

(*)

$(*)$

μ_{i}

$\mu_i$

ℓ (\hat{μ}; Y) = \sum_{i = 1}^{n} y_{i} (\log y_{i} - 1) - \sum_{i = 1}^{n} \log (y_{i}!) \neq 0

$\ell(\hat{\mu}; Y) = \sum_{i = 1}^n y_i(\log y_i - 1) -\sum_{i = 1}^n \log(y_i!) \neq 0$

y_{i}

$y_i$

Na stronie pomocy Rfunkcji glmpod pozycją deviancedokument wyjaśnia ten problem w następujący sposób:

deviance do stałej, minus dwukrotność maksymalnego prawdopodobieństwa logarytmu. Tam, gdzie jest to uzasadnione, wybrana jest stała, aby model nasycony miał zero odchyleń.

Zauważ, że wspomniano, że odchylenie zamiast logarytmu prawdopodobieństwa modelu nasyconego jest wybrane na zero.

Prawdopodobnie tak naprawdę chciałeś potwierdzić, że „ odchylenie nasyconego modelu jest zawsze podawane jako zero”, co jest prawdą, ponieważ odchylenie jest z definicji (patrz rozdział 4.5.1 analizy jakościowej danych (wydanie 2) Alana Agresti) to statystyka stosunku prawdopodobieństwa określonego GLM do modelu nasyconego. Wyżej constantwspomniane w dokumentacji R jest w rzeczywistości dwukrotnością maksymalnego logarytmu prawdopodobieństwa modelu nasyconego.

Jeśli chodzi o stwierdzenie „Jednak sposób, w jaki podano formułę dewiacji, sugeruje, że czasami ta liczba nie jest równa zero”, prawdopodobnie wynika to z nadużycia użycia terminu dewiacja . Na przykład, w R, statystyczny współczynnik prawdopodobieństwa porównywania dwóch dowolnych (zagnieżdżone) Modele i jest również określana jako odchyleniu, co będzie bardziej dokładnie określanej jako z różnicy pomiędzy odchyleniu od i odchyleniu od , jeżeli ściśle przestrzegał definicji podanej w książce Agresti. $M_1$ $M_2$ $M_1$ $M_2$

Wniosek

Prawdopodobieństwo logiczne nasyconego modelu jest zasadniczo niezerowe.
Odchylenie (w oryginalnej definicji) modelu nasyconego wynosi zero.
Odchyleniem wyjściowy oprogramowanie (takie jak R) jest na ogół tak niezerowych w rzeczywistości oznacza co innego (różnica pomiędzy odchyleń).

Poniżej przedstawiono pochodną dla ogólnego przypadku rodziny wykładniczej i inny konkretny przykład. Załóżmy, że dane pochodzą z rodziny wykładniczej (patrz Modern Applied Statistics with S , Rozdział ): gdzie są znanymi wcześniejszymi wagami, a to parametr dyspersji / skali (w wielu przypadkach, takich jak dwumianowy i Poisson, ten parametr jest znany, podczas gdy w innych przypadkach, takich jak normalny i Gamma, ten parametr jest nieznany). Następnie prawdopodobieństwo logarytmu podaje: $7$

\begin{matrix} (1) & f (y_{i}; θ_{i}, φ) = \exp [A_{i} (y_{i} θ_{i} - γ (θ_{i})) / φ + τ (y_{i}, φ / A_{i})] . \end{matrix}

$f(y_i; \theta_i, \varphi) = \exp[A_i(y_i\theta_i - \gamma(\theta_i))/\varphi + \tau(y_i, \varphi/A_i)]. \tag{1}$

A_{i}

$A_i$

φ

$\varphi$

ℓ (θ, φ; Y) = \sum_{i = 1}^{n} A_{i} (y_{i} θ_{i} - γ (θ_{i})) / φ + \sum_{i = 1}^{n} τ (y_{i}, φ / A_{i}) .

$\ell(\theta, \varphi; Y) = \sum_{i = 1}^n A_i(y_i \theta_i - \gamma(\theta_i))/\varphi + \sum_{i = 1}^n \tau(y_i, \varphi/A_i).$ Tak jak w przykładzie Poissona, parametry modelu nasyconego można oszacować, rozwiązując następującą funkcję punktacji :

0 = U (θ_{i}) = \frac{\partial ℓ (θ, φ; Y)}{\partial θ_{i}} = \frac{A_{i} (y_{i} - γ^{'} (θ_{i}))}{φ}

$0 = U(\theta_i) = \frac{\partial \ell(\theta, \varphi; Y)}{\partial \theta_i} = \frac{A_i(y_i - \gamma'(\theta_i))}{\varphi}$

Oznacz rozwiązanie powyższego równania wzorem , wtedy ogólna postać logarytmu prawdopodobieństwa modelu nasyconego (traktuj parametr skali jako stałą) to: $\hat{\theta}_i$

\begin{matrix} (* *) & ℓ (\hat{θ}, φ; Y) = \sum_{i = 1}^{n} A_{i} (y_{i} {\hat{θ}}_{i} - γ ({\hat{θ}}_{i})) / φ + \sum_{i = 1}^{n} τ (y_{i}, φ / A_{i}) . \end{matrix}

$\ell(\hat{\theta}, \varphi; Y) = \sum_{i = 1}^n A_i(y_i \hat{\theta}_i - \gamma(\hat{\theta}_i))/\varphi + \sum_{i = 1}^n \tau(y_i, \varphi/A_i). \tag{$**$}$

W mojej poprzedniej odpowiedzi błędnie stwierdziłem, że pierwszy termin po prawej stronie jest zawsze równy zero, powyższy przykład danych Poissona dowodzi, że jest on błędny. Aby uzyskać bardziej skomplikowany przykład, rozważ rozkład gamma podany w dodatku. $(**)$ $\Gamma(\alpha, \beta)$

Dowód pierwszego terminu w logarytmicznym prawdopodobieństwie nasycenia modelu gamma jest niezerowy : Biorąc pod uwagę musimy najpierw przeprowadzić ponowną parametryzację, aby miało wykładniczą postać rodziny . Można to zweryfikować, jeśli pozwalając wówczas ma reprezentację: gdzie

f (y; α, β) = \frac{β^{α}}{Γ (α)} e^{- β y} y^{α - 1}, y > 0, α > 0, β > 0,

$f(y; \alpha, \beta) = \frac{\beta^\alpha}{\Gamma(\alpha)}e^{-\beta y}y^{\alpha - 1}, \quad y > 0, \alpha > 0, \beta > 0,$

f

$f$

(1)

$(1)$

φ = \frac{1}{α}, θ = - \frac{β}{α},

$\varphi = \frac{1}{\alpha},\, \theta = -\frac{\beta}{\alpha},$

f

$f$

f (y; θ, φ) = \exp [\frac{θ y - (- \log (- θ))}{φ} + τ (y, φ)],

$f(y; \theta, \varphi) = \exp\left[\frac{\theta y - (-\log(-\theta))}{\varphi}+ \tau(y, \varphi)\right],$

τ (y, φ) = - \frac{\log φ}{φ} + (\frac{1}{φ} - 1) \log y - \log Γ (φ^{- 1}) .

$\tau(y, \varphi) = -\frac{\log \varphi}{\varphi} + \left(\frac{1}{\varphi} - 1\right)\log y - \log\Gamma(\varphi^{-1}).$ Dlatego MLE modelu nasyconego to . Stąd chyba że przyjmą bardzo specjalne wartości.

{\hat{θ}}_{i} = - \frac{1}{y_{i}}

$\hat{\theta}_i = -\frac{1}{y_i}$

\sum_{i = 1}^{n} \frac{1}{φ} [{\hat{θ}}_{i} y_{i} - (- \log (- {\hat{θ}}_{i}))] = \sum_{i = 1}^{n} \frac{1}{φ} [- 1 - \log (y_{i})] \neq 0,

$\sum_{i = 1}^n \frac{1}{\varphi}[\hat{\theta}_iy_i - (-\log(-\hat{\theta}_i))] = \sum_{i = 1}^n \frac{1}{\varphi}[-1 - \log(y_i)] \neq 0,$

y_{i}

$y_i$

— Zhanxiong
źródło

Czy logarytmiczność wynosi zero wtedy i tylko wtedy, gdy model może przypisać 100% prawdopodobieństwa każdemu z możliwych wyników?

— Alex

Nie do końca rozumiem, co miałeś na myśli. Ale z mojej pochodnej można wywnioskować, że wynosi wtedy i tylko wtedy, gdy jest identycznie i nie ma parametru dyspersji.

0

$0$

τ

$\tau$

0

$0$

— Zhanxiong,

Twoje pochodzenie jest bardzo dobre, ale formalny dowód jest obecnie nieco ponad moją głową. Dziękujemy za przykład z modelem Poissona. Z tego przykładu wyjąłem to, że model Poissona nie może przypisać 100% prawdopodobieństwa obserwowanemu wynikowi, biorąc pod uwagę jakąkolwiek wartość średniej Poissona, a zatem prawdopodobieństwo nie może wynosić zero.

— Alex

Stwierdzenie „model przypisuje prawdopodobieństwa obserwowanemu wynikowi” wydaje mi się dziwne. Czy masz na myśli to, że biorąc pod uwagę obserwacje , a jeśli jest zmienną losową Poissona, ?

100 %

$100\%$

y_{1}, \dots, y_{n}

$y_1, \ldots, y_n$

Y

$Y$

P (Y = y_{1}) + P (Y = y_{2}) + \dots + P (Y = y_{n}) < 1

$P(Y= y_1) + P(Y = y_2) + \cdots + P(Y = y_n) < 1$

— Zhanxiong,

Chodzi mi o to, że jeśli jest zmienną losową Poissona, to dla dowolnej wartości lub Poissona, dlatego nie można znaleźć żadnego parametru modelu, który dałby logarytmiczne prawdopodobieństwo zerowe dla obserwowanego . Może zupełnie nie rozumiem koncepcji nasyconego modelu.

Y

$Y$

P (Y = y_{i}) < 1

$P(Y = y_i) < 1$

i

$i$

— Alex

Odpowiedź Zhanxionga jest już świetna (+1), ale tutaj jest szybki pokaz, że prawdopodobieństwo logarytmu modelu nasyconego wynosi dla regresji logistycznej. Pomyślałem, że opublikuję, ponieważ nie widziałem tego TeX-a na tej stronie i ponieważ właśnie napisałem to na wykład. $0$

Prawdopodobieństwo to gdzie .

\begin{matrix} (1) & L (y; X, β) = \prod_{i = 1}^{n} f (y_{i}; x_{i}, β) = \prod_{i = 1}^{n} π_{i}^{y_{i}} (1 - π_{i})^{1 - y_{i}} = \prod_{i = 1}^{n} {(\frac{π_{i}}{1 - π_{i}})}^{y_{i}} (1 - π_{i}) \end{matrix}

$L(\mathbf{y} ; \mathbf{X}, \boldsymbol{\beta}) = \prod_{i=1}^n f(y_i ; \mathbf{x}_i, \boldsymbol{\beta}) = \prod_{i=1}^n \pi_i^{y_i}(1-\pi_i)^{1-y_i} = \prod_{i=1}^n\left( \frac{\pi_i}{1-\pi_i}\right)^{y_i} (1 - \pi_i) \tag{1}$

π_{i} = invlogit (x_{i}^{⊺} β)

$\pi_i = \text{invlogit}(\mathbf{x}_i^\intercal \boldsymbol{\beta} )$

Prawdopodobieństwo dziennika to

\begin{aligned} \log L (y; X, β) & = \sum_{i = 1}^{n} y_{i} \log (\frac{π_{i}}{1 - π_{i}}) + \log (1 - π_{i}) \\ = \sum_{i = 1}^{n} y_{i} logit (π_{i}) + \log (1 - π_{i}) \\ = \sum_{i = 1}^{n} y_{i} x_{i}^{⊺} β + \log (1 - invlogit (x_{i}^{⊺} β)) \\ = \sum_{i = 1}^{n} y_{i} x_{i}^{⊺} β + \log (invlogit (- x_{i}^{⊺} β)) \\ = \sum_{i = 1}^{n} y_{i} x_{i}^{⊺} β - \log (1 + \exp [x_{i}^{⊺} β])) \end{aligned}

$\begin{align*} \log L(\mathbf{y} ; \mathbf{X}, \boldsymbol{\beta}) &= \sum_{i=1}^n y_i \log \left( \frac{\pi_i}{1-\pi_i}\right) + \log(1-\pi_i) \\ &= \sum_{i=1}^n y_i \text{logit} \left( \pi_i \right) + \log(1-\pi_i) \\ &= \sum_{i=1}^n y_i \mathbf{x}_i^\intercal \boldsymbol{\beta} + \log( 1 - \text{invlogit}(\mathbf{x}_i^\intercal \boldsymbol{\beta} )) \\ &= \sum_{i=1}^n y_i \mathbf{x}_i^\intercal \boldsymbol{\beta} + \log( \text{invlogit}( - \mathbf{x}_i^\intercal \boldsymbol{\beta} )) \\ &= \sum_{i=1}^n y_i \mathbf{x}_i^\intercal \boldsymbol{\beta} - \log( 1 + \exp[ \mathbf{x}_i^\intercal \boldsymbol{\beta}] )) \end{align*}$

Jeśli weźmiesz pochodne w odniesieniu do wszystkich współczynników, otrzymasz

\begin{matrix} (2) & \nabla ℓ (β) = \sum_{i = 1}^{n} y_{i} x_{i} - \frac{\exp [x_{i}^{⊺} β]}{(1 + \exp [x_{i}^{⊺} β])} x_{i} . \end{matrix}

$\nabla \ell(\boldsymbol{\beta}) = \sum_{i=1}^n y_i \mathbf{x}_i - \frac{\exp[ \mathbf{x}_i^\intercal \boldsymbol{\beta}]}{( 1 + \exp[ \mathbf{x}_i^\intercal \boldsymbol{\beta}] ) }\mathbf{x}_i \tag{2}.$

Ustawienie tego wyrażenia na i rozwiązanie da ci odpowiedź. Zwykle nie można tego zrobić analitycznie, co tłumaczy popularność / konieczność stosowania algorytmów iteracyjnych w celu dopasowania do tego modelu, ale w przypadku modelu nasyconego jest to możliwe. $\mathbf{0}$ $\boldsymbol{\beta}$

Aby znaleźć model nasycony, podajemy każdy wiersz jego własny współczynnik. Więc a macierz projektowa razy wektor współczynnika to $\boldsymbol{\beta} \in \mathbb{R}^n$

X β = [\begin{matrix} 1 & 0 & \dots & 0 \\ 0 & 1 & \dots & 0 \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 0 & \dots & 1 \end{matrix}] [\begin{matrix} β_{1} \\ β_{2} \\ ⋮ \\ β_{n} \end{matrix}] .

$\mathbf{X}\boldsymbol{\beta} = \begin{bmatrix} 1 & 0 & \cdots & 0\\ 0 & 1 & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1\\ \end{bmatrix} \begin{bmatrix} \beta_1 \\ \beta_2 \\ \vdots \\ \beta_n \end{bmatrix}.$

Zwróć uwagę, że w szczególności . $\mathbf{x}_i^\intercal \boldsymbol{\beta} = \beta_i$

Biorąc ty wiersz równania (2) daje nam $j$

\sum_{i = 1}^{n} y_{i} x_{i, j} = \sum_{i = 1}^{n} \frac{\exp [x_{i}^{⊺} β]}{(1 + \exp [x_{i}^{⊺} β])} x_{i, j}

$\sum_{i=1}^n y_i x_{i,j} = \sum_{i=1}^n\frac{\exp[ \mathbf{x}_i^\intercal \boldsymbol{\beta}]}{( 1 + \exp[ \mathbf{x}_i^\intercal \boldsymbol{\beta}] ) }x_{i,j}$

który może być tylko prawda, jeśli dla każdej obserwacji : $i$

y_{i} = invlogit (β_{i})

$y_i = \text{invlogit}(\beta_i )$ lub innymi słowy każdy ma wartość plus lub minus nieskończoność (jeśli wynosi odpowiednio lub ). Możemy podłączyć te parametry z powrotem do (1), aby uzyskać maksymalne prawdopodobieństwo: Oczywiście log tego wynosi .

β_{i}

$\beta_i$

y_{i}

$y_i$

1

$1$

0

$0$

\prod_{i = 1}^{n} {\hat{π}}_{i}^{y_{i}} (1 - {\hat{π}}_{i})^{1 - y_{i}} = 1^{n} = 1.

$\prod_{i=1}^n \hat{\pi}_i^{y_i}(1-\hat{\pi}_i)^{1-y_i} = 1^n = 1.$

0

$0$

— Taylor
źródło

Ale zakłada to niepogrupowane dane . Jeśli masz grupy o (i tych samych wartościach zmiennych towarzyszących) (w R, na przykład przy użyciu formularza ), wówczas model nasycony nie ma logarytmu zero.

n_{i} > 1

$n_i>1$ glm( cbind(k, n-k) ~ x + ...

— kjetil b halvorsen

@kjetilbhalvorsen oh good point. Nigdy nie próbowałem, żebym mógł sprawdzić

— Taylor

@Alex: tak, zgadza się. przynajmniej dla dyskretnych rozkładów. dla ciągłych rozkładów sprowadzałoby się to do dopuszczenia, aby gęstość była równa 1, co niekoniecznie jest znaczące, a zatem nie jest rozsądną rzeczą do osiągnięcia. nieco bardziej ogólnie, logarytmiczne prawdopodobieństwo modelu nasyconego daje górną granicę wydajności dowolnego modelu, który jest zgodny z założeniem podstawowej rodziny dystrybucji. Innymi słowy, prawdopodobieństwo logarytmiczne nasyconego modelu dwumianowego jest „tak dobre, jak to możliwe” dla danego zestawu danych (X, Y) przy założeniu, że Y jest dwumianowy. Sensowne jest porównanie twojego modelu glm z tą górną granicą w przeciwieństwie do, powiedzmy, 100% (lub podobnego), ponieważ twój model jest z natury ograniczony przez twoje założenie dotyczące rozkładu odpowiedzi.

— bettmensch88
źródło