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

Question

Algunos de ustedes habrán leído este bonito artículo:

O'Hara RB, Kotze DJ (2010) Do not log-transform count data. Métodos en Ecología y Evolución 1:118-122. klick .

En mi campo de investigación (ecotoxicología) nos enfrentamos a experimentos poco replicados y los MLG no se utilizan mucho. Así que hice una simulación similar a la de O'Hara y Kotze (2010), pero imitando los datos ecotoxicológicos.

Simulaciones de potencia :

He simulado los datos de un diseño factorial con un grupo de control ( $\mu_c$ ) y 5 grupos de tratamiento ( $\mu_{1-5}$ ). La abundancia en el tratamiento 1 fue idéntica a la del control ( $\mu_1 = \mu_c$ ), la abundancia en los tratamientos 2-5 fue la mitad de la abundancia en el control ( $\mu_{2-5} = 0.5 \mu_c$ ). 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 una distribución binomial negativa con un parámetro de dispersión fijo ( $\theta = 3.91$ ). Se generaron 100 conjuntos de datos y se analizaron utilizando un MLG binomial negativo y un MLG gaussiano + datos transformados logarítmicamente.

Los resultados son los esperados: el MLG tiene mayor potencia, especialmente cuando no se muestrearon muchos animales. El código está aquí.

Error de tipo I :

A continuación me fijé en el error de tipo uno. Las simulaciones se hicieron como arriba, pero todos los grupos tenían la misma abundancia ( $\mu_c = \mu_{1-5}$ ).

Sin embargo, los resultados no son los esperados: El MLG binomial negativo mostró un mayor error de tipo I en comparación con el ML + transformación. Como era de esperar, la diferencia desaparecía al aumentar el tamaño de la muestra. El código está aquí.

Pregunta:

¿Por qué aumenta el error de tipo I en comparación con lm+transformación?

Si tenemos datos deficientes (tamaño de muestra pequeño, abundancia baja (muchos ceros)), ¿debemos utilizar entonces lm+transformación? Los tamaños de muestra pequeños (2-4 por tratamiento) son típicos de este tipo de experimentos y no pueden aumentarse fácilmente.

Aunque, el neg. bin. GLM puede justificarse como apropiado para estos datos, la transformación lm + puede evitarnos errores de tipo 1.

Answer 1

Se trata de un problema sumamente interesante. He revisado tu código y no encuentro ninguna errata evidente.

Me gustaría que volvieras a hacer esta simulación pero utilizando el prueba de máxima verosimilitud para hacer inferencias sobre la heterogeneidad entre grupos. Esto implicaría volver a ajustar un modelo nulo para poder obtener estimaciones de la $\theta$ s bajo la hipótesis nula de homogeneidad de tasas entre grupos. Creo que esto es necesario porque el modelo binomial negativo es no un modelo lineal (la tasa se parametriza linealmente, pero la $\theta$ s no lo son). Por lo tanto, no estoy convencido drop1 proporciona una inferencia correcta.

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 puntuación) y aproximar la amplitud (prueba de Wald) utilizando únicamente las estimaciones de los parámetros y la covarianza estimada bajo la hipótesis alternativa.

Dado que la binomial negativa no es lineal, creo que tendrá que ajustar un modelo nulo.

EDITAR:

He editado el código y he obtenido lo siguiente:

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

Answer 2

El artículo de O'Hara y Kotze (Methods in Ecology and Evolution 1:118-122) no es un buen punto de partida para el debate. Lo que más me preocupa es la afirmación del punto 4 del resumen:

Descubrimos que las transformaciones funcionaban mal, excepto... La modelos cuasi-Poisson y binomial negativo ... [mostraron] poco sesgo.

La media $\lambda$ para una distribución de Poisson o binomial negativa es para una distribución que, para valores de $\theta$ <= 2 y para el intervalo de valores de la media $\lambda$ que se investigó, tiene un sesgo altamente positivo. 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 se aproxima mucho más a la simetría que la distribución de y.

Las simulaciones de O'Hara y Kotze comparan E(log(y+c)], estimado por mean(log(y+c)), con log(E[y+c]). Pueden ser, y en los casos señalados lo son, muy diferentes. Sus gráficos no comparan una binomial negativa con un ajuste log(y+c), sino que comparan mean(log(y+c)] con log(E[y+c]). En el gráfico log( $\lambda$ ) mostrada en sus gráficos, ¡en realidad son los ajustes binomiales negativos los que están más sesgados!

