Virginian yliopiston kirjaston tutkimusaineistopalvelut + Sciences

Kun mallinnetaan lukumääriä (eli kokonaislukuja, jotka ovat suurempia kuin tai yhtä suuria kuin 0), aloitetaan usein Poissonin regressiosta. Tämä on yleistetty lineaarinen malli, jossa oletetaan, että vasteella on Poisson-jakauma, joka on riippuvainen ennustajien painotetusta summasta. Voimme esimerkiksi mallintaa NFL:n pelinrakentajien dokumentoitujen aivotärähdysten lukumäärää pelattujen otteluiden ja hyökkäyslinjan kokonaiskokemusvuosien funktiona. Poisson-regression yksi mahdollinen haittapuoli on kuitenkin se, että se ei välttämättä kuvaa tarkasti lukumäärien vaihtelua.

Poisson-jakaumaa parametrisoidaan \(\lambda\), joka sattuu olemaan sekä sen keskiarvo että varianssi. Vaikka se on kätevä muistaa, se ei useinkaan ole realistinen. Laskentajakaumalla on yleensä varianssi, joka ei ole yhtä suuri kuin sen keskiarvo. Kun näin tapahtuu tiedoissa, joiden oletamme (tai toivomme) olevan Poisson-jakaumaa, puhumme ali- tai ylihajonnasta riippuen siitä, onko varianssi pienempi vai suurempi kuin keskiarvo. Poisson-regression suorittaminen laskentatiedoille, joissa esiintyy tällaista käyttäytymistä, johtaa malliin, joka ei sovi hyvin.

Yksi lähestymistapa, joka ratkaisee tämän ongelman, on negatiivinen binomiregressio. Negatiivinen binomijakauma kuvaa Poisson-jakauman tavoin kokonaislukujen esiintymistodennäköisyyksiä, jotka ovat suurempia tai yhtä suuria kuin 0. Poisson-jakaumasta poiketen varianssi ja keskiarvo eivät ole ekvivalentteja. Tämä viittaa siihen, että se voi olla hyödyllinen approksimaatio sellaisten lukumäärien mallintamiseen, joiden vaihtelu poikkeaa keskiarvosta. Negatiivisen binomijakauman varianssi on keskiarvon funktio, ja sillä on lisäksi parametri k, jota kutsutaan hajontaparametriksi. Sanotaan, että laskentamme on satunnaismuuttuja Y negatiivisesta binomijakaumasta, jolloin Y:n varianssi on

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

Kun hajontaparametri kasvaa yhä suuremmaksi, varianssi konvergoituu samaan arvoon kuin keskiarvo, ja negatiivinen binomijakauma muuttuu Poisson-jakaumaksi.

Negatiivisen binomijakauman havainnollistamiseksi työskennellään muutamalla aineistolla Alan Agrestin kirjasta Categorical Data Analysis (2002). Tiedot on esitetty taulukossa 13.6 kohdassa 13.4.3. Aineisto on peräisin 1308 henkilön kyselytutkimuksesta, jossa heiltä kysyttiin, kuinka monta henkirikoksen uhria he tuntevat. Muuttujat ovat resp, vastaajan tuntemien uhrien lukumäärä, ja race, vastaajan rotu (musta tai valkoinen). Selittääkö rotu osaltaan sitä, kuinka monta henkirikoksen uhria henkilö tuntee? Aineisto on ensin syötettävä R:ään:

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

Ennen mallintamista tutkitaan aineistoa. Ensin huomaamme, että suurin osa vastaajista on valkoihoisia:

> table(race)racewhite black 1149 159 

Mustien keskiluku on suurempi kuin valkoihoisten:

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

Kunkin rodun osalta otoksen varianssi on noin kaksinkertainen keskiarvoon verrattuna. Näyttää siltä, että kyseessä on ylidispersio.

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

Viimeiseksi tarkastellaan lukumäärien jakaumaa roduittain.

> 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ääsemme mallin sovittamiseen. Kokeilemme ensin Poissonin regressiota glm()-funktiolla ja näytämme osan yhteenvetotuloksesta.

> # 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 ***

Rotu on erittäin merkitsevä. Näyttää siltä, että mustat tuntevat paljon todennäköisemmin jonkun, joka on joutunut henkirikoksen uhriksi. Mutta mitä tarkoittaa kerroin 1,73? Tässä yksinkertaisessa mallissa, jossa on yksi dikotominen ennustaja, se on log-odotettujen lukumäärien ero. Jos potensoimme kertoimen, saamme otoskeskiarvojen suhteen:

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

Tosiasiassa, jos teemme ennusteen tällä mallilla ja potensoimme tulokset, saamme otoskeskiarvot:

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

Tämä kertoo, että valkoihoisten tuntemien uhrien lukumäärä jakaantuu Poissonin kaltaisesti siten, että keskiarvo ja varianssi ovat yhtä suuret kuin 0.09, kun taas mustien tunnettujen uhrien määrä jakaantuu Poissonin tavoin keskiarvon ja varianssin ollessa 0,52. Tiedämme jo eksploratiivisesta analyysistämme, että havaitut varianssit olivat paljon suurempia, joten meidän ei pitäisi olla kovin tyytyväisiä mallin estimoimiin variansseihin. Jos tarkastelemme sovitettuja lukumääriä, näemme vieläkin enemmän todisteita sovituksen puutteesta:

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

