¿Cómo obtener el valor p (verificar significancia) de un efecto en un modelo mixto lme4?

56

Yo uso lme4 en R para ajustar el modelo mixto

lmer(value~status+(1|experiment)))

donde el valor es continuo, el estado y el experimento son factores, y obtengo

Linear mixed model fit by REML 
Formula: value ~ status + (1 | experiment) 
  AIC   BIC logLik deviance REMLdev
 29.1 46.98 -9.548    5.911    19.1
Random effects:
 Groups     Name        Variance Std.Dev.
 experiment (Intercept) 0.065526 0.25598 
 Residual               0.053029 0.23028 
Number of obs: 264, groups: experiment, 10

Fixed effects:
            Estimate Std. Error t value
(Intercept)  2.78004    0.08448   32.91
statusD      0.20493    0.03389    6.05
statusR      0.88690    0.03583   24.76

Correlation of Fixed Effects:
        (Intr) statsD
statusD -0.204       
statusR -0.193  0.476

¿Cómo puedo saber que el efecto del estado es significativo? R informa solo valores y no valores .ptp

ECII
fuente
1
Siguiendo las respuestas proporcionadas a esta pregunta, uno se pregunta qué es lo que realmente interesa a OP aquí: probar los coeficientes contra un valor nulo (la prueba vainilla se hace en una regresión lineal regular contra un valor nulo ), o pruebas para minimizar la varianza (la prueba que obtenemos de los muchos tipos de ANOVA). Esos dos apuntan a cosas diferentes. Aquí se encuentra una respuesta esclarecedora, aunque no se trata de modelos de efectos mixtos . H 0 : β = β nulo FtH0:β=βnullF
Firebug

Respuestas:

61

Hay mucha información sobre este tema en las preguntas frecuentes de GLMM . Sin embargo, en su caso particular, sugeriría usar

library(nlme)
m1 <- lme(value~status,random=~1|experiment,data=mydata)
anova(m1)

porque no necesita nada de lo que lmerofrece (mayor velocidad, manejo de efectos aleatorios cruzados, GLMM ...). lmedebe darle exactamente el coeficiente y la varianza mismas estimaciones, pero tendrán también df de cómputo y valores de p para usted (que hacen tiene sentido en un "clásico" de diseño como el que parecen tener). También es posible que desee considerar el término aleatorio ~status|experiment(que permite la variación de los efectos de estado entre bloques, o que incluye de manera equivalente una interacción estado por experimento). Los pósters anteriores también son correctos porque sus testadísticas son tan grandes que su valor p definitivamente será <0.05, pero me imagino que le gustaría valores p "reales".

Ben Bolker
fuente
3
No sé sobre esta respuesta. lmerpodría informar fácilmente los mismos tipos de valores p pero no por razones válidas. Supongo que es el comentario de que hay valores p "reales" aquí que me molestan. Podría argumentar que puede encontrar un posible límite y que se pasa cualquier límite razonable. Pero no puede argumentar que hay un valor p real.
John
11
Para un diseño clásico (equilibrado, anidado, etc.) creo que puedo argumentar que existe un valor p real, es decir, una probabilidad de obtener una estimación de beta de una magnitud observada o mayor si la hipótesis nula (beta = 0) eran falsos ... lme4 no proporciona este denominador df, creo, porque es más difícil de detectar en general a partir de una estructura de modelo lme4 cuando el modelo especificado es uno donde funcionaría alguna heurística para calcular un denominador clásico df ...
Ben Bolker
intente en su summary(m1)lugar (lo uso con el paquete nlme)
jena
36

Podrías usar el paquete lmerTest . Simplemente lo instala / carga y los modelos lmer se extienden. Entonces eg

library(lmerTest)
lmm <- lmer(value~status+(1|experiment)))
summary(lmm)
anova(lmm)

le daría resultados con valores p. Si los valores p son la indicación correcta, se cuestiona un poco, pero si desea tenerlos, esta es la forma de obtenerlos.

pbx101
fuente
28

Si puede manejar el abandono de los valores p ( y debería hacerlo ), puede calcular una razón de probabilidad que representaría el peso de la evidencia del efecto del estado a través de:

#compute a model where the effect of status is estimated
unrestricted_fit = lmer(
    formula = value ~ (1|experiment) + status
    , REML = F #because we want to compare models on likelihood
)
#next, compute a model where the effect of status is not estimated
restricted_fit = lmer(
    formula = value ~ (1|experiment)
    , REML = F #because we want to compare models on likelihood
)
#compute the AIC-corrected log-base-2 likelihood ratio (a.k.a. "bits" of evidence)
(AIC(restricted_fit)-AIC(unrestricted_fit))*log2(exp(1))
Mike Lawrence
fuente
16
Tenga en cuenta que las razones de probabilidad son asintóticas, es decir, no tenga en cuenta la incertidumbre en la estimación de la varianza residual ...
Ben Bolker
55
Estoy interesado en tu última línea. ¿Cuál es la interpretación del resultado? ¿Hay fuentes en las que pueda echarle un vistazo?
mguzmann
13

El problema es que el cálculo de los valores p para estos modelos no es trivial, consulte la discusión aquí, por lo que los autores del lme4paquete han elegido deliberadamente no incluir valores p en la salida. Puede encontrar un método para calcularlos, pero no necesariamente serán correctos.

