modello di effetti a trama mista in ggplot

Question

Sono nuovo con modelli a effetti misti e ho bisogno del tuo aiuto per favore. ho tracciato il grafico di seguito nella 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() 

Tuttavia, vorrei rappresentare un modello misto effetti invece di lm in geom_smooth, così posso comprendere SITE come un effetto casuale.

Il modello sarebbe il seguente:

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

Ho incluso TRTYEAR (anno di trattamento) perché mi interessa nei modelli di effetti, che possono aumentare o diminuire nel tempo per alcuni gruppi anche.

successivo è il mio miglior tentativo finora di estrarre le variabili di tracciato fuori del modello, ma estratto solo i valori per TRTYEAR = 5, 10 e 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 

Suggerimenti per un approccio completamente diverso per rappresentare questa analisi è gradita. Ho pensato che questa domanda fosse più adatta allo stackoverflow perché riguarda i tecnicismi di R piuttosto che le statistiche alla base. Grazie

Answer 1

È possibile rappresentare il modello in vari modi. Il modo più semplice è di tracciare i dati in base ai vari parametri usando diversi strumenti di tracciamento (colore, forma, tipo di linea, sfaccettatura), che è quello che hai fatto con il tuo esempio eccetto l'effetto casuale sito. I residui del modello possono anche essere tracciati per comunicare i risultati. Come commentato @MrFlick, dipende da cosa vuoi comunicare. Se si desidera aggiungere bande di confidenza/previsione intorno alle stime, è necessario scavare più a fondo e considerare problemi statistici più grandi (example1, example2).

Ecco un esempio che ti accompagna un po 'oltre.
Inoltre, nel tuo commento hai detto che non hai fornito un esempio riproducibile perché i dati non ti appartengono. Ciò non significa che non puoi fornire un esempio di dati inventati. Si prega di considerare che per i post futuri in modo da poter ottenere risposte più veloci.

#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 

Per inciso, il modello si adatta bene ai dati rispetto ai coefficienti di cui sopra:

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  

Adattamento tuo esempio per mostrare gli output del modello sovrapposti sui dati

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() 

Avviso tutto quello che è stato modificato il colore da Myc a sito e tipo di linea a Myc.

Spero che questo esempio dia alcune idee su come visualizzare il modello a effetti misti.