Pokud jde o modelování počtů (tj. celých čísel větších nebo rovných 0), často začínáme s Poissonovou regresí. Jedná se o zobecněný lineární model, kde se předpokládá, že odpověď má Poissonovo rozdělení podmíněné váženým součtem prediktorů. Například můžeme modelovat počet zdokumentovaných otřesů mozku u quarterbacků NFL jako funkci počtu odehraných zákroků a celkového počtu let zkušeností jeho útočné řady. Jednou z možných nevýhod Poissonovy regrese však je, že nemusí přesně popisovat variabilitu počtů.

Poissonovo rozdělení je parametrizováno pomocí \(\lambda\), což je shodou okolností jeho střední hodnota i rozptyl. I když je to pohodlné na zapamatování, není to často reálné. Rozdělení počtů bude mít obvykle rozptyl, který se nerovná jeho střední hodnotě. Když se to stane u dat, o kterých předpokládáme (nebo doufáme), že jsou Poissonovým rozdělením, říkáme, že máme pod- nebo nadměrný rozptyl, podle toho, zda je rozptyl menší nebo větší než průměr. Provedení Poissonovy regrese na počítaných datech, která vykazují toto chování, vede k modelu, který dobře nesedí.

Jeden z přístupů, který tento problém řeší, je záporná binomická regrese. Záporné binomické rozdělení, stejně jako Poissonovo rozdělení, popisuje pravděpodobnosti výskytu celých čísel větších nebo rovných 0. Na rozdíl od Poissonova rozdělení nejsou rozptyl a střední hodnota ekvivalentní. To naznačuje, že by mohlo sloužit jako užitečná aproximace pro modelování počtů s variabilitou odlišnou od jeho střední hodnoty. Rozptyl záporného binomického rozdělení je funkcí jeho střední hodnoty a má další parametr k, který se nazývá parametr rozptylu. Řekněme, že náš počet je náhodná veličina Y ze záporného binomického rozdělení, pak rozptyl Y je

$$ var(Y) = \mu + \mu^{2}/k $$

Když se parametr rozptylu zvětšuje, rozptyl konverguje ke stejné hodnotě jako střední hodnota a záporný binomický počet se mění na Poissonovo rozdělení.

Pro ilustraci záporného binomického rozdělení pracujme s daty z knihy Analýza kategoriálních dat od Alana Agrestiho (2002). Data jsou uvedena v tabulce 13.6 v kapitole 13.4.3. Data pocházejí z průzkumu, kterého se zúčastnilo 1308 lidí a v němž byli dotázáni, kolik znají obětí vražd. Proměnnými jsou resp, počet obětí, které respondent zná, a race, rasa respondenta (černoch nebo běloch). Pomáhá rasa vysvětlit, kolik obětí vražd daná osoba zná? Data je třeba nejprve zadat do programu R:

> # Table 13.6> # Agresti, p. 561> black <- c(119,16,12,7,3,2,0)> white <- c(1070,60,14,4,0,0,1)> resp <- c(rep(0:6,times=black), rep(0:6,times=white))> race <- factor(c(rep("black", sum(black)), rep("white", sum(white))),+ levels = c("white","black"))> victim <- data.frame(resp, race)

Než se pustíme do modelování, prozkoumejme data. Nejprve si všimneme, že většina respondentů jsou běloši:

> table(race)racewhite black 1149 159 

Černoši mají vyšší průměrný počet než běloši:

> with(victim, tapply(resp, race, mean)) white black 0.09225413 0.52201258 

Pro každou rasu je výběrový rozptyl zhruba dvakrát větší než průměr. Zdá se, že máme nadměrný rozptyl:

> with(victim, tapply(resp, race, var)) white black 0.1552448 1.1498288 

Nakonec se podíváme na rozdělení počtů podle ras.

> table(resp, race) raceresp white black 0 1070 119 1 60 16 2 14 12 3 4 7 4 0 3 5 0 2 6 1 0

Přejdeme k fitování modelu. Nejprve vyzkoušíme Poissonovu regresi pomocí funkce glm() a ukážeme část souhrnného výstupu.

