studio - parcela modelo de efectos mixtos en ggplot

ggplot2 tutorial

Question

Soy nuevo con modelos de efectos mixtos y necesito su ayuda por favor. He trazado el siguiente gráfico en ggplot:

ggplot(tempEf,aes(TRTYEAR,CO2effect,group=Myc,col=Myc)) + 
  facet_grid(~N) +
  geom_smooth(method="lm",se=T,size=1) +
  geom_point(alpha = 0.3) + 
  geom_hline(yintercept=0, linetype="dashed") +
  theme_bw()

Sin embargo, me gustaría representar un modelo de efectos mixtos en lugar de lm en geom_smooth , así que puedo incluir SITE como un efecto aleatorio.

El modelo sería el siguiente:

library(lme4)
tempEf$TRTYEAR <- as.numeric(tempEf$TRTYEAR)
mod <- lmer(r ~ Myc * N * TRTYEAR + (1|SITE), data=tempEf)

He incluido TRTYEAR (año de tratamiento) porque también estoy interesado en los patrones del efecto, que pueden aumentar o disminuir con el tiempo en algunos grupos.

El siguiente es mi mejor intento hasta ahora para extraer las variables de trazado del modelo, pero solo extrajo los valores para TRTYEAR = 5, 10 y 15.

library(effects)
ef <- effect("Myc:N:TRTYEAR", mod)
x <- as.data.frame(ef)
> x
   Myc     N TRTYEAR        fit         se       lower     upper
1   AM  Nlow       5 0.04100963 0.04049789 -0.03854476 0.1205640
2  ECM  Nlow       5 0.41727928 0.07342289  0.27304676 0.5615118
3   AM Nhigh       5 0.20562700 0.04060572  0.12586080 0.2853932
4  ECM Nhigh       5 0.24754017 0.27647151 -0.29556267 0.7906430
5   AM  Nlow      10 0.08913042 0.03751783  0.01543008 0.1628307
6  ECM  Nlow      10 0.42211957 0.15631679  0.11504963 0.7291895
7   AM Nhigh      10 0.30411129 0.03615213  0.23309376 0.3751288
8  ECM Nhigh      10 0.29540744 0.76966410 -1.21652689 1.8073418
9   AM  Nlow      15 0.13725120 0.06325159  0.01299927 0.2615031
10 ECM  Nlow      15 0.42695986 0.27301163 -0.10934636 0.9632661
11  AM Nhigh      15 0.40259559 0.05990085  0.28492587 0.5202653
12 ECM Nhigh      15 0.34327471 1.29676632 -2.20410343 2.8906529

Las sugerencias para un enfoque completamente diferente para representar este análisis son bienvenidas. Pensé que esta pregunta es más adecuada para el stackoverflow porque se trata de los tecnicismos en R en lugar de las estadísticas detrás. Gracias

Answer 1

Puedes representar tu modelo de una variedad de maneras diferentes. Lo más fácil es trazar los datos por los diversos parámetros utilizando diferentes herramientas de trazado (color, forma, tipo de línea, faceta), que es lo que hizo con su ejemplo, excepto el sitio de efectos aleatorios. Los residuos del modelo también se pueden trazar para comunicar los resultados. Como comentó @MrFlick, depende de lo que quieras comunicar. Si desea agregar bandas de confianza / predicción a sus estimaciones, tendrá que profundizar más y considerar problemas estadísticos más grandes ( example1 , example2 ).

Aquí hay un ejemplo llevando el tuyo un poco más lejos.
Además, en su comentario usted dijo que no proporcionó un ejemplo reproducible porque los datos no le pertenecen. Eso no significa que no pueda proporcionar un ejemplo de los datos inventados. Por favor, tenlo en cuenta para futuras publicaciones para que puedas obtener respuestas más rápidas

#Make up data:
tempEf <- data.frame(
  N = rep(c("Nlow", "Nhigh"), each=300),
  Myc = rep(c("AM", "ECM"), each=150, times=2),
  TRTYEAR = runif(600, 1, 15),
  site = rep(c("A","B","C","D","E"), each=10, times=12)   #5 sites
  )

# Make up some response data
tempEf$r <- 2*tempEf$TRTYEAR +                   
            -8*as.numeric(tempEf$Myc=="ECM") +
            4*as.numeric(tempEf$N=="Nlow") +
            0.1*tempEf$TRTYEAR * as.numeric(tempEf$N=="Nlow") +
            0.2*tempEf$TRTYEAR*as.numeric(tempEf$Myc=="ECM") +
           -11*as.numeric(tempEf$Myc=="ECM")*as.numeric(tempEf$N=="Nlow")+ 
            0.5*tempEf$TRTYEAR*as.numeric(tempEf$Myc=="ECM")*as.numeric(tempEf$N=="Nlow")+ 
           as.numeric(tempEf$site) +  #Random intercepts; intercepts will increase by 1
           tempEf$TRTYEAR/10*rnorm(600, mean=0, sd=2)    #Add some noise

library(lme4)
model <- lmer(r ~ Myc * N * TRTYEAR + (1|site), data=tempEf)
tempEf$fit <- predict(model)   #Add model fits to dataframe

Por cierto, el modelo se ajusta bien a los datos en comparación con los coeficientes anteriores:

model

#Linear mixed model fit by REML [''lmerMod'']
#Formula: r ~ Myc * N * TRTYEAR + (1 | site)
#   Data: tempEf
#REML criterion at convergence: 2461.705
#Random effects:
# Groups   Name        Std.Dev.
# site     (Intercept) 1.684   
# Residual             1.825   
#Number of obs: 600, groups:  site, 5
#Fixed Effects:
#         (Intercept)                MycECM                 NNlow               
#             3.03411              -7.92755               4.30380               
#             TRTYEAR          MycECM:NNlow        MycECM:TRTYEAR  
#             1.98889             -11.64218               0.18589  
#       NNlow:TRTYEAR  MycECM:NNlow:TRTYEAR  
#             0.07781               0.60224      

Adaptar su ejemplo para mostrar los resultados del modelo superpuestos en los datos

   library(ggplot2)
    ggplot(tempEf,aes(TRTYEAR, r, group=interaction(site, Myc), col=site, shape=Myc )) + 
      facet_grid(~N) +
      geom_line(aes(y=fit, lty=Myc), size=0.8) +
      geom_point(alpha = 0.3) + 
      geom_hline(yintercept=0, linetype="dashed") +
      theme_bw()

Note que todo lo que hice fue cambiar su color de Myc al sitio , y el tipo de línea a Myc .

Espero que este ejemplo le dé algunas ideas sobre cómo visualizar su modelo de efectos mixtos.