El siguiente código R ilustra este 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 prueba

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

La escala a la que se estiman los parámetros es muy importante.

Si se toman muestras de un Poisson, por supuesto se espera que el Poisson lo haga mejor, si se juzga por los criterios utilizados para ajustar el Poisson. Lo mismo ocurre con una binomial negativa. 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 se alejan de Poisson, la binomial negativa puede o no funcionar bien. Del mismo modo, especialmente si $\lambda$ es del orden de quizá 10 o más, para modelizar log(y+1) utilizando la teoría normal estándar.

Tenga en cuenta que los diagnósticos estándar funcionan mejor en una escala de log(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 Box-Cox, o la variante Yeo-Johnson de Box-Cox. [Yeo, I. y Johnson, R. (2000)]. Véase también la página de ayuda de powerTransform() en el paquete car de R. El paquete gamlss de R permite ajustar binomias negativas de tipo I (la variedad común) o II, u otras distribuciones que modelan la dispersión además de la media, con enlaces de transformación de potencia de 0 (=log, es decir, enlace logarítmico) o más. Es posible que los ajustes no siempre converjan.

Ejemplo: Muertes vs Daño Base Los datos corresponden a huracanes con nombre en el Atlántico que alcanzaron el territorio continental de Estados Unidos. Los datos están disponibles (nombre huracánNombrado ) de una versión reciente del paquete DAAG para R. La página de ayuda de los datos contiene más detalles.

El gráfico compara una línea ajustada obtenida utilizando un ajuste de modelo lineal robusto, con la curva obtenida transformando un ajuste binomial negativo con enlace log en la escala log(recuento+1) utilizada para el eje y en el gráfico. (Obsérvese que hay que utilizar algo parecido a una escala log(recuento+c), con c positivo, para mostrar los puntos y la "línea" ajustada a partir del ajuste binomial negativo en el mismo gráfico). Obsérvese el gran sesgo evidente para el ajuste binomial negativo en la escala logarítmica. El ajuste del modelo lineal robusto es mucho menos sesgado en esta escala, si se asume una distribución binomial negativa para los recuentos. El ajuste del modelo lineal sería insesgado según los supuestos de la teoría normal clásica. Me sorprendió el sesgo cuando creé por primera vez lo que era esencialmente el gráfico anterior. Una curva se ajustaría mejor a los datos, pero la diferencia está dentro de los límites de las normas habituales de variabilidad estadística. El modelo lineal robusto no se ajusta bien a 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 recuento procedentes de experimentos de expresión génica. El siguiente artículo compara el uso de un modelo lineal robusto, trabajando con log(count+1), con el uso de ajustes binomiales negativos (como en el paquete Bioconductor edgeR ). La mayoría de los recuentos, en la aplicación RNA-Seq que se tiene en mente principalmente, son lo suficientemente grandes como para que los ajustes de modelos log-lineales adecuadamente ponderados funcionen extremadamente bien.

Law, CW, Chen, Y, Shi, W, Smyth, GK (2014). Voom: precision weights unlock linear model analysis tools for RNA-seq read counts. Biología del genoma 15, R29. http://genomebiology.com/2014/15/2/R29

NB también el documento 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 ARN-seq y qué herramienta de expresión diferencial debe utilizar? [ ] http://www.rnajournal.org/cgi/doi/10.1261/rna.053959.115

Es interesante que el modelo lineal se ajuste utilizando el limma (como edgeR del grupo WEHI) resisten muy bien (en el sentido de mostrar pocas pruebas de sesgo), en relación con los resultados con muchas réplicas, a medida que se reduce el número de réplicas.

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"))

Answer 3

El post original refleja la ponencia 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 sesgadas, pero molesta su ignorancia del trasfondo: nos critica por mostrar que un método que todos estamos de acuerdo en que es defectuoso, es defectuoso. Muchos ecologistas hacen transformaciones logarítmicas a ciegas, y nosotros intentábamos señalar los problemas de hacer eso.

Aquí hay un debate más matizado: Warton (2016)

Ives, A. R. (2015), For testing the significance of regression coefficients, go ahead and log-transform count data. Methods Ecol Evol, 6: 828-835. doi:10.1111/2041-210X.12386

Warton, D. I., Lyons, M., Stoklosa, J. and Ives, A. R. (2016), Three points to consider when choosing a LM or GLM test for count data. Methods Ecol Evol. doi:10.1111/2041-210X.12552