GLM binomial negativo frente a transformación logarítmica para datos de recuento: tasa de error tipo I aumentada

Algunos de ustedes podrían haber leído este bonito artículo:

O'Hara RB, Kotze DJ (2010) No log-transforma los datos de recuento. Methods in Ecology and Evolution 1: 118–122. Klick .

En mi campo de investigación (ecotoxicología) estamos lidiando con experimentos pobremente replicados y los GLM no son ampliamente utilizados. Entonces hice una simulación similar a la de O'Hara y Kotze (2010), pero imité los datos ecotoxicológicos.

Simulaciones de poder :

Simulé datos de un diseño factorial con un grupo de control ( ) y 5 grupos de tratamiento ( ). La abundancia en el tratamiento 1 fue idéntica al control ( ), la abundancia en los tratamientos 2-5 fue la mitad de la abundancia en el control ( ). Para las simulaciones, varié el tamaño de la muestra (3,6,9,12) y la abundancia en el grupo de control (2, 4, 8, ..., 1024). Las abundancias se extrajeron de distribuciones binomiales negativas con un parámetro de dispersión fijo ( ). Se generaron y analizaron 100 conjuntos de datos utilizando un GLM binomial negativo y un GLG + gaussiano transformado de datos.μ 1 - 5 μ 1 = μ c μ 2 - 5 = 0.5 μ c θ = $\mu_c$$\mu_{1-5}$$\mu_1 = \mu_c$$\mu_{2-5} = 0.5 \mu_c$$\theta = 3.91$

Los resultados son los esperados: el GLM tiene mayor poder, especialmente cuando no se tomaron muestras de muchos animales. El código está aquí

Error tipo I :

Luego miré el error tipo uno. Las simulaciones se realizaron como anteriormente, sin embargo, todos los grupos tuvieron la misma abundancia ( ).$\mu_c = \mu_{1-5}$

Sin embargo, los resultados no son los esperados: la GLM binomial negativa mostró un mayor error de Tipo I en comparación con la transformación LM +. Como se esperaba, la diferencia se desvaneció al aumentar el tamaño de la muestra. El código está aquí

Pregunta:

¿Por qué hay un error de tipo I aumentado en comparación con la transformación lm +?

Si tenemos datos deficientes (tamaño de muestra pequeño, poca abundancia (muchos ceros)), ¿deberíamos usar la transformación lm +? Los tamaños de muestra pequeños (2-4 por tratamiento) son típicos para tales experimentos y no se pueden aumentar fácilmente.

Aunque, el neg. compartimiento. GLM puede justificarse como apropiado para estos datos, la transformación lm + puede evitar que cometamos errores de tipo 1.

Este es un problema extremadamente interesante. Revisé su código y no puedo encontrar ningún error tipográfico inmediatamente obvio.

Me gustaría verte rehacer esta simulación, pero usa la prueba de máxima verosimilitud para hacer inferencia sobre la heterogeneidad entre grupos. Esto implicaría reajustar un modelo nulo para que pueda obtener estimaciones de los s bajo la hipótesis nula de homogeneidad en las tasas entre grupos. Creo que esto es necesario porque el modelo binomial negativo no es un modelo lineal (la tasa se parametriza linealmente, pero los s no lo son). Por lo tanto, no estoy convencido de que el argumento proporcione la inferencia correcta.$\theta$$\theta$drop1

La mayoría de las pruebas para modelos lineales no requieren que vuelva a calcular el modelo bajo la hipótesis nula. Esto se debe a que puede calcular la pendiente geométrica (prueba de puntaje) y aproximar el ancho (prueba de Wald) utilizando estimaciones de parámetros y covarianza estimada solo bajo la hipótesis alternativa.

Como el binomio negativo no es lineal, creo que necesitará ajustar un modelo nulo.

EDITAR:

Edité el código y obtuve lo siguiente:

Código editado aquí: https://github.com/aomidpanah/simulations/blob/master/negativeBinomialML.r

El artículo de O'Hara y Kotze (Métodos en Ecología y Evolución 1: 118–122) no es un buen punto de partida para la discusión. Mi mayor preocupación es la afirmación en el punto 4 del resumen:

Encontramos que las transformaciones funcionaron mal, excepto. . .. Los modelos cuasi-Poisson y binomial negativo ... [mostraron] poco sesgo.

