AER дисперсия () противоречит отрицательной биномиальной дисперсии в R - PullRequest
Новая
       98

AER дисперсия () противоречит отрицательной биномиальной дисперсии в R

0 голосов
/ 17 апреля 2020

Я анализирую пуассоновскую регрессию подсчета данных. Пуассон требует, чтобы дисперсия и среднее значение были равны, поэтому я проверяю дисперсию, чтобы убедиться в этом. Для дисперсии я использую два метода:

  • disperiontest () для пакета AER.
  • проверить моделирование дисперсии как отрицательный бином с (glm.nb) 
> pm <- glm(myCounts ~ campaign, d, family = poisson)
> summary(pm)

Call:
glm(formula = myCounts ~ campaign, family = poisson, data = d)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-4.074  -1.599  -0.251   1.636   6.399  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  4.955870   0.032174  154.03   <2e-16 ***
campaign    -0.025879   0.001716  -15.08   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 428.04  on 35  degrees of freedom
Residual deviance: 195.81  on 34  degrees of freedom
AIC: 426.37

Number of Fisher Scoring iterations: 4

> dispersiontest(pm)

    Overdispersion test

data:  pm
z = 3.1933, p-value = 0.0007032
alternative hypothesis: true dispersion is greater than 1
sample estimates:
dispersion 
   5.53987 

> # Calculate dispersion with Negative Binomial
> nb_reg <- glm.nb(myCounts ~ campaign, data=d)
> summary.glm(nb_reg)

Call:
glm.nb(formula = myCounts ~ campaign, data = d, init.theta = 22.0750109, 
    link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9235  -0.7083  -0.1776   0.6707   2.4495  

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.914728   0.082327  59.697  < 2e-16 ***
campaign    -0.023471   0.003965  -5.919  1.1e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(22.075) family taken to be 1.069362)

    Null deviance: 76.887  on 35  degrees of freedom
Residual deviance: 35.534  on 34  degrees of freedom
AIC: 325.76

Number of Fisher Scoring iterations: 1

Как видите, NB обеспечивает дисперсию 1,069362. Тем не менее, Dispertest () приводит к 5,5 с явной чрезмерной дисперсии. Если я не ошибаюсь, AER не является параметрическим c тестом, поэтому мы можем только знать, есть ли избыточная / недостаточная дисперсия или нет. Тем не менее оба метода противоречат друг другу.

Кто-нибудь знает почему?

1 Ответ

1 голос
/ 17 апреля 2020

В glm.nb () дисперсия параметризована как ? + ? ^ 2 / ?, где ? - ваше среднее значение (подробнее см. в этом обсуждении ), а ? - тета, тогда как в пуассоне оно равно ϕ * ?, где ϕ - дисперсия 5.53987, которую вы видите.

В отрицательном биноме дисперсия 1.069362 не будет иметь смысла, вам нужно взглянуть на тета внутри отрицательного бинома (), который в вашем случае равен 22,075. У меня нет ваших данных, но я использую ваш перехват как приблизительную оценку среднего значения: 

mu = exp(4.914728)
theta = 22.0750109
variance = mu + (mu^2)/theta
variance
977.6339
variance / mu
[1] 7.173598

Что дает вам нечто похожее на дисперсию, которую вы. Вы должны заметить, что ваша дисперсия оценивается по полной модели, тогда как я просто угадал ее по вашему перехвату.

Итог - результаты не противоречат. Вы можете использовать отрицательный бином для моделирования ваших данных.

Ниже приведен пример, который иллюстрирует вышеуказанное соотношение. Мы моделируем чрезмерно рассредоточенные данные с использованием отрицательного бинома (это самое простое): 

y = c(rnbinom(100,mu=100,size=22),rnbinom(100,mu=200,size=22))
x = rep(0:1,each=100)

AER::dispersiontest(glm(y~x,family=poisson))

    Overdispersion test

data:  glm(y ~ x, family = poisson)
z = 8.0606, p-value = 3.795e-16
alternative hypothesis: true dispersion is greater than 1
sample estimates:
dispersion 
  8.200214

Грубо говоря, это получается путем деления дисперсии в каждой группе на среднее значение в каждой группе: 

mean(tapply(y,x,var)/tapply(y,x,mean))
[1] 8.283044

И Вы можете видеть, что дисперсия показывает 1, когда на самом деле ваши данные слишком раздроблены: 

    summary(MASS::glm.nb(y~x))

Call:
MASS::glm.nb(formula = y ~ x, init.theta = 21.5193965, link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-4.0413  -0.6999  -0.1275   0.5800   2.4774  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  4.66551    0.02364  197.36   <2e-16 ***
x            0.66682    0.03274   20.37   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(21.5194) family taken to be 1)
...