Kiedy Poisson i ujemne regresje dwumianowe pasują do tych samych współczynników?

Zauważyłem, że w regresji R, Poissona i regresji dwumianowej ujemnej (NB) zawsze wydaje się pasować do tych samych współczynników dla predyktorów jakościowych, ale nie ciągłych.

Na przykład oto regresja z predyktorem jakościowym:

data(warpbreaks)
library(MASS)

rs1 = glm(breaks ~ tension, data=warpbreaks, family="poisson")
rs2 = glm.nb(breaks ~ tension, data=warpbreaks)

#compare coefficients
cbind("Poisson"=coef(rs1), "NB"=coef(rs2))

Oto przykład z ciągłym predyktorem, w którym Poisson i NB pasują do różnych współczynników:

data(cars)
rs1 = glm(dist ~ speed, data=cars, family="poisson")
rs2 = glm.nb(dist ~ speed, data=cars)

#compare coefficients
cbind("Poisson"=coef(rs1), "NB"=coef(rs2))

(Oczywiście nie są to dane zliczające, a modele nie mają znaczenia ...)

Następnie przekoduję predyktor na czynnik, a dwa modele ponownie będą pasować do tych samych współczynników:

library(Hmisc)
speedCat = cut2(cars$speed, g=5) 
#you can change g to get a different number of bins

rs1 = glm(cars$dist ~ speedCat, family="poisson")
rs2 = glm.nb(cars$dist ~ speedCat)

#compare coefficients
cbind("Poisson"=coef(rs1), "NB"=coef(rs2))

Jednak ujemna regresja dwumianowa Josepha Hilbe daje przykład (6.3, str. 118-119), w którym predyktor kategoryczny, płeć, ma nieco inne współczynniki według Poissona ( ) i NB ( ). Mówi: „Wskaźniki częstości występowania między modelami Poissona i NB są dość podobne. Nie jest to zaskakujące, biorąc pod uwagę bliskość [odpowiadającego w R] do zera. ”b = 0,881 α 1 /$b=0.883$$b=0.881$$\alpha$$1/\theta$

Rozumiem to, ale w powyższych przykładach summary(rs2)mówi nam, że jest szacowana odpowiednio na 9,16 i 7,93.$\theta$

Dlaczego więc współczynniki są dokładnie takie same? A dlaczego tylko dla predyktorów jakościowych?

---

Edytuj nr 1

Oto przykład z dwoma nieortogonalnymi predyktorami. Rzeczywiście współczynniki nie są już takie same:

data(cars)

#make random categorical predictor
set.seed(1); randomCats1 = sample( c("A","B","C"), length(cars$dist), replace=T)
set.seed(2); randomCats2 = sample( c("C","D","E"), length(cars$dist), replace=T)

rs1 = glm(dist ~ randomCats1 + randomCats2, data=cars, family="poisson")
rs2 = glm.nb(dist ~ randomCats1 + randomCats2, data=cars)

#compare coefficients
cbind("Poisson"=coef(rs1), "NB"=coef(rs2))

I włączenie innego predyktora powoduje, że modele pasują do różnych współczynników, nawet gdy nowy predyktor jest ciągły. Czy ma to coś wspólnego z ortogonalnością zmiennych zastępczych, które utworzyłem w moim oryginalnym przykładzie?

rs1 = glm(dist ~ randomCats1 + speed, data=cars, family="poisson")
rs2 = glm.nb(dist ~ randomCats1 + speed, data=cars)

#compare coefficients
cbind("Poisson"=coef(rs1), "NB"=coef(rs2))

Odkryłeś intymną, ale ogólną właściwość GLM pasującą według maksymalnego prawdopodobieństwa . Wynik zanika, gdy weźmie się pod uwagę najprostszy ze wszystkich przypadków: dopasowanie jednego parametru do jednej obserwacji!

Odpowiedź w jednym zdaniu : Jeśli zależy nam tylko na oddzielnych środkach do rozłączania podzbiorów naszej próbki, wówczas GLM zawsze da dla każdego podzestawu , więc rzeczywista struktura błędu i parametryzacja gęstości stają się oba nie ma znaczenia dla (punktowego) oszacowania!$\hat\mu_j = \bar y_j$$j$

