我试图仅绘制来自lmer对象的随机效应的几个术语。我将借用oshun here发布的数据集。
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
)
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)
plot_model(model, type = "re")
plot_model(model, type = "re", terms = c("A", "E"))
由于某些原因,术语选项不起作用。应该反映给sjplot的作者。这是两种解决方法:
1)使用ggplot2手动定义术语:
plot_model(model, type = "re") + scale_x_discrete(limits=c("A","D"))
由于抛出数据而收到以下警告
'x'的刻度已经存在。为“ x”添加另一个比例,将取代现有的规模。警告消息:1:已删除3行包含缺失值(geom_point)。 2:删除了3行,其中包含缺少值(geom_errorbar)。
2)从get_model()绘制它>
df = get_model_data(model,type="re")
df = subset(df,term %in% c("A","D"))
ggplot(df,aes(x=term,y=estimate,col=group)) +
geom_point(show.legend=FALSE) +
geom_segment(aes(xend=term,y=conf.low,yend=conf.high),show.legend=FALSE)+
scale_colour_brewer(palette = "Set1")+
coord_flip()+ggtitle("Random Effects")