Es esta una forma aceptable para analizar los modelos de efectos mixtos con lme4 en R?

Question

Tengo un desbalance de medidas repetidas conjunto de datos a analizar, y he leído que la forma en que la mayoría de los paquetes estadísticos de manejar esto con ANOVA (es decir, el tipo III suma de cuadrados) está equivocado. Por lo tanto, me gustaría utilizar un modelo de efectos mixtos para el análisis de estos datos. He leído mucho acerca de los modelos mixtos en R, pero todavía soy muy nuevo en R y modelos de efectos mixtos y no muy seguro de que estoy haciendo las cosas bien. Tenga en cuenta que yo no puedo todavía completamente divorcio a mí mismo de métodos "tradicionales", y todavía necesita $p$-valores y post hoc de ensayos.

Me gustaría saber si el siguiente enfoque tiene sentido, o si estoy haciendo algo terriblemente mal. Aquí está mi código:

# load packages
library(lme4)
library(languageR)
library(LMERConvenienceFunctions)
library(coda)
library(pbkrtest)

# import data
my.data <- read.csv("data.csv")

# create separate data frames for each DV & remove NAs
region.data <- na.omit(data.frame(time=my.data$time, subject=my.data$subject, dv=my.data$dv1))

# output summary of data
data.summary <- summary(region.data)

# fit model
# "time" is a factor with three levels ("t1", "t2", "t3")
region.lmer <- lmer(dv ~ time + (1|subject), data=region.data)

# check model assumptions
mcp.fnc(region.lmer)

# remove outliers (over 2.5 standard deviations)
rm.outliers <- romr.fnc(region.lmer, region.data, trim=2.5)
region.data <- rm.outliers$data
region.lmer <- update(region.lmer)

# re-check model assumptions
mcp.fnc(region.lmer)

# compare model to null model
region.lmer.null <- lmer(dv ~ 1 + (1|subject), data=region.data)
region.krtest <- KRmodcomp(region.lmer, region.lmer.null)

# output lmer summary
region.lmer.summary <- summary(region.lmer)

# run post hoc tests
t1.pvals <- pvals.fnc(region.lmer, ndigits=10, withMCMC=TRUE)

region.lmer <- lmer(dv ~ relevel(time,ref="t2") + (1|subject), data=region.data)
t2.pvals <- pvals.fnc(region.lmer, ndigits=10, withMCMC=TRUE)

region.lmer <- lmer(dv ~ relevel(time,ref="t3") + (1|subject), data=region.data)
t3.pvals <- pvals.fnc(region.lmer, ndigits=10, withMCMC=TRUE)

# Get mcmc mean and 50/95% HPD confidence intervals for graphs
# repeated three times and stored in a matrix (not shown here for brevity)
as.numeric(t1.pvals$fixed$MCMCmean)
as.numeric(t1.pvals$fixed$HPD95lower)
as.numeric(t1.pvals$fixed$HPD95upper)
HPDinterval(as.mcmc(t1.pvals$mcmc),prob=0.5)
    HPDinterval(as.mcmc(t1.pvals$mcmc),prob=0.5)

Específicamente algunas de las preguntas que tengo:

1. Es esta una forma válida de análisis de modelos de efectos mixtos? Si no, ¿qué debo estar haciendo en su lugar.
2. Son las críticas parcelas de salida por mcp.fnc lo suficientemente bueno para la verificación de los supuestos del modelo, o debo de tomar medidas adicionales.
3. Entiendo que para los modelos mixtos para ser válidos, los datos deben respetar los supuestos de normalidad y homoscedasticity. Cómo puedo juzgar lo que es "aproximadamente normal" y qué no lo es por mirar el la crítica gráficos generados por mcp.fnc? Qué necesito para tener una idea para esto, o es su una forma establecida de hacer las cosas? Cómo robusto son modelos mixtos en relación a estos supuestos?
4. Necesito para evaluar las diferencias entre los tres puntos de tiempo de ~20 características (biomarcadores) de los sujetos en mi muestra. Es montaje y pruebas de modelos separados para cada aceptable tanto tiempo como yo informe llevado a cabo pruebas (significativas o no), o necesito alguna formulario de corrección para comparaciones múltiples.

Para ser un poco más precisos en lo que respecta a la experiencia, aquí hay algunos detalles más. Seguimos con un número de participantes longitudinalmente como se sometieron a un tratamiento. Hemos medido una serie de biomarcadores antes del inicio del tratamiento y en dos puntos de tiempo después. Lo que me gustaría ver es si hay diferencia en estos biomarcadores entre los tres puntos en el tiempo.

Me estoy basando la mayor parte de lo que yo estoy haciendo aquí, en este tutorial, pero hizo algunos cambios en base a mis necesidades y cosas que he leído. Los cambios que he hecho son:

1. relevel el factor "tiempo" para obtener t1-t2, t2-t3 y t1-t3 comparaciones con pvals.fnc (desde el languageR paquete)
2. comparar mi modelo mixto para el modelo nulo uso de un aproximado de F-test, basado en un Kenward-Roger el enfoque de la (utilizando el pbkrtest paquete) en lugar de una prueba de razón de verosimilitud (porque he leído, que Kenward-Roger es mejor considerado ahora)
3. El uso de la LMERConvenienceFunctions paquete para comprobar hipótesis y eliminar los valores atípicos (porque he leído que los modelos mixtos son muy sensibles a los valores extremos)

Answer 1

Su pregunta(s) es un poco "grande", así que voy a empezar con algunos comentarios generales y consejos.

Algunos antecedentes de la lectura y de paquetes útiles

