在ggplot中绘制具有固定效应相互作用的混合效应模型
plot mixed effect model with interaction in fixed effects in ggplot
我正在绘制基于 lmer()
对象的混合效应模型中固定效应的相互作用。为此,我根据我的模型预测新值。这工作正常,除了由于我生成它们的方式,预测延伸到整个可能的 x 轴范围。我现在可以通过基于循环定义 new.dat(根据分组变量 "Variety" 更改最大值和最小值)等,将预测的回归线限制在各自分组变量的范围内,但是 - 是有更优雅/更简单的解决方案来绘制这个吗?我是否错过了什么(我对 R 比较陌生)?
数据:
library(datasets)
data("Oats")
# manipulate data so it resembles more my actual data
Oats <- Oats %>%
filter((Variety == "Golden Rain" & nitro>=0.2) | (Variety == "Marvellous" & nitro <=0.4) | (Variety == "Victory" & nitro<=0.4 & nitro>=0.2)) #%>%
模型和绘图:
mod2 <- lmer(yield ~ nitro * Variety + (1| Variety), data=Oats)
new.dat <- data.frame(nitro=seq(min(Oats$nitro),max(Oats$nitro), length.out = 48), Variety= Oats$Variety)
new.dat$pred<-predict(mod2,newdata=new.dat,re.form=~0)
ggplot(data=Oats, aes(x=nitro, y=yield, col = Variety)) +
geom_point() +
geom_line(data=new.dat, aes(y=pred)) +
geom_point(data=new.dat, aes(y=pred))
非常感谢您的每一个提示!
您可以通过计算每个组的 min/max 然后按组计算序列来获得它。继续使用 tidyverse,因为您的代码已经在使用它:
library(tidyverse)
library(pairwiseCI)
#> Loading required package: MCPAN
#> Loading required package: coin
#> Loading required package: survival
library(lme4)
#> Loading required package: Matrix
#>
#> Attaching package: 'Matrix'
#> The following objects are masked from 'package:tidyr':
#>
#> expand, pack, unpack
data("Oats")
## manipulate data so it resembles more my actual data
Oats <-
Oats %>%
filter((Variety == "Golden Rain" & nitro>=0.2) | (Variety == "Marvellous" & nitro <=0.4) | (Variety == "Victory" & nitro<=0.4 & nitro>=0.2)) #%>%
mod2 <- lmer(yield ~ nitro * Variety + (1| Variety), data=Oats)
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
#> unable to evaluate scaled gradient
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
#> Hessian is numerically singular: parameters are not uniquely determined
## Calculate min/max by group
all_vals <-
Oats %>%
group_by(Variety) %>%
summarize(min_nitro = min(nitro),
max_nitro = max(nitro))
## Calculate sequence for each group
new.dat <-
all_vals %>%
group_split(Variety) %>%
map_dfr(~ data.frame(Variety = .x$Variety, nitro = seq(.x$min_nitro, .x$max_nitro, length.out = 20)))
new.dat$pred<-predict(mod2,newdata=new.dat,re.form=~0)
ggplot(data=Oats, aes(x=nitro, y=yield, col = Variety)) +
geom_point() +
geom_line(data=new.dat, aes(y=pred)) +
geom_point(data=new.dat, aes(y=pred))
我正在绘制基于 lmer()
对象的混合效应模型中固定效应的相互作用。为此,我根据我的模型预测新值。这工作正常,除了由于我生成它们的方式,预测延伸到整个可能的 x 轴范围。我现在可以通过基于循环定义 new.dat(根据分组变量 "Variety" 更改最大值和最小值)等,将预测的回归线限制在各自分组变量的范围内,但是 - 是有更优雅/更简单的解决方案来绘制这个吗?我是否错过了什么(我对 R 比较陌生)?
数据:
library(datasets)
data("Oats")
# manipulate data so it resembles more my actual data
Oats <- Oats %>%
filter((Variety == "Golden Rain" & nitro>=0.2) | (Variety == "Marvellous" & nitro <=0.4) | (Variety == "Victory" & nitro<=0.4 & nitro>=0.2)) #%>%
模型和绘图:
mod2 <- lmer(yield ~ nitro * Variety + (1| Variety), data=Oats)
new.dat <- data.frame(nitro=seq(min(Oats$nitro),max(Oats$nitro), length.out = 48), Variety= Oats$Variety)
new.dat$pred<-predict(mod2,newdata=new.dat,re.form=~0)
ggplot(data=Oats, aes(x=nitro, y=yield, col = Variety)) +
geom_point() +
geom_line(data=new.dat, aes(y=pred)) +
geom_point(data=new.dat, aes(y=pred))
非常感谢您的每一个提示!
您可以通过计算每个组的 min/max 然后按组计算序列来获得它。继续使用 tidyverse,因为您的代码已经在使用它:
library(tidyverse)
library(pairwiseCI)
#> Loading required package: MCPAN
#> Loading required package: coin
#> Loading required package: survival
library(lme4)
#> Loading required package: Matrix
#>
#> Attaching package: 'Matrix'
#> The following objects are masked from 'package:tidyr':
#>
#> expand, pack, unpack
data("Oats")
## manipulate data so it resembles more my actual data
Oats <-
Oats %>%
filter((Variety == "Golden Rain" & nitro>=0.2) | (Variety == "Marvellous" & nitro <=0.4) | (Variety == "Victory" & nitro<=0.4 & nitro>=0.2)) #%>%
mod2 <- lmer(yield ~ nitro * Variety + (1| Variety), data=Oats)
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
#> unable to evaluate scaled gradient
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
#> Hessian is numerically singular: parameters are not uniquely determined
## Calculate min/max by group
all_vals <-
Oats %>%
group_by(Variety) %>%
summarize(min_nitro = min(nitro),
max_nitro = max(nitro))
## Calculate sequence for each group
new.dat <-
all_vals %>%
group_split(Variety) %>%
map_dfr(~ data.frame(Variety = .x$Variety, nitro = seq(.x$min_nitro, .x$max_nitro, length.out = 20)))
new.dat$pred<-predict(mod2,newdata=new.dat,re.form=~0)
ggplot(data=Oats, aes(x=nitro, y=yield, col = Variety)) +
geom_point() +
geom_line(data=new.dat, aes(y=pred)) +
geom_point(data=new.dat, aes(y=pred))