Website do GR

4.A Modelos para Dados de Contagem Superdispersos

Usamos dados de Long (1990) sobre o número de publicações produzidas por bioquímicos Ph.D. para ilustrar a aplicação de Poisson, Poisson superdispersos, binômio negativo e modelos de Poisson com inflação zero.

As variáveis no conjunto de dados são

phd

Variável	Descrição
art	artigos nos últimos três anos de Ph.D.
fem	codificado para mulheres
mar	codificado se casadas
kid5	número de crianças menores de seis anos
prestígio de Ph.D. programa
ment	artigos por mentor nos últimos três anos

Estes dados também foram analisados por Long e Freese (2001), e estão disponíveis no site da Stata:

> library(foreign)> ab names(ab) "art" "fem" "mar" "kid5" "phd" "ment"> r c(mean=r, var=r, ratio=r/r) mean var ratio 1.692896 3.709742 2.191358 

O número médio de artigos é 1,69 e a variância é 3,71, um pouco mais que o dobro da média. Os dados estão superdispersos, mas é claro que ainda não consideramos nenhum covariante.

Um Modelo Poisson

Deixe-nos ajustar o modelo usado por Long e Freese(2001), um modelo aditivo simples usando todos os cinco preditores.

> mp summary(mp)Call:glm(formula = art ~ fem + mar + kid5 + phd + ment, family = poisson, data = ab)Deviance Residuals: Min 1Q Median 3Q Max -3.5672 -1.5398 -0.3660 0.5722 5.4467 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.304617 0.102981 2.958 0.0031 ** femWomen -0.224594 0.054613 -4.112 3.92e-05 ***marMarried 0.155243 0.061374 2.529 0.0114 * kid5 -0.184883 0.040127 -4.607 4.08e-06 ***phd 0.012823 0.026397 0.486 0.6271 ment 0.025543 0.002006 12.733 
We see that the model obviously doesn't fit the data. The five-percent critical value for a chi-squared with 909 d.f. is
>qchisq(0.95, df.residual(mp)) 980.2518> deviance(mp) 1634.371> pr sum(pr^2) 1662.547
e o desvio e o qui-quadrado de Pearson estão ambos no 1600s.
 4.A.1. Variação Extra-Poisson
Agora assumimos que a variância é proporcional ao invés de igual à média, e estimamos o parâmetro da escala φ dividindo o qui-quadrado de Pearson pelo seu d.f.:
> phi round(c(phi,sqrt(phi)),4) 1.8290 1.3524
Vemos que a variância é cerca de 83% maior do que a média. Isto significa que devemos ajustar os erros padrão multiplicando por 1,35, a raiz quadrada de 1,83.
R pode fazer este cálculo para nós se usarmos a quasipoissonfamily:
> mq summary(mq)Call:glm(formula = art ~ fem + mar + kid5 + phd + ment, family = quasipoisson, data = ab)Deviance Residuals: Min 1Q Median 3Q Max -3.5672 -1.5398 -0.3660 0.5722 5.4467 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.304617 0.139273 2.187 0.028983 * femWomen -0.224594 0.073860 -3.041 0.002427 ** marMarried 0.155243 0.083003 1.870 0.061759 . kid5 -0.184883 0.054268 -3.407 0.000686 ***phd 0.012823 0.035700 0.359 0.719544 ment 0.025543 0.002713 9.415 
The estimates are exactly the same as before, but the standard errorsare about 35% larger. We can verify this fact easily. First we writea useful function to extract standard errors and then use it onour fits:
> se round(data.frame(p=coef(mp), q=coef(mq),+ se.p=se(mp), se.q=se(mq), ratio=se(mq)/se(mp)), 4) p q se.p se.q ratio(Intercept) 0.3046 0.3046 0.1030 0.1393 1.3524femWomen -0.2246 -0.2246 0.0546 0.0739 1.3524marMarried 0.1552 0.1552 0.0614 0.0830 1.3524kid5 -0.1849 -0.1849 0.0401 0.0543 1.3524phd 0.0128 0.0128 0.0264 0.0357 1.3524ment 0.0255 0.0255 0.0020 0.0027 1.3524
Uma abordagem alternativa é encaixar um modelo de Poisson e usar therobust ou estimador de sanduíche dos erros padrão. Isto normalmente dá resultados muito semelhantes aos do modelo Poisson sobre-disperso.
 4.A.2 Regressão Binomial Negativa