Trochę więcej : dopasowanie ortogonalnych czynników kategorialnych według maksymalnego prawdopodobieństwa jest równoważne dopasowaniu osobnych środków do rozłącznych podzbiorów naszej próbki, więc wyjaśnia to, dlaczego Poisson i ujemne dwumianowe GLM dają te same oszacowania parametrów. Rzeczywiście, to samo dzieje się bez względu na to, czy stosujemy regresję Poissona, negbina, gaussa, odwrotną regresję gaussowską czy gamma (patrz poniżej). W przypadku Poissona i Negbina domyślną funkcją łącza jest link , ale jest to czerwony śledź; chociaż daje to takie same surowe oszacowania parametrów, zobaczymy poniżej, że ta właściwość naprawdę nie ma nic wspólnego z funkcją link.$\log$

Gdy interesuje nas parametryzacja o większej strukturze lub zależna od ciągłych predyktorów, wówczas przyjęta struktura błędu staje się istotna ze względu na relację średniej wariancji rozkładu, ponieważ odnosi się do parametrów i funkcji nieliniowej zastosowanej do modelowania warunkowego znaczy.

GLM i rodziny dyspersji wykładniczej: kurs zderzeniowy

Wykładniczy rodziny dyspersja w naturalnej postaci jest taki, że gęstość dziennika ma postać

$$\log f(y;\,\theta,\nu) = \frac{\theta y - b(\theta)}{\nu} + a(y,\nu) \>.$$

Tutaj jest parametrem naturalnym, a jest parametrem dyspersji . Gdyby były znane, byłaby to tylko standardowa jednoparametrowa rodzina wykładnicza. Wszystkie GLM rozważane poniżej zakładają model błędu z tej rodziny.$\theta$$\nu$$\nu$

Rozważ próbkę pojedynczej obserwacji z tej rodziny. Jeśli dopasujemy według maksymalnego prawdopodobieństwa, otrzymamy, że , niezależnie od wartości . Rozciąga się to na przypadek próbki iid, ponieważ prawdopodobieństwa dziennika są dodawane, dając .$\theta$$y = b'(\hat\theta)$$\nu$$\bar y = b'(\hat\theta)$

Ale wiemy również, ze względu na ładną regularność gęstości logów w funkcji , że Tak więc w rzeczywistości .$\theta$

$$\frac{\partial}{\partial \theta} \mathbb E \log f(Y;\theta,\nu) = \mathbb E \frac{\partial}{\partial \theta} \log f(Y;\theta,\nu) = 0 \>.$$ $b'(\theta) = \mathbb E Y = \mu$

Ponieważ szacunki maksymalnego prawdopodobieństwa są niezmienne podczas transformacji, oznacza to, że dla tej rodziny gęstości.$\bar y = \hat\mu$

Teraz w GLM modelujemy jako gdzie jest funkcją łącza. Ale jeśli jest wektorem wszystkich zer, z wyjątkiem pojedynczego 1 w pozycji , to . Prawdopodobieństwo GLM jest następnie rozkładane na czynniki według i postępujemy jak wyżej. Dokładnie tak jest w przypadku czynników ortogonalnych.$\mu_i$$\mu_i = g^{-1}(\mathbf x_i^T \beta)$$g$$\mathbf x_i$$j$$\mu_i = g(\beta_j)$$\beta_j$

Co takiego różni się w ciągłych predyktorach?

Gdy predyktory są ciągłe lub mają charakter kategoryczny, ale nie można ich sprowadzić do postaci ortogonalnej, wówczas prawdopodobieństwo nie uwzględnia już poszczególnych terminów z oddzielną średnią w zależności od osobnego parametru. W tym momencie, struktura i funkcja błędu Link nie wchodzą w grę.

Jeśli przejdziemy przez (nużącą) algebrę, równania prawdopodobieństwa staną się dla wszystkich gdzie . Tutaj parametry i wchodzą domyślnie poprzez relację łącza i wariancję .j = 1 , … , p λ i = x T i β β ν μ i = g ( λ i ) = g ( x T i β ) σ 2

$$\sum_{i=1}^n \frac{(y_i - \mu_i)x_{ij}}{\sigma_i^2}\frac{\partial \mu_i}{\partial \lambda_i} = 0\>,$$ $j = 1,\ldots,p$$\lambda_i = \mathbf x_i^T \beta$$\beta$$\nu$$\mu_i = g(\lambda_i) = g(\mathbf x_i^T \beta)$$\sigma_i^2$

W ten sposób funkcja łącza i przyjęty model błędu stają się istotne dla oszacowania.

Przykład: model błędu (prawie) nie ma znaczenia

W poniższym przykładzie generujemy ujemne losowe dane dwumianowe w zależności od trzech czynników kategorialnych. Każda obserwacja pochodzi z jednej kategorii i zastosowano ten sam parametr dyspersji ( ).$k = 6$

