Cómo conseguir el efecto global para el modelo lineal mixto en lme4 en r?

Question

Me gustaría obtener el valor de significación y tamaño del efecto de la variable independiente en general, en lugar de la salida normal de lme4 en R. es igual a la cosa que la gente se informe cuando se ejecuta ANOVA. ¿Tiene usted alguna idea de como puedo conseguir esto?

Answer 1

Tanto de los conceptos que usted menciona (valores de p y de los tamaños del efecto lineal de los modelos mixtos) tienen problemas inherentes. Con respecto al tamaño del efecto, citando a Doug Bates, el autor original de la lme4,

Suponiendo que uno quiere definir un $R^2$ medida, creo que un argumento podría hacerse para el tratamiento de la penalizado suma residual de los cuadrados de un modelo lineal mixto de la misma manera que consideramos la suma residual de plazas a partir de un modelo lineal. O uno podría utilizar sólo la suma residual de plazas sin que la sanción o el mínimo de la suma de cuadrados residual obtener a partir de un conjunto dado de condiciones, que corresponde a un infinito la precisión de la matriz. No sé, la verdad. Depende de lo que se tratando de caracterizar.

Para más información, puedes mirar en este hilo, este hilo, y este mensaje. Básicamente, el problema es que no hay un acuerdo sobre el método para la inclusión y la descomposición de la varianza de los efectos aleatorios en el modelo. Sin embargo, hay algunas normas que se utilizan. Si usted tiene un vistazo a la Wiki configurar para/por el r-sig-mixto-modelos de lista de correo, hay un par de enfoques mencionados.

Uno de los métodos sugeridos se ve en la correlación entre los armarios y los valores observados. Esto puede ser implementado en R como se sugiere por Jarrett Byrnes en uno de esos hilos:

r2.corr.mer <- function(m) {
  lmfit <-  lm(model.response(model.frame(m)) ~ fitted(m))
  summary(lmfit)$r.squared
}

Así, por ejemplo, decir que estima el siguiente modelo lineal mixto:

set.seed(1)
d <- data.frame(y = rnorm(250), x = rnorm(250), z = rnorm(250),
                g = sample(letters[1:4], 250, replace=T)       )
library(lme4)
summary(fm1 <- lmer(y ~ x + (z | g), data=d))
# Linear mixed model fit by REML ['lmerMod']
# Formula: y ~ x + (z | g)
#    Data: d
# REML criterion at convergence: 744.4
# 
# Scaled residuals: 
#     Min      1Q  Median      3Q     Max 
# -2.7808 -0.6123 -0.0244  0.6330  3.5374 
# 
# Random effects:
#  Groups   Name        Variance Std.Dev. Corr 
#  g        (Intercept) 0.006218 0.07885       
#           z           0.001318 0.03631  -1.00
#  Residual             1.121439 1.05898       
# Number of obs: 250, groups: g, 4
# 
# Fixed effects:
#             Estimate Std. Error t value
# (Intercept)  0.02180    0.07795   0.280
# x            0.04446    0.06980   0.637
# 
# Correlation of Fixed Effects:
#   (Intr)
# x -0.005

Se puede calcular el tamaño del efecto mediante la función definida anteriormente:

r2.corr.mer(fm1)
# [1] 0.0160841

Una alternativa similar es recomendado en un papel por Ronghui Xu, que se conoce como $\Omega^{2}_{0}$, y puede ser calculado en R simplemente:

1-var(residuals(fm1))/(var(model.response(model.frame(fm1))))
# [1] 0.01173721  # Usually, it would be even closer to the value above

Con respecto a los valores de p, esto es mucho más polémico tema (al menos en el R/lme4 de la comunidad). Ver las discusiones en las preguntas aquí, aquí, y aquí , entre muchos otros. Hace referencia a la página de la Wiki de nuevo, hay un par de enfoques para poner a prueba hipótesis sobre los efectos en los lineales de los modelos mixtos. Listado de "peor a mejor" (según los autores de la página de la Wiki que creo que incluye Doug Bates así como Ben Bolker que contribuye mucho por aquí):

• Wald Z-pruebas
• Para equilibrado, anidada LMMs donde df puede ser calculada: Wald pruebas t
• Prueba de razón de verosimilitud, ya sea mediante la creación del modelo, por lo que el parámetro puede ser aislado/caída (via anova o drop1), o a través de la computación probabilidad de perfiles
• MCMC o paramétricas intervalos de confianza bootstrap

Se recomienda la cadena de Markov de Monte Carlo de muestreo de enfoque y también la lista de una serie de posibilidades para implementar este de pseudo y totalmente Bayesiano enfoques, se enumeran a continuación.

Pseudo-Bayesiano:

• Post-hoc de muestreo, normalmente (1) suponiendo que el plano de los priores y (2) a partir de la MLE, posiblemente utilizando aproximado de varianza-covarianza de la estimación para elegir a un candidato de distribución
• A través de mcmcsamp (si está disponible para su problema: es decir, LMMs con sencillo de efectos aleatorios - no GLMMs o complejo de efectos aleatorios)
  A través de pvals.fnc en la languageR paquete, un contenedor para mcmcsamp)
• En AD Model Builder, posiblemente a través de la glmmADMB paquete (el uso de la mcmc=TRUE opción) o el R2admb paquete (escriba su propia definición de modelo en AD Model Builder), o fuera de R
• A través de la sim función de la arm paquete (simula la parte posterior sólo para la versión beta (fijo-efecto de los coeficientes de

Totalmente Bayesiano enfoques:

• A través de la MCMCglmm paquete
• El uso de glmmBUGS (un WinBUGS contenedor/R de la interfaz)
• El uso de PUNTAS/WinBUGS/OpenBUGS etc., a través de la rjags/r2jags/R2WinBUGS/BRugs paquetes

Por el bien de la ilustración para mostrar lo que podría parecer, a continuación es una MCMCglmm estimado de uso de la MCMCglmm paquete que usted verá los rendimientos de los resultados parecidos a los del modelo anterior y tiene algún tipo de Bayesiana p-valores:

library(MCMCglmm)
summary(fm2 <- MCMCglmm(y ~ x, random=~us(z):g, data=d))
# Iterations = 3001:12991
# Thinning interval  = 10
#  Sample size  = 1000 
# 
#  DIC: 697.7438 
# 
#  G-structure:  ~us(z):g
# 
#       post.mean  l-95% CI u-95% CI eff.samp
# z:z.g 0.0004363 1.586e-17 0.001268    397.6
# 
#  R-structure:  ~units
# 
#       post.mean l-95% CI u-95% CI eff.samp
# units    0.9466   0.7926    1.123     1000
# 
#  Location effects: y ~ x 
# 
#             post.mean l-95% CI u-95% CI eff.samp pMCMC
# (Intercept)  -0.04936 -0.17176  0.07502     1000 0.424
# x            -0.07955 -0.19648  0.05811     1000 0.214

Espero que esto ayude un poco. Creo que el mejor consejo para alguien que empieza con modelos mixtos lineales y tratando de estimar en R es leer la Wiki de preguntas frecuentes de donde la mayor parte de esta información fue extraída. Es un recurso excelente para todo tipo de efectos variados temas, desde básico a avanzado y, a partir de la elaboración de modelos de trazado.

Answer 2

Yo uso el lmerTest paquete. Este hotel incluye una estimación de la p-valor en la anova() de salida para mi negocio MULTINIVEL análisis, pero no se da un tamaño del efecto por las razones dadas en otros posts aquí.