¿Glm.nb (en R) da valores p inflados cuando se prueba bajo el nulo?

Question

Estoy tratando de hacer una simple prueba de glm.nb en R.

Simulo resultados de una binomial negativa utilizando rnegbin. Tengo una variable caso/control 0/1 y me gustaría comprobar su significación, pero me encuentro con que hay demasiados valores p significativos cuando compruebo datos nulos:

pvalues = c()
for (blah in 1:20000) {
  outcomes = rnegbin(20,mu=30,theta=5)
  casecontrol =  c(rep(0,10),rep(1,10))
  model = glm.nb(outcomes~casecontrol,maxit=1000)
  pvalues = c(pvalues,summary(model)$coefficients[2,4])
}
hist(pvalues,40)

¿Qué estoy haciendo mal? Tengo datos binomiales negativos y los estoy probando con regresión binomial negativa, ¿quizás hay algo fundamental que no estoy entendiendo? Estaría enormemente agradecido si alguien me puede apuntar en la dirección correcta.

Answer 1

Los valores p habituales basados en métodos de máxima verosimilitud suelen implicar estadísticas t construidas dividiendo los parámetros estimados por los errores estándar (Wald) y comparando después los resultados con la distribución normal estándar. Aquí hay dos fuentes de error: La primera es que, como en el modelo de regresión habitual, hay variabilidad en el error estándar que hace que la distribución t sea más apropiada que la distribución z; y el segundo error es que la distribución de la estimación es sólo aproximadamente normal. Ambos problemas disminuyen con tamaños de muestra mayores. A continuación se muestra un código modificado que aborda tanto los problemas de la distribución t frente a la distribución z como los del tamaño de las muestras.

pvalues = c()
tvalues = c()
ndiv2 = 10
for (blah in 1:20000) {
outcomes = rnegbin(2*ndiv2,mu=30,theta=5)
casecontrol =  c(rep(0,ndiv2),rep(1,ndiv2))
model = glm.nb(outcomes~casecontrol,maxit=1000)
pvalues = c(pvalues,summary(model)$coefficients[2,4])
tvalues = c(tvalues,summary(model)$coefficients[2,3])
}
pvalues1 = 2*(1 - pt(abs(tvalues), 2*(ndiv2-1)))
hist(pvalues,40)
hist(pvalues1, 40)
mean(pvalues <=0.05)  # should be close to .05
mean(pvalues  <= 0.01)  # should be close to .01
mean(pvalues  <=0.005)  # should be close to .005
mean(pvalues <= 0.001)  # should be close to .001
mean(pvalues1 <=0.05)  # should be close to .05
mean(pvalues1  <= 0.01)  # should be close to .01
mean(pvalues1  <=0.005)  # should be close to .005
mean(pvalues1 <= 0.001)  # should be close to .001

Incluso con tu pequeña muestra de 20 (y para ser claros, ese es el problema con tus resultados), los resultados basados en t parecen mucho mejores. Y si aumenta el tamaño de la muestra de 20 a 200 (cambiando ndiv2 a 100), los resultados son aún mejores.

Answer 2

La respuesta corta es que la prueba generará el valor p uniforme sólo si estima theta correctamente o si proporciona theta.

El tamaño de la muestra es pequeño, lo que hace que la estimación del parámetro de dispersión (theta) sea menos estable e imprecisa.

Puedes escribir algo que te permita rastrear lo que da los valores p pequeños:

library(MASS)
res = lapply(1:20000,function(i){
  set.seed(i)
  outcomes = rnegbin(20,mu=30,theta=5)
  casecontrol =  c(rep(0,10),rep(1,10))
  model = glm.nb(outcomes~casecontrol,maxit=1000)
  data.frame(
  seed=i,
  p= summary(model)$coefficients[2,4],
  converged = model$converged,
  theta = model$theta
  )
})

res = do.call(rbind,res)

head(res)
  seed         p converged     theta
1    1 0.5062888      TRUE  7.361704
2    2 0.5727485      TRUE  4.116351
3    3 0.6651575      TRUE 10.457000
4    4 0.9183633      TRUE  7.348471
5    5 0.1878434      TRUE  8.519955
6    6 0.3917041      TRUE  3.897681

hist(res$p,br=40)

Es más de 0,05 y 1,5 veces más de lo esperado:

mean(res$p<0.05)
[1] 0.0815

Podemos fijarnos en los significativos:

head(res[order(res$p),],10)
       seed            p converged     theta
10794 10794 8.936969e-09      TRUE 15.504781
18191 18191 3.835794e-07      TRUE 12.724549
8409   8409 6.447190e-07      TRUE 33.455136
6371   6371 6.618804e-07      TRUE 93.952097
496     496 7.851968e-07      TRUE 13.578130
5600   5600 1.606424e-06      TRUE  9.295402
8531   8531 3.123901e-06      TRUE  8.908264
9109   9109 3.126698e-06      TRUE 24.742166
1470   1470 4.151136e-06      TRUE 18.737336
17462 17462 4.298971e-06      TRUE 16.478784

Puede ver que las estimaciones de theta son mucho mayores, lo que indica que el modelo está estimando una dispersión menor (dispersión = 1/theta). Podemos ver cómo es la simulación:

set.seed(10794)
outcomes = rnegbin(20,mu=30,theta=5)
grp = c(rep(0,10),rep(1,10))
boxplot(outcomes ~ grp ,horizontal = TRUE)
rug(outcomes[grp==0],col="blue")
rug(outcomes[grp==1],col="red")

Lo que parece bastante decente si el modelo no sabe qué esperar para theta.

Vamos a proporcionar theta como simulado y realizar la prueba:

sim_p = sapply(1:20000,function(i){
  set.seed(i)
  outcomes = rnegbin(20,mu=30,theta=5)
  casecontrol =  c(rep(0,10),rep(1,10))
  model = glm(outcomes~casecontrol,maxit=1000,family=neg.bin(5))
  summary(model)$coefficients[2,4]
})

hist(sim_p,br=40)

mean(sim_p<0.05)
[1] 0.05425

Si va a utilizar este glm para realizar pruebas, piense si dispone de muestras suficientes para estimar theta correctamente.