Probablemente, usted debe echar un vistazo a algunos de los tutoriales de introducción a uso de modelos mixtos, yo recomiendo empezar con Baayen et al (2008) -- Baayen es el autor de languageR. Barr et al (2013) discutir algunos problemas con los efectos aleatorios de la estructura, y Ben Bolker es uno de los desarrolladores actuales de lme4. Baayen et al es especialmente bueno para sus preguntas, porque pasan mucho tiempo discutiendo el uso de bootstrap / pruebas de permutación (la materia detrás de mcp.fnc).

La literatura

• Baayen, R. H., D. Davidson und D. Bates (2008). De efectos mixtos de modelado con cruzado de efectos aleatorios para los sujetos y objetos. Diario de la Memoria y el Lenguaje, 59(390-412).
• Barr, Dale J, R. Levy, C. Scheepers und H. J. Tily (2013). Los efectos aleatorios de la estructura de confirmación de resultados de pruebas de hipótesis: Mantener la máxima. Diario de la Memoria y el Lenguaje, 68:255– 278.
• Bolker, Benjamin M, M. E. Brooks, C. J. Clark, S. W. Geange, J. R. Poulsen, M. H. H. Stevens und J.-S. S. White (2009). Los modelos mixtos lineales generalizados: una guía práctica para la ecología y la evolución. Tendencias Ecol Evol, 24(3):127-35.

Florian Jaeger también tiene un montón de cosas en el lado práctico de modelos mixtos en su laboratorio del blog.

También hay una serie de útiles R paquetes para la visualización y prueba de los modelos mixtos, como lmerTest y effects. El effects paquete es especialmente agradable, ya que permite trazar la regresión lineal y los intervalos de confianza de ir detrás de las escenas.

Bondad de ajuste y comparación de modelos

$p$-valores son realmente una pésima manera de comparar los modelos mixtos (y generalmente no son un gran método para nada). Hay una razón por qué Doug Bates (el autor original de lme4) siente que no es necesario incluir. Si realmente quieres ir de esa manera (y os ruego que no), la mencionada lmerTest proporciona algunas facilidades adicionales para el cálculo y tratar con ellos.

Mejores formas de comparación de modelos log-verosimilitud o la teoría de la información criterios como la AIC y BIC. De la AIC y BIC, la regla general es "más pequeño es mejor", pero tienes que tener cuidado ya que ambos penalizar más el modelo de grados de libertad y que no hay "absoluta" bueno o malo. Para obtener un buen resumen de la AIC y la log-verosimilitud de los modelos, puede utilizar el anova() función, que ha sido sobrecargado para aceptar mer objetos. Por favor, tenga en cuenta que el $\chi^2$ valores dados son válidos sólo para las comparaciones entre modelos anidados. Sin embargo, el resultado es bastante bueno para ser de manera tabular, incluso para otras comparaciones. Para modelos anidados, usted tiene un buen $\chi^2$ prueba de que no necesita de $p$-valores de comparar directamente los dos modelos. La desventaja de esto es que no es inmediatamente claro de lo bueno que es su ajuste es.

Para observar los efectos individuales (es decir, las cosas que se ven en una tradicional ANOVA), usted debe buscar en la $t$-valores de los efectos fijos en los modelos (que son parte de la summary() si no me equivoco). En general, nada de lo $|t|>2$ se considera buena (más detalles en Baayen et al). También puede acceder a los efectos fijos directamente con la función auxiliar fixef().

También debe asegurarse de que ninguno de sus efectos fijos son demasiado fuertemente correlacionada -- multicolinealidad viola los supuestos del modelo. Florian Jaeger ha escrito un poco sobre esto, así como un par de posibles soluciones.

Comparaciones múltiples

Voy a tratar de su pregunta #4 un poco más directamente. La respuesta es básicamente el mismo como una buena práctica con el tradicional análisis de la VARIANZA, por desgracia, este parece ser un lugar donde hay una gran cantidad de incertidumbre para la mayoría de los investigadores. (Es el mismo que el tradicional análisis de la VARIANZA debido a que ambos modelos de efectos mixtos lineal y análisis de la VARIANZA se basa en el modelo lineal general, es sólo que los modelos mixtos tienen un plazo adicional para los efectos al azar.) Si usted está asumiendo que las ventanas de tiempo de hacer una diferencia y quiero comparar los efectos del tiempo, se deberán incluir en el modelo. Esta, por cierto, también proporciona una forma conveniente para usted para juzgar si vez que hizo una diferencia: hay una principal (fijo) efecto por el tiempo? La desventaja de ir por esta ruta es que su modelo MUCHO más complejo y el single "súper modelo" con el tiempo como un paramater probablemente tomará más tiempo para calcular de tres modelos más pequeños sin tiempo como un paramater. De hecho, el clásico ejemplo de este tutorial para los modelos mixtos es sleepstudy que utiliza el tiempo como un paramater.

Yo no podía saber si las diferentes características están destinados a ser predictores o variables dependientes. Si están destinados a ser predictores, se puede tirar a todos dentro de un modelo y ver cuál tiene el mayor $t$-valor, pero este modelo sería increíblemente complejo, incluso si usted no permiten interacciones, y tomar un largo tiempo para calcular. El más rápido, y, potencialmente, la manera más fácil, sería para calcular los diferentes modelos para cada predictor. Me gustaría recomendar el uso de foreach bucle para paralelizar a través de tantos núcleos como usted tiene - lme4 aún no hacen su matriz de operaciones en paralelo, de manera que se puede calcular un modelo diferente en cada núcleo en paralelo. (Ver breve introducción aquí) a Continuación, puede comparar la AIC y BIC de cada modelo para encontrar cual es el mejor predictor. (No están anidados en este caso, por lo que'$\chi^2$ estadística.)

Si sus características son la variable dependiente, entonces usted tendrá que calcular los diferentes modelos de todos modos, y entonces usted puede utilizar el AIC y BIC para comparar los resultados.