Построение модели смешанных эффектов в ggplot

Question

Построение модели смешанных эффектов в 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()

Тем не менее, я хотел бы представить модель смешанных эффектов вместо 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

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

27

r ggplot2 lmer

Источник

user2122078 26 июн '15 в 14:10

1 ответ

Решение

Другие вопросы по тегам r ggplot2 lmer

user4718512 28 июн '15 в 01:30 2015-06-28 01:30 · Accepted Answer · 2015-06-28 01:30

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

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

#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 на сайт и тип линии на Myc. Имер со случайными эффектами

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