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.