Quando se trata de modelagem conta (ou seja, números inteiros maiores ou iguais a 0), muitas vezes começamos com a regressão de Poisson. Este é um modelo linear generalizado onde uma resposta é suposta ter uma distribuição de Poisson condicional a uma soma ponderada de preditores. Por exemplo, podemos modelar o número de concussões documentadas para quarterbacks da NFL em função de snaps jogados e da experiência total de anos da sua linha ofensiva. No entanto, uma desvantagem potencial da regressão de Poisson é que ela pode não descrever com precisão a variabilidade das contagens.

Uma distribuição de Poisson é parametrizada por {lambda}, que por acaso é tanto a sua média como a sua variância. Embora conveniente de lembrar, muitas vezes não é realista. Uma distribuição de contagens geralmente terá uma variância que não é igual à sua média. Quando vemos isso acontecer com dados que assumimos (ou esperamos) ser Poisson distribuídos, dizemos que temos sub ou superdispersão, dependendo se a variância é menor ou maior do que a média. A realização da regressão de Poisson em dados de contagem que exibem esse comportamento resulta em um modelo que não se encaixa bem.

Uma abordagem que aborda esta questão é a Regressão Binomial Negativa. A distribuição binomial negativa, como a distribuição de Poisson, descreve as probabilidades da ocorrência de números inteiros maiores ou iguais a 0. Ao contrário da distribuição de Poisson, a variância e a média não são equivalentes. Isto sugere que ela pode servir como uma aproximação útil para a modelagem conta com variabilidade diferente da sua média. A variância de uma distribuição binomial negativa é uma função da sua média e tem um parâmetro adicional, k, chamado parâmetro de dispersão. Digamos que nossa contagem é a variável aleatória Y de uma distribuição binomial negativa, então a variância de Y é

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

> Quando o parâmetro de dispersão fica cada vez maior, a variância converge para o mesmo valor da média, e o binômio negativo se transforma em uma distribuição de Poisson.

Para ilustrar a distribuição binomial negativa, vamos trabalhar com alguns dados do livro, Categorical Data Analysis, de Alan Agresti (2002). Os dados são apresentados na Tabela 13.6 na seção 13.4.3. Os dados são de um inquérito a 1308 pessoas em que lhes foi perguntado quantas vítimas de homicídio conhecem. As variáveis são resp, o número de vítimas que o respondente conhece, e race, a raça do respondente (preto ou branco). A raça ajuda a explicar quantas vítimas de homicídios uma pessoa conhece? Os dados primeiro precisam ser inseridos em 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)

Antes de chegarmos à modelagem, vamos explorar os dados. Primeiro notamos que a maioria dos entrevistados são brancos:

> table(race)racewhite black 1149 159 

Os pretos têm uma contagem média mais alta do que os brancos:

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

Para cada raça a variância da amostra é aproximadamente o dobro da média. Parece que temos sobredispersão.

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

Finalmente, olhamos para a distribuição das contagens por raça.

> 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

Em relação à montagem do modelo. Primeiro tentamos a regressão de Poisson usando a função glm() e mostramos uma parte da saída resumida.

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

Raça é muito significativa. Parece que os negros são muito mais propensos a conhecer alguém que tenha sido vítima de um homicídio. Mas o que significa o coeficiente 1,73? Neste modelo simples com um preditor dicotómico, é a diferença em log de contagens esperadas. Se exponenciarmos o coeficiente obtemos uma razão de médias amostrais:

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

De facto, se fizermos uma previsão com este modelo e exponenciarmos os resultados, obtemos as médias amostrais:

> 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 

Isto diz que a contagem de vítimas conhecidas para brancos é distribuída como um Poisson com média e variância igual a 0.09, enquanto a contagem de vítimas conhecidas para negros é distribuída como um Poisson com média e variância igual a 0,52. Já sabemos pela nossa análise exploratória que as variâncias observadas eram muito maiores, portanto não devemos estar muito satisfeitos com as variâncias estimadas do modelo. Se examinarmos as contagens ajustadas, veremos ainda mais evidências da falta de ajuste:

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

