在ggplot中绘制混合效果模型 [英] plot mixed effects model in 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 ，所以我可以包含 SITE 作为随机效应。
该模型如下：
<$ p $ lt; code> library（lme4）
tempEf $ TRTYEAR< - as.numeric（tempEf $ TRTYEAR）
mod <-lmer（r_Myc * N * TRTYEAR +（ 1 | SITE），data = tempEf）
我已经包含 TRTYEAR （治疗年），因为我也对效果的模式感兴趣，可能会增加或减少随着时间的推移，一些团体。
接下来是我迄今为止最好的尝试从模型中提取绘图变量，但只提取 TRTYEAR = 5,10和15。
库（特效）
ef< ; - 效果（Myc：N：TRTYEAR，mod）
x< - as.data.frame（ef）
> x
Myc N TRTYEAR适合低位
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
对完全不同的repr方法的建议如有这种分析，欢迎。我认为这个问题更适合于stackoverflow，因为它是关于R的技术而不是背后的统计。谢谢
您可以用各种不同的方式来表示您的模型。最简单的方法是使用不同的绘图工具（颜色，形状，线型，刻面）通过各种参数绘制数据，除了随机效果站点之外，这是您对示例所做的操作。模型残差也可以被绘制以传达结果。就像@MrFlick评论的那样，它取决于你想要沟通的东西。如果你想围绕你的估计值增加置信度/预测带，你就必须深入挖掘并考虑更大的统计问题（ example1 ， example2 ）。
以下是一个仅以你为例进一步。
另外，在你的评论中，你说你没有提供一个可重复的例子，因为数据不属于你。这并不意味着你无法提供一个数据构成的例子。
＃制作数据：
tempEf< - data.frame（
N = rep（c（Nlow，Nhigh），each = 300），
Myc = rep（c（AM，ECM）， ，时间= 2），
TRTYEAR = runif（600,1,15），
site = rep（c（A，B，C，D，E ），每个= 10，次数= 12）＃5站点
）
＃组成一些响应数据
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截取;截距增加1
tempEf $ TRTYEAR / 10 * rnorm（600，mean = 0，sd = 2）＃增加一些噪音
library（lme4）
model< - lmer（r〜Myc * N * TRTYEAR +（1 | site），data = tempEf）
tempEf $ fit< - 预测（模型）＃增加模型拟合数据框
model
＃线性混合模型适合REML ['lmerMod']
#Formula：r〜Myc * N * TRTYEAR +（1 | site）
＃Data：tempEf
收敛时的#REML标准：2461.705
＃随机效果：
＃组名称Std.Dev。
＃网站（拦截）1.684
＃残值1.825
＃观看次数：600，群组：网站，5
＃固定效果：
＃（截距）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）+
me_bw（）
注意我所做的只是将您的颜色从 Myc 更改为站点，线型为 Myc 。
I希望这个例子给出了一些想法，如何可视化您的混合效果模型。
I am new with mixed effect models and I need your help please.
I have plotted the below graph in 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()
However, I would like to represent a mixed effects model instead of lmin geom_smooth, so I can include SITEas a random effect.
The model would be the following:
library(lme4)
tempEf$TRTYEAR <- as.numeric(tempEf$TRTYEAR)
mod <- lmer(r ~ Myc * N * TRTYEAR + (1|SITE), data=tempEf)
I have included TRTYEAR(year of treatment) because I am also interested in the patterns of the effect, that may increase or decrease over time for some groups.
Next is my best attempt so far to extract the plotting variables out of the model, but only extracted the values for TRTYEAR= 5, 10 and 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
Suggestions to a completely different approach to represent this analysis are welcome. I thought this question is better suited for stackoverflow because it’s about the technicalities in R rather than the statistics behind. Thanks
解决方案 You can represent your model a variety of different ways. The easiest is to plot data by the various parameters using different plotting tools (color, shape, line type, facet), which is what you did with your example except for the random effect site. Model residuals can also be plotted to communicate results. Like @MrFlick commented, it depends on what you want to communicate. If you want to add confidence/prediction bands around your estimates, you'll have to dig deeper and consider bigger statistical issues (example1, example2).
Here's an example taking yours just a bit further.
Also, in your comment you said you didn't provide a reproducible example because the data do not belong to you. That doesn't mean you can't provide an example out of made up data. Please consider that for future posts so you can get faster answers.
#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
Incidentally, the model fit the data well compared to the coefficients above:
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
Adapting your example to show the model outputs overlaid on the data
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()
Notice all I did was change your colour from Myc to site, and linetype to Myc.
I hope this example gives some ideas how to visualize your mixed effects model.
这篇关于在ggplot中绘制混合效果模型的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