> # Poisson model> pGLM <- glm(resp ~ race, data=victim, family = poisson)> summary(pGLM)Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.38321 0.09713 -24.54 <2e-16 ***raceblack 1.73314 0.14657 11.82 <2e-16 ***

Rasa je velmi významná. Ukazuje se, že černoši mnohem častěji znají někoho, kdo se stal obětí vraždy. Co však znamená koeficient 1,73? V tomto jednoduchém modelu s jedním dichotomickým prediktorem je to rozdíl v logaritmech očekávaných počtů. Pokud koeficient exponenciálně vynásobíme, dostaneme poměr výběrových průměrů:

> exp(coef(pGLM))raceblack 5.658419 > # same thing> mean(victim$resp)/mean(victim$resp) 5.658419

Ve skutečnosti, pokud provedeme predikci s tímto modelem a výsledky exponenciálně vynásobíme, dostaneme výběrové průměry:

> exp(predict(pGLM, newdata = data.frame(race=c("white","black")))) 1 2 0.09225413 0.52201258 > # same thing> with(victim, tapply(resp, race, mean)) white black 0.09225413 0.52201258 

To říká, že počet známých obětí pro bělochy je rozdělen jako Poissonovo rozdělení se střední hodnotou a rozptylem rovným nule.09, zatímco počet známých obětí pro černochy je rozdělen jako Poissonovo rozdělení se střední hodnotou a rozptylem rovným 0,52. Z naší průzkumné analýzy již víme, že pozorované rozptyly byly mnohem větší, takže bychom neměli být příliš spokojeni s rozptyly odhadnutými modelem. Pokud prozkoumáme fitované počty, uvidíme ještě více důkazů o nedostatečném přizpůsobení:

> # fitted counts for Poisson GLM:> fmeans <- exp(predict(pGLM, newdata = data.frame(race = c("white","black"))))> fmeans 1 2 0.09225413 0.52201258 > fittedW <- dpois(0:6,lambda = fmeans) * sum(victim$race=="white") > fittedB <- dpois(0:6,lambda = fmeans) * sum(victim$race=="black") > data.frame(Response=0:6,BlackObs=black, BlackFit=round(fittedB,1), + WhiteObs=white, WhiteFit=round(fittedW,1)) Response BlackObs BlackFit WhiteObs WhiteFit1 0 119 94.3 1070 1047.72 1 16 49.2 60 96.73 2 12 12.9 14 4.54 3 7 2.2 4 0.15 4 3 0.3 0 0.06 5 2 0.0 0 0.07 6 0 0.0 1 0.0

Výše jsme nejprve uložili předpovězené průměry do objektu s názvem fmeans. Poté jsme pomocí funkce dpois spolu s odhadnutými průměry vytvořili fitované počty, abychom předpověděli pravděpodobnost získání hodnot 0 až 6. Poté jsme tyto pravděpodobnosti vynásobili počtem respondentů, abychom získali fitované počty. Nakonec jsme vše spojili do datového rámce, abychom mohli snadno porovnat pozorované a fitované hodnoty. Dvě z nejdramatičtějších věcí, kterých je třeba si všimnout, je to, že jsme nedostatečně přizpůsobili počty 0 a nadměrně přizpůsobili počty 1.

K vizualizaci fitování regresního modelu počtů můžeme použít rootogram. To usnadňuje funkce rootogram() v balíčku countreg.

> countreg::rootogram(pGLM)

Červená zakřivená čára je teoretický Poissonův fit. Z každého bodu na červené přímce “visí” sloupec, jehož výška představuje rozdíl mezi očekávanými a pozorovanými počty. Sloupec visící pod 0 znamená nedostatečné přizpůsobení. Sloupec visící nad 0 znamená nadměrné přizpůsobení. Počty byly transformovány pomocí odmocninové transformace, aby se zabránilo tomu, že menší počty budou zastíněny a přehlušeny většími počty. Vidíme velkou míru podhodnocení pro počty 2 a vyšší a masivní nadhodnocení pro počet 1.

Nyní zkusíme fitovat negativní binomický model. Všimli jsme si, že variabilita počtů byla větší pro obě rasy. Zdá se, že záporné binomické rozdělení by lépe aproximovalo rozdělení počtů.