La media para una distribución binomial de Poisson o negativa es para una distribución que, para los valores de <= 2 y para el rango de valores de la media que se investigó, es muy sesgada. Las medias de las distribuciones normales ajustadas están en una escala de log (y + c) (c es el desplazamiento), y estiman E (log (y + c)]. Esta distribución está mucho más cerca de ser simétrica que la distribución de y .θ $\lambda$$\theta$$\lambda$

Las simulaciones de O'Hara y Kotze comparan E (log (y + c)], según lo estimado por la media (log (y + c)), con log (E [y + c]). Pueden ser, y en los casos notados son muy diferentes. Sus gráficos no comparan un binomio negativo con un ajuste log (y + c), sino que comparan la media (log (y + c)] con log (E [y + c]). En el registro ( ) que se muestran en sus gráficos, ¡en realidad son los ajustes binomiales negativos los que están más sesgados! $\lambda$

El siguiente código R ilustra el punto:

x <- rnbinom(10000, 0.5, mu=2)  
## NB: Above, this 'mu' was our lambda. Confusing, is'nt it?
log(mean(x+1))
[1] 1.09631
log(2+1)  ## Check that this is about right
[1] 1.098612

mean(log(x+1))
[1] 0.7317908

O tratar

log(mean(x+.5))
[1] 0.9135269
mean(log(x+.5))
[1] 0.3270837

¡La escala en la que se estiman los parámetros es muy importante!

Si se toma una muestra de un Poisson, por supuesto, se espera que el Poisson funcione mejor, si se juzga por los criterios utilizados para ajustarse al Poisson. Lo mismo para un binomio negativo. La diferencia puede no ser tan grande, si la comparación es justa. Los datos reales (por ejemplo, tal vez, en algunos contextos genéticos) a veces pueden estar muy cerca de Poisson. Cuando parten de Poisson, el binomio negativo puede o no funcionar bien. Del mismo modo, especialmente si es del orden de quizás 10 o más, para modelar log (y + 1) utilizando la teoría normal estándar.$\lambda$

Tenga en cuenta que los diagnósticos estándar funcionan mejor en una escala de registro (x + c). La elección de c puede no importar demasiado; a menudo 0.5 o 1.0 tienen sentido. También es un mejor punto de partida para investigar las transformaciones de Box-Cox, o la variante Yeo-Johnson de Box-Cox. [Yeo, I. y Johnson, R. (2000)]. Vea además la página de ayuda para powerTransform () en el paquete de auto de R. El paquete gamlss de R permite ajustar los tipos binomiales negativos I (la variedad común) o II, u otras distribuciones que modelan la dispersión y la media, con enlaces de transformación de potencia de 0 (= log, es decir, enlace de registro) o más . Los ajustes pueden no siempre converger.

Ejemplo: Los datos de muertes vs daños base son para huracanes atlánticos con nombre que llegaron al territorio continental de EE. UU. Los datos están disponibles (nombre hurricNamed ) de una versión reciente del paquete DAAG para R. La página de ayuda para los datos tiene detalles.

El gráfico compara una línea ajustada obtenida usando un ajuste de modelo lineal robusto, con la curva obtenida transformando un ajuste binomial negativo con enlace de registro en la escala de registro (conteo + 1) utilizada para el eje y en el gráfico. (Tenga en cuenta que uno debe usar algo similar a una escala logarítmica (conteo + c), con c positiva, para mostrar los puntos y la "línea" ajustada del ajuste binomial negativo en el mismo gráfico). Tenga en cuenta el gran sesgo que es evidente por el ajuste binomial negativo en la escala logarítmica. El ajuste robusto del modelo lineal es mucho menos sesgado en esta escala, si se supone una distribución binomial negativa para los recuentos. Un ajuste de modelo lineal sería imparcial bajo los supuestos de la teoría normal clásica. ¡Encontré el sesgo sorprendente cuando creé por primera vez lo que era esencialmente el gráfico anterior! Una curva encajaría mejor con los datos, pero la diferencia está dentro de los límites de los estándares habituales de variabilidad estadística. El ajuste robusto del modelo lineal hace un mal trabajo para los recuentos en el extremo inferior de la escala.