Acima salvamos primeiro as médias previstas em um objeto chamado fmeans. Em seguida, geramos contagens ajustadas usando a função dpois juntamente com os meios estimados para prever a probabilidade de se obter de 0 a 6. Em seguida, multiplicamos essas probabilidades pelo número de respondentes para obter contagens ajustadas. Finalmente, combinamos tudo em um quadro de dados para comparar facilmente os valores observados e ajustados. Duas das coisas mais dramáticas a notar é que estamos subajustando as contagens 0 e superajustando as contagens 1.

Podemos usar um rootograma para visualizar o ajuste de um modelo de regressão de contagem. A função rootogram() no pacote countreg torna isto fácil.

> countreg::rootogram(pGLM)

A linha curva vermelha é o ajuste teórico de Poisson. “Pendurado” a partir de cada ponto da linha vermelha é uma barra, cuja altura representa a diferença entre as contagens esperadas e observadas. Uma barra pendurada abaixo de 0 indica um encaixe inferior. Uma barra pendurada acima de 0 indica um ajuste excessivo. As contagens foram transformadas com uma transformação da raiz quadrada para evitar que contagens menores sejam obscurecidas e sobrecarregadas por contagens maiores. Vemos uma grande quantidade de subajustamento para as contagens 2 e superiores e um enorme sobreajustamento para a contagem 1.

Agora vamos tentar encaixar um modelo binomial negativo. Notamos que a variabilidade das contagens era maior para as duas raças. Parece que a distribuição binomial negativa se aproximaria melhor da distribuição das contagens.

Para encaixar um modelo binomial negativo em R voltamos para a função glm.nb() no pacote MASS (um pacote que vem instalado com R). Novamente mostramos apenas parte do resumo da saída:

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

Primeiro notar que os coeficientes são os mesmos de antes. Mais uma vez podemos exponenciar o coeficiente de corrida para obter uma razão de meios de amostra e fazer previsões para obter os meios de amostra originais.

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

Mas observe que o erro padrão para o coeficiente da raça é maior, indicando mais incerteza em nossa estimativa (0,24 versus 0,15). Isto faz sentido, dada a variabilidade observada em nossas contagens. Repare também na estimativa de Theta. Esse é o nosso parâmetro de dispersão. Podemos acessá-lo diretamente do nosso objeto modelo da seguinte forma:

> nbGLM$theta 0.2023119

E podemos usá-lo para obter as variâncias estimadas para as contagens:

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

Estas estão muito mais próximas das variâncias observadas do que as dadas pelo modelo de Poisson.

Após novamente visualizamos o ajuste usando um rootograma:

> countreg::rootogram(nbGLM)

Esta parece muito melhor do que o rootograma do modelo de Poisson. Há um ligeiro sub-conjunto/sobre-conjunto para as contagens de 1 a 3, mas de resto parece muito bom.

Para obter mais informações sobre o nosso modelo binomial negativo, vamos usar seus parâmetros para simular dados e comparar os dados simulados com os dados observados. Abaixo nós carregamos o pacote magrittr para acesso ao operador %>% que nos permite funções de “corrente”.

Primeiro tabulamos as contagens e criamos um barplot para os participantes brancos e pretos, respectivamente. Depois usamos os parâmetros do modelo para simular os dados de uma distribuição binomial negativa. Os dois parâmetros são mu e tamanho (ou seja, parâmetro de dispersão). Note que usamos a função coef() para extrair os coeficientes apropriados para cada prova. Para branco é apenas a intercepção, para preto é a intercepção e a inclinação (assim somamo-los). Depois exponenciamos para converter da escala de log para a escala original. Simulamos o mesmo número de observações que temos em nossos dados originais.

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)

Os dados simulados são muito parecidos com os dados observados, dando-nos novamente confiança na escolha da regressão binomial negativa para modelar estes dados. Estas parcelas também demonstram a natureza condicional do nosso modelo. A distribuição binomial negativa das contagens depende, ou é condicionada, da raça. Cada raça tem uma média diferente, mas um parâmetro de dispersão comum.

• 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

Para perguntas ou esclarecimentos sobre este artigo, entre em contato com a Biblioteca UVA StatLab: [email protected]

Veja toda a coleção de artigos da Biblioteca UVA StatLab.

Clay Ford
Statistical Research Consultant
University of Virginia Library