Agora encaixamos um modelo binomial negativo com os mesmos preditores.Para isso precisamos da função glm.nb() no pacoteMASS.
> library(MASS)> mnb summary(mnb)Call:glm.nb(formula = art ~ fem + mar + kid5 + phd + ment, data = ab, init.theta = 2.264387695, link = log)Deviance Residuals: Min 1Q Median 3Q Max -2.1678 -1.3617 -0.2806 0.4476 3.4524 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.256144 0.137348 1.865 0.062191 . femWomen -0.216418 0.072636 -2.979 0.002887 ** marMarried 0.150489 0.082097 1.833 0.066791 . kid5 -0.176415 0.052813 -3.340 0.000837 ***phd 0.015271 0.035873 0.426 0.670326 ment 0.029082 0.003214 9.048 1/mnb$theta 0.4416205> -2*(logLik(mp)-logLik(mnb))'log Lik.' 180.196 (df=6)
R's theta é a precisão do efeito multiplicativador aleatório, e corresponde a 1/σ2 em thenotes. A estimativa corresponde a uma variância estimada de 0,44 e é altamente significativa.
Para testar a significância deste parâmetro você pode pensar em computar duas vezes a diferença em log-likelihoods entre este modelo e o modelo de Poisson, 180,2, e tratá-lo como um qui-quadrado com um d.f. As assimptóticas usuais não se aplicam, no entanto, porque a hipótese nula está num limite do espaço do parâmetro.
Há algum trabalho mostrando que uma melhor aproximação é tratar a estatística como uma mistura 50:50 de zero e um qui-quadrado com um d.f. Alternativamente, tratar a estatística como um qui-quadrado com um d.f. dá um teste conservador. De qualquer forma, temos evidência esmagadora de dispersão excessiva.
Para testar hipóteses sobre os coeficientes de regressão, podemos usar testes de Wald ou testes de razão de verossimilhança, que são possíveis porque fizemos suposições distributivas completas.
Há também uma família negative.binomial para glm e isto pode ser usado desde que os parâmetrostheta sejam dados. (Isto é baseado no resultado que o binômio negativo está na família glm para a variância fixa). Aqui está uma verificação rápida:
> mnbg all(abs(coef(mnb)==coef(mnbg)) 
The standard errors would differ, however, because  allows for the fact that was estimated, whereas  does not.
Unobserved Heterogeneity
R has a function  to compute the density of a gamma distribution with given shape and scale (or its reciprocal the rate). In particular, when the random effect hasvariance  the density is. This can beused to draw the density.
> v = 1/mnb$theta> gd = data.frame(x = seq(0.0, 3, 0.1), g = dgamma(x, shape = 1/v, scale = v))> ggplot(gd, aes(x, g)) + geom_line()> ggsave("gamdenr.png", width=500/72, height=400/72, dpi=72)

 Também podemos computar quantis usandoqgamma(p, shape, rate = 1, scale = 1/rate, lower.tail = TRUE). Em particular, quando o efeito aleatório tem variância v os quartis são qgamma((1:3)/4, shap e= 1/v, scale = v).
> qgamma((1:3)/4, shape = 1/v, scale = v) 0.5114167 0.8572697 1.3347651
Biochemista no Q1 da distribuição de heterogeneidade não observada publica 49% menos artigos do que o esperado de suas características observadas, enquanto os da mediana publicam 14% menos e os do Q3 publicam 33% mais do que o esperado.
Estimativas de comparação e erros-padrão
Deixe-nos comparar estimativas de parâmetros e erros-padrão sob o modelo de Poisson, Poisson sobre-disperso e modelos binomiais negativos:
> round(data.frame(+ p=coef(mp),q=coef(mq),nb=coef(mnb),+ se.p=se(mp),se.q=se(mq),se.nb=se(mnb)),4) p q nb se.p se.q se.nb(Intercept) 0.3046 0.3046 0.2561 0.1030 0.1393 0.1373femWomen -0.2246 -0.2246 -0.2164 0.0546 0.0739 0.0726marMarried 0.1552 0.1552 0.1505 0.0614 0.0830 0.0821kid5 -0.1849 -0.1849 -0.1764 0.0401 0.0543 0.0528phd 0.0128 0.0128 0.0153 0.0264 0.0357 0.0359ment 0.0255 0.0255 0.0291 0.0020 0.0027 0.0032
As estimativas binomiais negativas não são muito diferentes daquelas baseadas no modelo de Poisson, e ambos os conjuntos levariam às mesmas conclusões. 
Olhando para os erros padrão, vemos que ambas as abordagens de sobre-dispersão levam a erros padrão estimados muito semelhantes, e que a regressão normal de Poisson subestima os erros padrão.
Bom do ajuste
Podemos avaliar a bondade do ajuste do modelo binomial negativo usando o desvio
> deviance(mnbg) 1004.281
O modelo binomial negativo encaixa melhor que o Poisson, mas ainda tem um desvio acima do valor crítico de 5%980,25,
A função de desvio
Os modelos sobre-dispersão de Poisson e binomial negativo têm diferentes funções de desvio. Uma maneira de verificar qual pode ser a mais apropriada é criar grupos baseados no preditor linear, calcular a média e variância de cada grupo e, finalmente, plotar a relação média-variância.
Aqui estão os grupos baseados no preditor linear binomial negativo, criados usando cut com quebras no percentil 5(5)95 produzem vinte grupos de tamanho aproximadamente igual. Note que predict() computa um preditor linear, a menos que especificado de outra forma. Para predizer na escala da respostaadd a opção type="response". 
> xb g m v png("c4afig1.png",width=500,height=400)> plot(m, v, xlab="Mean", ylab="Variance", + main="Mean-Variance Relationship")> mtext("Articles Published by Ph.D. Biochemists",padj=-0.5)> x lines(x, x*phi, lty="dashed")> lines(x, x*(1+x/mnb$theta))> legend("topleft", lty=c("dashed","solid"), + legend=c("Q. Poisson","Neg. Binom."), inset=0.05)> dev.off();null device 1 

O gráfico traça a média versus a variância e sobrepõe as curvas correspondentes ao modelo de Poisson sobre-disperso, onde a variância é φμ, e o modelo de binomialmodel negativo, onde a variância é μ(1+μσ2).
A função de variância de Poisson faz um bom trabalho para a maior parte dos dados, mas não captura as altas variâncias dos estudiosos mais produtivos. A função de variância binomial negativa não é muito diferente mas, sendo um quadrático, pode subir mais rápido e faz um trabalho melhor na parte alta. Concluímos que o modelo binomial negativo fornece uma melhor descrição dos dados do que o modelo de Poisson sobre-disperso.
 4.A.3 Modelos de Poisson Inflacionados a Zero
 Uma ocorrência frequente com dados de contagem é um excesso de zeros em relação ao que é esperado sob um modelo de Poisson. Isto é na verdade um problema com os nossos dados: 
> zobs zpoi c(obs=mean(zobs), poi=mean(zpoi)) obs poi 0.3005464 0.2092071 
Vemos que 30,0% dos cientistas da amostra não publicaram artigos nos últimos três anos do seu Ph.D., mas o modelo Poisson prevê que apenas 20,9% não teriam publicações. Claramente o modelo subestima a probabilidade de zero contagens.
 Uma forma de modelar este tipo de situação é assumir que os dados provêm de uma mistura de duas populações, uma onde as contagens são sempre zero, e outra onde a contagem tem uma distribuição de Poisson com média μ. Neste modelo as contagens zero podem vir de qualquer população, enquanto as contagens positivas vêm apenas da segunda.
No contexto das publicações dos bioquímicos Ph.D., podemos imaginar que alguns tinham em mente trabalhos onde as publicações não seriam importantes, enquanto outros visavam trabalhos acadêmicos onde um registro de publicações era esperado. Os membros do primeiro grupo publicariam zero artigos, enquanto os membros do segundo grupo publicariam 0,1,2,..., uma contagem que pode ser assumida como tendo uma distribuição de Poisson.
A distribuição do resultado pode então ser modelada em termos de dois parâmetros, π a probabilidade de 'sempre zero', e μ, o número médio de publicações para aqueles que não fazem parte do grupo 'sempre zero'. Uma forma natural de introduzir covariates é modelar o logit da probabilidade π de sempre zero e o log da média μ para aqueles que não estão na classe sempre zero.
Este tipo de modelo pode ser ajustado em R usando o zeroinfl()função no pacote pscl. A fórmula do modelo pode ser especificada como de costume se as mesmas variáveis devem ser incluídas em ambas as equações. Caso contrário, pode-se fornecer dois conjuntos de preditores separados por uma barra vertical, digite ?zeroinfl para aprender mais.
Aqui está um modelo com inflação zero com todas as covariáveis em ambas as equações:
> library(pscl)> mzip summary(mzip)Call:zeroinfl(formula = art ~ fem + mar + kid5 + phd + ment, data = ab)Pearson residuals: Min 1Q Median 3Q Max -2.3253 -0.8652 -0.2826 0.5404 7.2976 Count model coefficients (poisson with log link): Estimate Std. Error z value Pr(>|z|) (Intercept) 0.640838 0.121307 5.283 1.27e-07 ***femWomen -0.209145 0.063405 -3.299 0.000972 ***marMarried 0.103751 0.071111 1.459 0.144565 kid5 -0.143320 0.047429 -3.022 0.002513 ** phd -0.006166 0.031008 -0.199 0.842378 ment 0.018098 0.002294 7.888 3.07e-15 ***Zero-inflation model coefficients (binomial with logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) -0.577059 0.509386 -1.133 0.25728 femWomen 0.109746 0.280082 0.392 0.69518 marMarried -0.354014 0.317611 -1.115 0.26502 kid5 0.217097 0.196482 1.105 0.26919 phd 0.001274 0.145263 0.009 0.99300 ment -0.134114 0.045243 -2.964 0.00303 **---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Number of iterations in BFGS optimization: 21 Log-likelihood: -1605 on 12 Df
Olhando para a equação inflacionada vemos que o único preditor significativo de estar na classe 'sempre zero' é o número de artigos publicados pelo mentor, com cada artigo do mentor associado a 12,6% de probabilidade menor de nunca publicar.
Locando na equação do número médio ou artigos entre aqueles que não estão na classe sempre zero, encontramos desvantagens significativas para as mulheres e cientistas com crianças menores de cinco anos, e um grande efeito positivo do número de publicações do mentor, com cada artigo associado a um aumento de 1,8% no número esperado de publicações.
Para verificar que o modelo resolve o problema do excesso de zeros, prevemos π e μ, e calculamos a probabilidade combinada de nenhuma publicação. Existem opções na função predict()função chamada "zero" e "count" para obtê-las. Há também uma opção "prob" para calcular a densidade prevista, mas isso é exagerado, pois só queremos a probabilidade de zero.
> pr mu zip mean(zip) 0.2985679
Então o modelo resolve o problema do excesso de zeros, prevendo que 29,9% dos bioquímicos não publicarão artigos, muito mais próximos do valor observado de 30,0%.
Existe também um modelo binomial negativo inflado de zero, que usa binomial negativo para a contagem na classe 'nem sempre zero'. Este modelo pode ser ajustado usando o parâmetro zeroinfl() com o parâmetro dist="negbin". Links alternativos para a equação inflacionar incluem a probit,que pode ser especificada usando o parâmetro link="probit". 
Comparação do modelo com AIC
Como acontece, para estes dados o binômio negativo resolve o problema também! Aqui está a probablidade de zero artigos no binômio negativo. Procedemos a partir dos primeiros princípios, mas também se poderia usar a densidade do binômio negativo incorporado
> munb theta znb # also dnbinom(0, mu=munb, size=theta)> mean(znb) 0.3035957
O modelo prevê que 30,4% dos bioquímicos não publicariam artigos nos últimos três anos do seu Ph.D., muito próximo do valor observado de 30,0%.
Para escolher entre o binômio negativo e os modelos sem insuflação, precisamos recorrer a outros critérios. Uma maneira muito simples de comparar modelos com diferentes números de parâmetros é calcular o critério de informação da Akaike (AIC), que definimos como
AIC = -2logL + 2p
 onde p é o número de parâmetros no modelo. O primeiro termo é essencialmente o desvio e o segundo uma penalização para o número de parâmetros, e pode ser computado para a maioria dos modelos usando a função AIC. Eu também o obtenho "à mão" para que se possa verificar o cálculo. Para os nossos dados
> c(nb=AIC(mnb), zip=AIC(mzip)) nb zip 3135.917 3233.546 > mzip$rank aic sapply(list(mnb,mzip), aic) 3133.917 3233.546
 Para este conjunto de dados o modelo binomial negativo é um claro vencedor em termos de parcimónia e bondade de ajuste. 
Outros critérios diagnósticos que podemos observar são a distribuição marginal das contagens previstas e observadas e as funções de variância. 
Modelos de zero-truncado e Hurdle
Outros modelos que não cobrimos são o Poisson zero-truncado e o binômio negativo, projetados para dados que não incluem zeros. Um exemplo comum é o tempo de permanência em um hospital, que é de pelo menos um dia. Uma abordagem sensata é encaixar um modelo de Poisson ou binômio negativo que exclui zero e redimensiona as outras probabilidades para somar a um. Deve-se ter cuidado na interpretação desses modelos porque μ não é o resultado esperado, mas a média de uma distribuição subjacente que inclui os zeros. Estes modelos são implementados na função vglm() no pacote VGAM, usando as famílias pospoisson e posnegbinomial.
 Uma abordagem alternativa ao excesso (ou falta) de zeros é usar um processo de dois estágios, com um modelo logit para distinguir entre zero e contagens positivas e, em seguida, um modelo Poisson zero-truncado ou modelo binomial negativo para as contagens positivas. No nosso exemplo poderíamos usar um modelo logit para diferenciar os que publicam dos que não o fazem, e depois um modelo Poisson truncado ou modelo binomial negativo para o número de artigos daqueles que publicam pelo menos um. Estes modelos são muitas vezes chamados de hurdle models. Eles podem ser ajustados em R usando o logit separado e modelos binomiais truncados zero ou modelos binomiais negativos e simplesmente adicionando os log-likelihoods, ou usando a função hurdle() no pacote pscl.
Comparando os modelos de hurdle e zero-inflated, eu acho a distinção entre zero e um ou mais para ser mais clara com os modelos hurdle, mas a interpretação da média é mais clara com os modelos zero-inflated.
O seguinte artigo do JSS tem uma discussão útil de todos estes modelos: Modelos de Regressão para Dados de Contagem em R.