Модель смешанных эффектов в ggplot

Я новичок в смешанных моделях, и мне нужна ваша помощь. Я построил график ниже в 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()

enter image description here

Однако я хотел бы представить модель смешанных эффектов вместо lm в geom_smooth, поэтому я могу включить SITE как случайный эффект.

Модель будет следующей:

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

Я включил TRTYEAR (год лечения), потому что меня также интересуют шаблоны эффекта, которые могут увеличиваться или уменьшаться со временем для некоторых групп.

Далее моя лучшая попытка до сих пор извлекать переменные построения из модели, но только извлекала значения для TRTYEAR= 5, 10 и 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

Приветствуются предложения относительно совершенно другого подхода к представлению этого анализа. Я подумал, что этот вопрос лучше подходит для stackoverflow, потому что его технические характеристики в R, а не статистика. Благодаря

Ответы

Ответ 1

Вы можете представить свою модель различными способами. Самый простой способ состоит в том, чтобы отображать данные по различным параметрам с использованием разных инструментов построения (цвет, форма, тип линии, фасет), что и было сделано с вашим примером, за исключением сайта случайного эффекта. Модификации модели также могут быть построены для передачи результатов. Как и комментарий @MrFlick, это зависит от того, что вы хотите сообщить. Если вы хотите добавить группы доверия/прогнозов вокруг своих оценок, вам придется глубже копать и рассматривать более крупные статистические проблемы (example1, example2).

Вот пример, в котором вы немного больше.
Кроме того, в своем комментарии вы сказали, что не предоставили воспроизводимый пример, потому что данные не принадлежат вам. Это не значит, что вы не можете представить пример из сделанных данных. Пожалуйста, учтите, что для будущих должностей вы можете получить более быстрые ответы.

#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

Кстати, модель хорошо подходит для данных по сравнению с приведенными выше коэффициентами:

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

Адаптация вашего примера для отображения выходов модели, наложенных на данные

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

Обратите внимание, что все, что я сделал, это изменить цвет с Myc на сайт и тип linceype на Myc. lmer with random effects

Надеюсь, что этот пример дает некоторые идеи о том, как визуализировать вашу модель смешанных эффектов.