Michelle
fuente
9

Considera lo que estás preguntando. Si solo desea saber si el valor p general para el efecto del estado pasa algún tipo de valor de corte arbitrario, como 0.05, entonces eso es fácil. Primero, desea averiguar el efecto general. Podrías obtener eso de anova.

m <- lmer(...) #just run your lmer command but save the model
anova(m)

Ahora tienes un valor F. Puede tomar eso y buscarlo en algunas tablas F. Simplemente elija el mínimo denom. grados de libertad. El límite allí será alrededor de 20. Su F puede ser mayor que eso, pero podría estar equivocado. Incluso si no es así, mire el número de grados de libertad de un cálculo ANOVA convencional aquí usando el número de experimentos que tiene. Pegar ese valor en usted tiene alrededor de 5 para un límite. Ahora lo pasa fácilmente en su estudio. El "verdadero" df para su modelo será algo mayor que eso porque está modelando cada punto de datos en lugar de los valores agregados que modelaría un ANOVA.

Si realmente desea un valor p exacto, no existe tal cosa a menos que esté dispuesto a hacer una declaración teórica al respecto. Si lees Pinheiro y Bates (2001, y tal vez algunos libros más sobre el tema ... mira otros enlaces en estas respuestas) y te sale un argumento para un df específico, entonces podrías usar eso. Pero de todos modos no estás buscando un valor p exacto. Menciono esto porque, por lo tanto, no debe informar un valor p exacto, solo que se pasa su límite.

Realmente debería considerar la respuesta de Mike Lawrence porque la idea de quedarse con un punto de paso para los valores p como la información final y más importante para extraer de sus datos generalmente es errónea (pero puede que no sea así en su caso, ya que nosotros no ' Realmente no tengo suficiente información para saber). Mike está usando una versión para mascotas del cálculo de LR que es interesante, pero puede ser difícil encontrar mucha documentación al respecto. Si observa la selección e interpretación de modelos utilizando AIC, puede que le guste.

John
fuente
9

Editar: este método ya no es compatible con las versiones más recientes de lme4. Use el paquete lmerTest como se sugiere en esta respuesta por pbx101 .

Hay una publicación en la lista R del autor de lme4 sobre por qué no se muestran los valores p. Sugiere usar muestras MCMC en su lugar, lo que haces usando pvals.fnc del paquete languageR:

library("lme4")
library("languageR")
model=lmer(...)
pvals.fnc(model)

Consulte http://www2.hawaii.edu/~kdrager/MixedEffectsModels.pdf para obtener un ejemplo y detalles.

Jeff
fuente
3
lme4 ya no es compatible con esto. Esta publicación podría actualizarse para evitar que las personas tengan que descubrir esto como lo acabo de hacer.
timothy.s.lau
5

¿Le interesa saber si el efecto combinado de statustiene un efecto significativo value? Si es así, puede usar la Anovafunción en el carpaquete (no debe confundirse con la anovafunción en la base R).

dat <- data.frame(
  experiment = sample(c("A","B","C","D"), 264, replace=TRUE), 
  status = sample(c("D","R","A"), 264, replace=TRUE), 
  value = runif(264)   
)
require(lme4)
(fm <- lmer(value~status+(1|experiment), data=dat))

require(car)
Anova(fm)

Echa un vistazo ?Anovadespués de cargar el carpaquete.

smillig
fuente
¿Alguna idea de cómo car::Anova()evita los problemas difíciles que rodean el cálculo de los valores p que Michelle vincula?
Mike Lawrence
No lo hago, pero supongo que evita los problemas difíciles al ignorarlos. Después de releer la publicación original, siento que podría haber entendido mal la pregunta. Si el OP quiere valores p exactos para los parámetros de efectos fijos, está en problemas. Pero si el OP solo quiere saber si son significativos, creo que los valores t son mayores que cualquier incertidumbre sobre cómo se calcularía el valor p exacto. (En otras palabras, son importantes.)
smillig
1
Creo que definitivamente fue una buena idea redirigir a un cálculo ANOVA para averiguar el efecto general de las estadísticas, pero no estoy seguro de que el ajuste de los valores p sea bueno. El anovacomando regular te dará F's.
John
Creo que esto es un poco más pegajoso que aparente. Ejecutar ANOVA es válido cuando desea minimizar la varianza, pero a partir de la redacción de la pregunta, creo que OP desea establecer el efecto marginal de las variables, es decir , los coeficientes de prueba contra un valor nulo.
Firebug
0

La función pvals.fncya no es compatible con lme4. Al usar el paquete lmerTest, es posible usar otro método para calcular el valor p, como las aproximaciones de Kenward-Roger

model=lmer(value~status+1|experiment)
anova(model, ddf="Kenward-Roger")
Stefano Cacciatore
fuente
0

Simplemente cargando el paquete afex imprimirá los valores p en la salida de la función lmer del paquete lme4 (no es necesario que esté usando el afex; simplemente cárguelo):

library(lme4)  #for mixed model
library(afex)  #for p-values

Esto agregará automáticamente una columna de valor p a la salida del lmer (yourmodel) para los efectos fijos.

Maryam Nasseri
fuente