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

38

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))

wprowadź opis zdjęcia tutaj

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))

wprowadź opis zdjęcia tutaj

(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))

wprowadź opis zdjęcia tutaj

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.883b=0.881α1/θ

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

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))

wprowadź opis zdjęcia tutaj

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))

wprowadź opis zdjęcia tutaj

półprzejście
źródło
6
(+1) Ładne pytanie. Spróbuję coś dla ciebie napisać za kilka godzin. W międzyczasie możesz spróbować zobaczyć, co dzieje się z wieloma predyktorami jakościowymi, które nie są ortogonalne (wskazówka!).
kardynał
1
Intryga! Na pewno tego spróbuję. I dziękuję, z niecierpliwością czekam na twoją odpowiedź.
półfinał
Przepraszam @ half-pass; Nie zapomniałem o tym i spróbuję coś wstać w ciągu jednego dnia. (Mam połowę odpowiedzi ułożoną razem, ale pociągnęły mnie inne zadania.) Mam nadzieję, że nagroda przyciągnie także inne odpowiedzi. Twoje zdrowie. :-)
kardynał
Nie martw się, @cardinal! Wiem, że ty i wszyscy inni niesamowici guru mają tutaj życie poza CV :)
półfinał

Odpowiedzi:

29

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!μ^j=y¯jj

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ć

logf(y;θ,ν)=θyb(θ)ν+a(y,ν).

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.θνν

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 .θy=b(θ^)νy¯=b(θ^)

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

θElogf(Y;θ,ν)=Eθlogf(Y;θ,ν)=0.
b(θ)=EY=μ

Ponieważ szacunki maksymalnego prawdopodobieństwa są niezmienne podczas transformacji, oznacza to, że dla tej rodziny gęstości.y¯=μ^

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.μiμi=g1(xiTβ)gxijμi=g(βj)β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 i

i=1n(yiμi)xijσi2μiλi=0,
j=1,,pλi=xiTββνμi=g(λi)=g(xiTβ)σi2

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 ) βloglog(β^)log(β^2)β^

> 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))
kardynał
źródło
+1, jasne i wyczerpujące, lepsza odpowiedź niż to, co dałem.
mpr
@mpr, wolałbym uważać to za komplementarne do twojego. Byłem bardzo zadowolony, kiedy zobaczyłem twój post; jasno i zwięźle opisałeś (i pokazałeś), co się dzieje. Moje posty czasem czytają trochę; Obawiam się, że to jeden z nich.
kardynał
Oboje jesteście niesamowici. Dziękuję bardzo za wyjaśnienie mi tego z taką jasnością i surowością. Muszę teraz poświęcić trochę więcej czasu na trawienie :)
półfinał
@cardinal Skąd masz 2 w „ family=negative.binomial(theta=2)”?
timothy.s.lau,
23

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. λ

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.λ

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ą:λθθλ

λ=ˉx

L(X,λ,θ)=(θλ+θ)θΓ(θ+xi)xi!Γ(θ)(λλ+θ)xilogL(X,λ,θ)=θ(logθlog(λ+θ))+xi(logλlog(λ+θ))+log(Γ(θ+xi)xi!Γ(θ))ddλlogL(X,λ,θ)=xiλθ+xiλ+θ=n(x¯λx¯+θλ+θ),
więc maksimum jest osiągane, gdy .λ=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.λlog(λ)λ

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 λ 11λλ251055=25λ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.

mpr
źródło
5
(+1) Dobra odpowiedź. Jedną rzeczą, którą myślę, że można tu wyjaśnić bardziej wyraźnie, jest to, że tak naprawdę nie ma to nic wspólnego z jakimkolwiek związkiem między Poissonem a ujemnym dwumianem, ale raczej z bardziej podstawowymi faktami na temat dopasowania GLM z maksymalnym prawdopodobieństwem.
kardynał
Słuszna uwaga. Jedyną prawdziwą rzeczą, z którą ma związek Poisson i ujemny dwumian, jest to, że są one określone w logarytmie parametru średniego. Na przykład, jeśli wykonałeś zwykłą regresję najmniejszych kwadratów, współczynniki byłyby nominalnie różne, ponieważ wówczas parametr byłby rzeczywistą średnią, a nie logarytmiczną.
mpr
2
To prawda, ale myślę, że wykracza poza to: dopasuj dowolny GLM za pomocą dowolnej funkcji łącza (według ML i z bardzo małymi zastrzeżeniami), a ponieważ dopasowane środki będą pasować do średnich próbek, wówczas oszacowania parametrów są identyczne do nieliniowych transformacja między różnymi linkami. Konkretny model błędu jest nieistotny poza faktem, że pochodzi on z wykładniczej rodziny dyspersyjnej.
kardynał
To dobra uwaga, której nie omówiłem. Podszedłem do tego pytania z bardziej ogólnego punktu widzenia oszacowania ML, a nie GLM. Istnieje wiele naturalnie występujących modeli, w których ML produkuje dopasowane środki, które różnią się od średnich próbek (np. Lognormalne). Ale w przypadku GLM twoja obserwacja prowadzi do bardziej zwięzłego i bardziej ogólnego wyjaśnienia.
mpr
2
@ half-pass: Dopasuj wszystkie swoje ortogonalne modele jakościowe y~X+0i spróbuj ponownie. :-)
kardynał