Nota --- Estudios con datos de RNA-Seq: La comparación de los dos estilos de modelo ha sido de interés para el análisis de datos de conteo de experimentos de expresión génica. El siguiente artículo compara el uso de un modelo lineal robusto, trabajando con log (cuenta + 1), con el uso de ajustes binomiales negativos (como en el paquete de bioconductores edgeR ). La mayoría de los recuentos, en la aplicación RNA-Seq que se tiene en cuenta principalmente, son lo suficientemente grandes como para que el modelo log-lineal pesado adecuadamente funcione extremadamente bien.

Law, CW, Chen, Y, Shi, W, Smyth, GK (2014). Voom: los pesos de precisión desbloquean herramientas de análisis de modelos lineales para conteos de lectura de RNA-seq. Genome Biology 15, R29. http://genomebiology.com/2014/15/2/R29

NB también el artículo reciente:

Schurch NJ, Schofield P, Gierliński M, Cole C, Sherstnev A, Singh V, Wrobel N, Gharbi K, Simpson GG, Owen-Hughes T, Blaxter M, Barton GJ (2016). ¿Cuántas réplicas biológicas se necesitan en un experimento de RNA-seq y qué herramienta de expresión diferencial debe usar? ARN http://www.rnajournal.org/cgi/doi/10.1261/rna.053959.115

Es interesante que el modelo lineal se ajuste utilizando el paquete de limma (como edgeR , del grupo WEHI) se mantiene extremadamente bien (en el sentido de mostrar poca evidencia de sesgo), en relación con los resultados con muchas réplicas, ya que el número de réplicas es reducido.

Código R para el gráfico anterior:

library(latticeExtra, quietly=TRUE)
hurricNamed <- DAAG::hurricNamed
ytxt <- c(0, 1, 3, 10, 30, 100, 300, 1000)
xtxt <- c(1,10, 100, 1000, 10000, 100000, 1000000 )
funy <- function(y)log(y+1)
gph <- xyplot(funy(deaths) ~ log(BaseDam2014), groups= mf, data=hurricNamed,
   scales=list(y=list(at=funy(ytxt), labels=paste(ytxt)),
           x=list(at=log(xtxt), labels=paste(xtxt))),
   xlab = "Base Damage (millions of 2014 US$); log transformed scale",
   ylab="Deaths; log transformed; offset=1",
   auto.key=list(columns=2),
   par.settings=simpleTheme(col=c("red","blue"), pch=16))
gph2 <- gph + layer(panel.text(x[c(13,84)], y[c(13,84)],
           labels=hurricNamed[c(13,84), "Name"], pos=3,
           col="gray30", cex=0.8),
        panel.text(x[c(13,84)], y[c(13,84)],
           labels=hurricNamed[c(13,84), "Year"], pos=1, 
           col="gray30", cex=0.8))
ab <- coef(MASS::rlm(funy(deaths) ~ log(BaseDam2014), data=hurricNamed))

gph3 <- gph2+layer(panel.abline(ab[1], b=ab[2], col="gray30", alpha=0.4))
## 100 points that are evenly spread on a log(BaseDam2014) scale
x <- with(hurricNamed, pretty(log(BaseDam2014),100))
df <- data.frame(BaseDam2014=exp(x[x>0])) 
hurr.nb <- MASS::glm.nb(deaths~log(BaseDam2014), data=hurricNamed[-c(13,84),])
df[,'hatnb'] <- funy(predict(hurr.nb, newdata=df, type='response'))
gph3 + latticeExtra::layer(data=df,
       panel.lines(log(BaseDam2014), hatnb, lwd=2, lty=2, 
           alpha=0.5, col="gray30"))    

La publicación original refleja el artículo de Tony Ives: Ives (2015) . Está claro que las pruebas de significación dan resultados diferentes a la estimación de parámetros.

John Maindonald explica por qué las estimaciones son parciales, pero su ignorancia de los antecedentes es molesto: nos critica por mostrar que un método que todos estamos de acuerdo es erróneo. Muchos ecologistas registran ciegamente la transformación, y estábamos tratando de señalar los problemas para hacerlo.

Aquí hay una discusión más matizada: Warton (2016)

Ives, AR (2015), para probar la importancia de los coeficientes de regresión, continúe y log-transforme los datos de conteo. Métodos Ecol Evol, 6: 828–835. doi: 10.1111 / 2041-210X.12386

Warton, DI, Lyons, M., Stoklosa, J. e Ives, AR (2016), Tres puntos a considerar al elegir una prueba LM o GLM para los datos de conteo. Métodos Ecol Evol. doi: 10.1111 / 2041-210X.12552