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

введите описание изображения здесь

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


r ggplot2 lmer
person fede_luppi    schedule 26.06.2015    source источник
comment
Если у вас есть такой случайный эффект, вы больше не получите красивых простых линий. Как вы ожидаете увидеть сюжет? Кроме того, при обращении за помощью по программированию вы должны включить воспроизводимый пример, в котором есть образцы входных данных, чтобы мы могли запустить ваш код и протестировать возможные решения.   -  person MrFlick    schedule 26.06.2015
comment
Спасибо @MrFlick. Возможно, я ожидал бы построить КЭ, но у меня нет опыта, поэтому я не знаю, каким может быть ожидаемый результат с точки зрения графика. Что касается данных, я хотел точно представить проблемы и тип анализа, который мне нужен, но, конечно, настоящие данные мне не принадлежат, поэтому мне не разрешено публиковать их в Интернете.   -  person fede_luppi    schedule 26.06.2015
comment
@MrFlick Для публикации, не могли бы вы предложить использовать график, аналогичный приведенному выше, с lm для его визуализации и использовать lmer для статистического анализа?   -  person fede_luppi    schedule 26.06.2015

Ответы (1)


arrow_upward
31
arrow_downward

Вы можете представить свою модель разными способами. Самый простой - это построить график данных по различным параметрам с использованием различных инструментов построения графиков (цвет, форма, тип линии, фасет), что вы и сделали в своем примере, за исключением случайного эффекта site. Остатки модели также могут быть нанесены на график для представления результатов. Как прокомментировал @MrFlick, это зависит от того, о чем вы хотите сообщить. Если вы хотите добавить к своим оценкам доверительные интервалы / диапазоны прогнозов, вам придется копнуть глубже и рассмотреть более серьезные статистические проблемы (example1, пример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 на site и тип линии на Myc. lmer со случайными эффектами

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

person oshun    schedule 28.06.2015
comment
Спасибо за ответ. Ваш исчерпывающий ответ позволил мне понять различные потенциальные результаты анализа и то, что мне действительно нужно. Спасибо - person fede_luppi; 29.06.2015