¿Cómo obtener un valor p "general" y un tamaño de efecto para un factor categórico en un modelo mixto (lme4)?

Me gustaría obtener un valor p y un tamaño de efecto de una variable categórica independiente (con varios niveles), que es "general" y no para cada nivel por separado, como es la salida normal de lme4R. Es igual que lo que la gente informa cuando ejecuta un ANOVA.

¿Cómo puedo conseguir esto?

Los dos conceptos que menciona (valores p y tamaños de efectos de modelos lineales mixtos) tienen problemas inherentes. Con respecto al tamaño del efecto , citando a Doug Bates, el autor original de lme4,

Suponiendo que uno quiera definir una medida de , creo que podría hacerse un argumento para tratar la suma residual de cuadrados penalizada de un modelo mixto lineal de la misma manera que consideramos la suma residual de cuadrados de un modelo lineal. O se podría usar solo la suma residual de cuadrados sin la penalización o la suma mínima residual de cuadrados obtenible de un conjunto dado de términos, que corresponde a una matriz de precisión infinita. No lo se, de verdad. Depende de lo que estés tratando de caracterizar.$R^2$

Para obtener más información, puede ver este hilo , este hilo y este mensaje . Básicamente, el problema es que no existe un método acordado para la inclusión y descomposición de la varianza de los efectos aleatorios en el modelo. Sin embargo, hay algunos estándares que se utilizan. Si echa un vistazo a la Wiki creada para / por la lista de correo r-sig-mixed-models , hay un par de enfoques enumerados.

Uno de los métodos sugeridos analiza la correlación entre los valores ajustados y los observados. Esto se puede implementar en R como lo sugiere 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
}

Entonces, por ejemplo, supongamos que estimamos el siguiente modelo mixto lineal:

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

Podemos calcular el tamaño del efecto utilizando la función definida anteriormente:

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

$\Omega^{2}_{0}$

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 p , este es un tema mucho más polémico (al menos en la comunidad R / lme4). Vea las discusiones en las preguntas aquí , aquí y aquí, entre muchos otros. Haciendo referencia nuevamente a la página Wiki, hay algunos enfoques para probar hipótesis sobre los efectos en modelos lineales mixtos. Enumerado de "peor a mejor" (según los autores de la página Wiki, que creo que incluye a Doug Bates y a Ben Bolker, que contribuye mucho aquí):

• Pruebas Z de Wald
• Para LMM anidados y equilibrados donde se puede calcular df: pruebas t de Wald
• Prueba de razón de verosimilitud, ya sea configurando el modelo para que el parámetro pueda aislarse / descartarse (mediante anovao drop1), o mediante el cálculo de perfiles de probabilidad
• MCMC o intervalos de confianza de bootstrap paramétricos

Recomiendan el enfoque de muestreo Monte Carlo de la cadena de Markov y también enumeran una serie de posibilidades para implementar esto desde enfoques pseudo y completamente bayesianos, que se enumeran a continuación.

Pseudo-Bayesian:

• Muestreo post-hoc, típicamente (1) suponiendo anteriores planas y (2) comenzando desde el MLE, posiblemente usando la estimación aproximada de varianza-covarianza para elegir una distribución candidata
• Vía mcmcsamp(si está disponible para su problema: es decir, LMM con efectos aleatorios simples, no GLMM o efectos aleatorios complejos)
  Vía pvals.fncen el languageRpaquete, un contenedor para mcmcsamp)
• En AD Model Builder, posiblemente a través del glmmADMBpaquete (use la mcmc=TRUEopción) o el R2admbpaquete (escriba su propia definición de modelo en AD Model Builder), o fuera de R
• A través de la simfunción del armpaquete (simula el posterior solo para los coeficientes beta (efecto fijo)

Enfoques completamente bayesianos:

• A través del MCMCglmmpaquete
• Uso glmmBUGS(una interfaz WinBUGS wrapper / R )
• Usando JAGS / WinBUGS / OpenBUGS, etc., a través de los paquetes rjags/ r2jags/ R2WinBUGS/BRugs

En aras de la ilustración para mostrar cómo se vería esto, a continuación se MCMCglmmestima que utilizando el MCMCglmmpaquete, verá resultados similares que el modelo anterior y tiene algún tipo de valores p bayesianos:

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 comienza con modelos lineales mixtos y trata de estimarlos en R es leer las preguntas frecuentes de Wiki de donde se extrajo la mayor parte de esta información. Es un recurso excelente para todo tipo de temas de efectos mixtos, desde básico hasta avanzado y desde modelado hasta trazado.

Con respecto al cálculo de los valores de significación ( p ), Luke (2016) Evaluando la significación en modelos lineales de efectos mixtos en R informa que el método óptimo es la aproximación de Kenward-Roger o Satterthwaite para grados de libertad (disponible en R con paquetes como lmerTesto afex)

Abstracto

Los modelos de efectos mixtos se utilizan cada vez con más frecuencia en el análisis de datos experimentales. Sin embargo, en el paquete lme4 en R, los estándares para evaluar la importancia de los efectos fijos en estos modelos (es decir, obtener valores p) son algo vagos. Hay buenas razones para esto, pero como los investigadores que usan estos modelos están obligados en muchos casos a informar valores p, se necesita algún método para evaluar la importancia del resultado del modelo. Este artículo informa los resultados de las simulaciones que muestran que los dos métodos más comunes para evaluar la significancia, usando pruebas de razón de probabilidad y aplicando la distribución z a los valores t de Wald del resultado del modelo (t-as-z), son algo anticonservadores, especialmente para muestras más pequeñas. Otros métodos para evaluar la importancia,Los resultados de estas simulaciones sugieren que las tasas de error Tipo 1 son más cercanas a .05 cuando los modelos se ajustan usando REML y los valores p se derivan usando las aproximaciones de Kenward-Roger o Satterthwaite, ya que ambas aproximaciones produjeron tasas de error Tipo 1 aceptables incluso para los más pequeños. muestras

(énfasis añadido)

Yo uso el lmerTestpaquete. Esto incluye convenientemente una estimación del valor p en la anova()salida para mis análisis MLM, pero no da un tamaño del efecto por las razones dadas en otras publicaciones aquí.