Następnie dopasowujemy się do tych danych przy użyciu pięciu różnych GLM, każdy z linkiem : ( a ) dwumian ujemny, ( b ) Poisson, ( c ) gaussowski, ( d ) odwrotny gaussowski i ( e ) gamma GLM. Wszystkie są przykładami wykładniczych rodzin dyspersyjnych.$\log$

Z tabeli wynika, że ​​oszacowania parametrów są identyczne , nawet jeśli niektóre z tych GLM dotyczą danych dyskretnych, a inne ciągłe, a niektóre są danymi nieujemnymi, a inne nie.

      negbin  poisson gaussian invgauss    gamma
XX1 4.234107 4.234107 4.234107 4.234107 4.234107
XX2 4.790820 4.790820 4.790820 4.790820 4.790820
XX3 4.841033 4.841033 4.841033 4.841033 4.841033

Zastrzeżenie w nagłówku wynika z faktu, że procedura dopasowania zakończy się niepowodzeniem, jeśli obserwacje nie mieszczą się w zakresie określonej gęstości. Na przykład, gdybyśmy mieli liczb losowo generowanych w powyższych danych, wówczas Gamma GLM nie zbiegałby się, ponieważ Gamma GLM wymagają ściśle dodatnich danych.$0$

Przykład: Funkcja link (prawie) nie ma znaczenia

Korzystając z tych samych danych, powtarzamy procedurę dopasowywania danych do Poissona GLM z trzema różnymi funkcjami łącza: ( a ) link, ( b ) łącze tożsamości i ( c ) łącze pierwiastka kwadratowego. Poniższa tabela pokazuje szacunki współczynników po konwersji z powrotem do parametryzacji dziennika. (Tak więc druga kolumna pokazuje a trzecia pokazuje przy użyciu surowego z każdego pasowania). Ponownie, szacunki są identyczne.log ( β ) log ( β 2 ) $\log$$\log(\hat \beta)$$\log(\hat \beta^2)$$\hat\beta$

> coefs.po
         log       id     sqrt
XX1 4.234107 4.234107 4.234107
XX2 4.790820 4.790820 4.790820
XX3 4.841033 4.841033 4.841033

Zastrzeżenie w nagłówku odnosi się po prostu do faktu, że surowe oszacowania będą się różnić w zależności od funkcji łącza, ale implikowane oszacowania parametru średniego nie będą.

Kod R.

# Warning! This code is a bit simplified for compactness.
library(MASS)
n <- 5
m <- 3
set.seed(17)
b <- exp(5+rnorm(m))
k <- 6

# Random negbin data; orthogonal factors
y <- rnbinom(m*n, size=k, mu=rep(b,each=n))
X <- factor(paste("X",rep(1:m,each=n),sep=""))

# Fit a bunch of GLMs with a log link
con <- glm.control(maxit=100)
mnb <- glm(y~X+0, family=negative.binomial(theta=2))
mpo <- glm(y~X+0, family="poisson")
mga <- glm(y~X+0, family=gaussian(link=log), start=rep(1,m), control=con)
miv <- glm(y~X+0, family=inverse.gaussian(link=log), start=rep(2,m), control=con)
mgm <- glm(y~X+0, family=Gamma(link=log), start=rep(1,m), control=con)    
coefs <- cbind(negbin=mnb$coef, poisson=mpo$coef, gaussian=mga$coef
                   invgauss=miv$coef, gamma=mgm$coef)

# Fit a bunch of Poisson GLMs with different links.
mpo.log  <- glm(y~X+0, family=poisson(link="log"))
mpo.id   <- glm(y~X+0, family=poisson(link="identity"))
mpo.sqrt <- glm(y~X+0, family=poisson(link="sqrt"))   
coefs.po <- cbind(log=mpo$coef, id=log(mpo.id$coef), sqrt=log(mpo.sqrt$coef^2))

Aby zobaczyć, co się tutaj dzieje, warto najpierw wykonać regresję bez przechwytywania, ponieważ przechwytywanie w regresji kategorycznej za pomocą tylko jednego predyktora jest bez znaczenia:

> rs1 = glm(breaks ~ tension-1, data=warpbreaks, family="poisson")
> rs2 = glm.nb(breaks ~ tension-1, data=warpbreaks)

Ponieważ Poisson i ujemne regresje dwumianowe określają log parametru średniego, to w przypadku regresji kategorialnej wykładnik współczynników da rzeczywisty parametr średni dla każdej kategorii:

>  exp(cbind("Poisson"=coef(rs1), "NB"=coef(rs2)))
          Poisson       NB
tensionL 36.38889 36.38889
tensionM 26.38889 26.38889
tensionH 21.66667 21.66667

Te parametry odpowiadają faktycznym średnim wartościom różnych kategorii:

> with(warpbreaks,tapply(breaks,tension,mean))
       L        M        H 
36.38889 26.38889 21.66667 

Tak więc dzieje się tak, że średni parametr w każdym przypadku, który maksymalizuje prawdopodobieństwo, jest równy średniej próbki dla każdej kategorii. $\lambda$

W przypadku rozkładu Poissona jasne jest, dlaczego tak się dzieje. Jest tylko jeden parametr do dopasowania, a zatem maksymalizacja ogólnego prawdopodobieństwa modelu z jednym predyktorem jakościowym jest równoważna niezależnemu znalezieniu która maksymalizuje prawdopodobieństwo obserwacji w każdej konkretnej kategorii. Estymator największego prawdopodobieństwa dla rozkładu Poissona jest po prostu średnią próbki, dlatego współczynniki regresji są dokładnie (logami) średnich próbek dla każdej kategorii.$\lambda$

W przypadku ujemnego dwumianowego nie jest to tak proste, ponieważ pasują do niego dwa parametry: i parametr kształtu . Co więcej, regresja pasuje do pojedynczego który obejmuje cały zestaw danych, więc w tej sytuacji regresja kategoryczna nie jest po prostu równoważna dopasowaniu całkowicie oddzielnego modelu dla każdej kategorii. Jednak po przeanalizowaniu funkcji wiarygodności możemy zauważyć, że dla każdej danej theta funkcja wiarygodności jest ponownie maksymalizowana poprzez ustawienie na próbkę średnią:$\lambda$$\theta$$\theta$$\lambda$

λ=ˉ

$$\begin{align} L(X,\lambda,\theta) &= \prod \left(\frac{\theta}{\lambda+\theta}\right)^\theta \frac{\Gamma(\theta + x_i)}{x_i!\Gamma(\theta)}\left(\frac{\lambda}{\lambda+\theta}\right)^{x_i}\\ \log L(X,\lambda,\theta) &= \sum\theta\left(\text{log}\theta-\text{log}(\lambda+\theta)\right) +x_i\left(\text{log}\lambda-\text{log}(\lambda+\theta)\right) +\log\left(\frac{\Gamma(\theta + x_i)}{x_i!\Gamma(\theta)}\right)\\ \frac{d}{d\lambda}\log L(X,\lambda,\theta) &= \sum \frac{x_i}{\lambda}-\frac{\theta+x_i}{\lambda+\theta} =n\left(\frac{\bar{x}}{\lambda}-\frac{\bar{x}+\theta}{\lambda+\theta}\right), \end{align}$$ więc maksimum jest osiągane, gdy .$\lambda=\bar{x}$

Powodem, dla którego nie otrzymujesz takich samych współczynników dla ciągłych danych, jest to, że w regresji ciągłej nie będzie już składową stałą funkcji zmiennych predykcyjnych, ale liniową. Maksymalizacja funkcji prawdopodobieństwa w tym przypadku nie ograniczy się do samodzielnego dopasowania wartości dla rozłącznych podzbiorów danych, ale będzie raczej nietrywialnym problemem, który jest rozwiązywany numerycznie i może przynieść różne wyniki dla różnych funkcji wiarygodności.$\text{log}(\lambda)$$\lambda$

Podobnie, jeśli masz wiele predyktorów jakościowych, pomimo faktu, że dopasowany model ostatecznie określi jako funkcję cząstkową stałej, ogólnie rzecz biorąc, nie będzie wystarczającej liczby stopni swobody, aby umożliwić niezależne określenie dla każdego stałego segmentu. Załóżmy na przykład, że masz predyktory z kategoriami. W takim przypadku Twój model ma stopni swobody, podczas gdy istnieje unikalnych różnych kombinacji kategorii, z których każda będzie miała własną dopasowaną wartość . Zakładając, że przecięcia tych kategorii nie są puste (a przynajmniej żeλ 2 5 10 5 ∗ 5 = 25 λ $\lambda$$\lambda$$2$$5$$10$$5*5=25$$\lambda$$11$ z nich są niepuste), problem maksymalizacji prawdopodobieństwa ponownie staje się nietrywialny i generalnie generuje różne wyniki dla Poissona w porównaniu z ujemnym dwumianowym lub innym rozkładem.