Pro fitování záporného binomického modelu v R se obrátíme na funkci glm.nb() v balíčku MASS (balíček, který se instaluje s R). Opět si ukážeme pouze část souhrnného výstupu:

> # Fit a negative binomial model> library(MASS)> nbGLM <- glm.nb(resp ~ race, data=victim)> summary(nbGLM)Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -2.3832 0.1172 -20.335 < 2e-16 ***raceblack 1.7331 0.2385 7.268 3.66e-13 *** Theta: 0.2023 

Nejprve si všimněte, že koeficienty jsou stejné jako dříve. Opět můžeme exponovat koeficient závodu, abychom získali poměr výběrových průměrů, a provést predikci, abychom získali původní výběrové průměry.

> # fitted counts for Negative Binomial GLM:> fmeans <- exp(predict(pGLM, newdata = data.frame(race = c("white","black"))))> fmeans # same as pGLM 1 2 0.09225413 0.52201258 

Všimněte si však, že standardní chyba pro koeficient rasy je větší, což naznačuje větší nejistotu našeho odhadu (0,24 oproti 0,15). To dává smysl vzhledem k pozorované variabilitě našich počtů. Všimněte si také odhadu Theta. To je náš parametr rozptylu. Můžeme k němu přistupovat přímo z našeho objektu modelu takto:

> nbGLM$theta 0.2023119

A můžeme jej použít k získání odhadovaných rozptylů pro počty:

> fmeans + fmeans^2 * (1/nbGLM$theta) 1 2 0.134322 1.868928 

Ty jsou mnohem blíže pozorovaným rozptylům než rozptyly dané Poissonovým modelem.

Znovu vizualizujeme shodu pomocí rootogramu:

> countreg::rootogram(nbGLM)

Ten vypadá mnohem lépe než rootogram Poissonova modelu. U počtů 1 až 3 dochází k mírnému podhodnocení/přehodnocení, ale jinak to vypadá docela dobře.

Abychom získali další poznatky k našemu negativnímu binomickému modelu, použijeme jeho parametry k simulaci dat a porovnáme simulovaná data s pozorovanými daty. Níže načteme balíček magrittr pro přístup k operátoru %>%, který nám umožňuje “řetězit” funkce.

Nejprve sestavíme tabulku počtů a vytvoříme sloupcový graf pro bílé, resp. černé účastníky. Poté použijeme parametry modelu k simulaci dat ze záporného binomického rozdělení. Těmito dvěma parametry jsou mu a velikost (tj. parametr rozptylu). Všimněte si, že k extrakci příslušných koeficientů pro jednotlivé rasy používáme funkci coef(). Pro bělochy je to pouze intercept, pro černochy je to intercept a sklon (tedy je sečteme). Poté provedeme exponenciální převod z logaritmického měřítka na původní měřítko. Simulujeme stejný počet pozorování, jaký máme v původních datech.

library(magrittr) # for %>% operatorop <- par(mfrow=c(1,2))set.seed(1)victim$resp %>% `)) %>% table() %>% barplot(main = "Simulated White")victim$resp %>% `[`(victim$race=="black") %>% table() %>% barplot(main = "Observed Black")rnbinom(n = 159, size = nbGLM$theta, mu = exp(sum(coef(nbGLM)))) %>% table() %>% barplot(main = "Simulated Black")par(op)

Simulovaná data jsou velmi podobná pozorovaným datům, což nám opět dává důvěru ve volbu negativní binomické regrese pro modelování těchto dat. Tyto grafy také ukazují podmíněnost našeho modelu. Záporné binomické rozdělení počtů závisí nebo je podmíněno rasou. Každá rasa má jiný průměr, ale společný parametr rozptylu.

• Agresti, Alan (2002), Categorical Data Analysis, Wiley.
• Kleiber, Christian & Zeileis, Achim (2016): Visualizing Count Data Regressions Using Rootograms, The American Statistician, DOI: 10.1080/00031305.2016.1173590

V případě dotazů nebo vysvětlení týkajících se tohoto článku kontaktujte StatLab knihovny UVA: [email protected]

Podívejte se na celou kolekci článků StatLab knihovny UVA.

Clay Ford
Statistical Research Consultant
University of Virginia Library
.