fabuła efektów mieszanych model w ggplot

Question

Jestem nowy z mieszanym efektem i potrzebuję twojej pomocy proszę. Mam kreślone na wykresie poniżej w 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() 

Jednak chciałbym reprezentować mieszane efekty modelu zamiast lm w geom_smooth, więc mogę to SITE jako efekt losowy.

model byłby następujący:

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

mam włączone TRTYEAR (rok leczenia), ponieważ jestem również zainteresowany wzorców efekt, który może zwiększać lub zmniejszać w czasie dla niektórych grup.

Dalej jest moja najlepsza próba tak daleko, aby wyodrębnić zmienne kreślenia z modelem, ale tylko wyodrębnione wartości TRTYEAR = 5, 10 i 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 

sugestie do zupełnie innego podejścia do reprezentowania ta analiza jest mile widziana. Pomyślałem, że to pytanie lepiej pasuje do stackoverflow, ponieważ chodzi o techniczne szczegóły w R zamiast statystyk. Dzięki

Answer 1

Możesz reprezentować swój model na wiele różnych sposobów. Najłatwiej jest narysować dane za pomocą różnych parametrów za pomocą różnych narzędzi do kreślenia (kolor, kształt, typ linii, aspekt), co zostało zrobione w przykładzie z wyjątkiem efektu losowego strona. Resztki modelu można również nanieść w celu przekazania wyników. Jak komentarz @MrFlick skomentował, to zależy od tego, co chcesz komunikować. Jeśli chcesz dodać przedziały ufności/przewidywania wokół swoich oszacowań, będziesz musiał głębiej kopać i rozważyć większe problemy statystyczne (example1, example2).

Oto przykład nieco dalej.
Ponadto w swoim komentarzu podałeś, że nie dostarczyłeś odtwarzalnego przykładu, ponieważ dane nie należą do Ciebie. To nie znaczy, że nie możesz podać przykładu z wymyślonych danych. Zastanów się, że w przypadku przyszłych postów możesz uzyskać szybsze odpowiedzi.

#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 

Nawiasem mówiąc, model pasuje do danych również w stosunku do współczynników powyżej:

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  

Adaptacja przykład pokazanie wyników modelu nałożone na dane

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

uwaga Wszystkie I zmieniono kolor z Myc na strona, a rodzaj linii na Myc.

Mam nadzieję, że ten przykład daje kilka pomysłów, jak zobrazować swój mieszany model efektów.