plot mixed effects modell in ggplot

Question

Ich bin neu mit Mixed-Effekt-Modelle und ich brauche Ihre Hilfe bitte. Ich habe die unten Graph in ggplot aufgetragen:

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

Allerdings würde ich statt lm in geom_smooth ein gemischtes Effects-Modell darstellen mag, so kann ich SITE als Zufallseffekt enthalten.

Das Modell wäre die folgende:

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

Ich habe TRTYEAR (Behandlungsjahr) aufgenommen, weil ich auch in den Mustern der Wirkung interessiert bin, die für einige Gruppen im Laufe der Zeit steigen oder fallen kann.

Weiter ist mein bester Versuch so weit die Plotten Variablen aus dem Modell zu extrahieren, aber nur extrahiert die Werte für TRTYEAR = 5, 10 und 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 

Vorschläge zu einem völlig anderen Ansatz zu vertreten Diese Analyse ist willkommen. Ich dachte, dass diese Frage besser für stackoverflow geeignet ist, da es sich um die technischen Details in R handelt und nicht um die Statistik dahinter. Danke

Answer 1

Sie können Ihr Modell auf verschiedene Arten darstellen. Am einfachsten ist es, die Daten anhand der verschiedenen Parameter mit verschiedenen Plotwerkzeugen (Farbe, Form, Linientyp, Facette) zu plotten, was Sie mit Ihrem Beispiel gemacht haben, mit Ausnahme des Random-Effekts site. Modellresiduen können auch grafisch dargestellt werden, um Ergebnisse zu kommunizieren. Wie @MrFlick kommentiert, hängt es davon ab, was Sie kommunizieren möchten. Wenn Sie Vertrauens-/Vorhersagebereiche um Ihre Schätzungen herum hinzufügen möchten, müssen Sie tiefer einsteigen und größere statistische Probleme in Betracht ziehen (example1, example2).

Hier ist ein Beispiel, das Ihnen nur ein bisschen weiter geht.
In Ihrem Kommentar haben Sie auch gesagt, dass Sie kein reproduzierbares Beispiel angegeben haben, weil die Daten nicht zu Ihnen gehören. Das bedeutet nicht, dass Sie kein Beispiel aus erfundenen Daten liefern können. Bitte beachte das für zukünftige Beiträge, damit du schnellere Antworten bekommst.

#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 

Übrigens passen das Modell auch die Daten im Vergleich zu den Koeffizienten oben:

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  

Ihr Beispiel Anpassung der Modellergebnisse auf den Daten

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

Hinweis alles, was ich überlagert zu zeigen, Did war Ihre Farbe von Myc zu Website, und Linientyp zu Myc.

Ich hoffe, dieses Beispiel gibt einige Ideen, wie Sie Ihr Mixed-Effects-Modell visualisieren können.