Tallensimme edellä ensin ennustetut keskiarvot objektiin nimeltä fmeans. Sen jälkeen generoimme sovitettuja lukumääriä käyttämällä dpois-funktiota yhdessä estimoitujen keskiarvojen kanssa ennustamaan todennäköisyyttä saada 0-6. Sitten kerroimme nämä todennäköisyydet vastaajien lukumäärällä saadaksemme sovitetut lukumäärät. Lopuksi yhdistimme kaiken datakehykseen, jotta havaittuja ja sovitettuja arvoja oli helppo verrata. Kaksi dramaattisempaa huomiota herättävää asiaa on se, että 0-lukuja ali- ja 1-lukuja ylisovitetaan.

Voimme käyttää juurikuvausta visualisoidaksemme lukumäärien regressiomallin sovitusta. Countreg-paketin rootogram()-funktio tekee tästä helppoa.

> countreg::rootogram(pGLM)

Punainen kaareva viiva on teoreettinen Poissonin sovitus. Punaisen viivan jokaisesta pisteestä “roikkuu” palkki, jonka korkeus edustaa odotettujen ja havaittujen lukumäärien välistä eroa. Alle 0:n roikkuva palkki osoittaa vajaasovitusta. Pylväs, joka roikkuu yli 0:n, osoittaa ylisovitusta. Lukumäärät on muunnettu neliöjuurimuunnoksella, jotta pienemmät lukumäärät eivät peittyisi ja hukkuisi suurempiin lukumääriin. Näemme paljon alisovitusta lukumäärille 2 ja sitä suuremmille lukumäärille ja massiivista ylisovitusta lukumäärälle 1.

Kokeillaan nyt negatiivisen binomialueen mallin sovittamista. Huomasimme, että laskentojen vaihtelu oli suurempaa molemmilla roduilla. Näyttäisi siltä, että negatiivinen binomijakauma lähestyisi paremmin laskentojen jakaumaa.

Negatiivisen binomimallin sovittamiseksi R:ssä käännymme MASS-paketin (R:n mukana asennettu paketti) glm.nb()-funktion puoleen. Näytämme jälleen vain osan yhteenvetotulosteesta:

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

Huomaa ensin, että kertoimet ovat samat kuin aiemmin. Jälleen kerran voimme potensoida rotukertoimen saadaksemme otoskeskiarvojen suhteen ja tehdä ennusteita saadaksemme alkuperäiset otoskeskiarvot.

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

Mutta huomaa, että rotukertoimen keskivirhe on suurempi, mikä viittaa suurempaan epävarmuuteen estimaatissamme (0,24 vs. 0,15). Tämä on järkevää, kun otetaan huomioon laskuissamme havaittu vaihtelu. Huomaa myös Thetan estimaatti. Se on hajontaparametrimme. Voimme käyttää sitä suoraan malliobjektistamme seuraavasti:

> nbGLM$theta 0.2023119

Ja voimme käyttää sitä saadaksemme estimoidut varianssit laskennoille:

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

Nämä ovat paljon lähempänä havaittuja variansseja kuin Poisson-mallin antamat varianssit.

Jälleen kerran visualisoimme sovitusta juurimallin avulla:

> countreg::rootogram(nbGLM)

Tämä näyttää paljon paremmalta kuin Poisson-mallin juurimallin juurikuva. Laskentojen 1-3 kohdalla on pientä ali- tai ylisovittumista, mutta muuten se näyttää melko hyvältä.

Voidaksemme saada lisää tietoa negatiivisesta binomimallistamme, käytetään sen parametreja datan simulointiin ja verrataan simuloitua dataa havaittuun dataan. Alla lataamme magrittr-paketin, jotta pääsemme käsiksi %>%-operaattoriin, jonka avulla voimme “ketjuttaa” funktioita.

Aluksi taulukoimme lukumäärät ja luomme pylväsdiagrammin valkoisille ja mustille osallistujille. Sitten käytämme malliparametreja simuloidaksemme tietoja negatiivisesta binomijakaumasta. Nämä kaksi parametria ovat mu ja size (eli hajontaparametri). Huomaa, että käytämme coef()-funktiota poimimaan kullekin rodulle sopivat kertoimet. Valkoisille se on vain leikkauskerroin, mustille se on leikkauskerroin ja kaltevuus (summaamme ne siis yhteen). Tämän jälkeen eksponentit muunnetaan logaritmisesta asteikosta alkuperäiseen asteikkoon. Simuloimme saman määrän havaintoja kuin mitä meillä on alkuperäisessä aineistossamme.

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)

Simuloitu aineisto on hyvin samankaltainen kuin havaittu aineisto, mikä taas antaa meille luottamusta negatiivisen binomiregression valitsemiseen tämän aineiston mallintamiseen. Nämä kuvaajat osoittavat myös mallimme ehdollisen luonteen. Lukumäärien negatiivinen binomijakauma riippuu rodusta tai on riippuvainen siitä. Jokaisella rodulla on erilainen keskiarvo mutta yhteinen hajontaparametri.

• 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

Jos sinulla on tätä artikkelia koskevia kysymyksiä tai selvennyksiä, ota yhteyttä UVA:n kirjaston StatLab:iin: [email protected]

Katsele koko UVA:n kirjaston StatLab-artikkelikokoelmaa.

Clay Ford
Statistical Research Consultant
Virginiassa sijaitsevan yliopiston kirjasto
.