每个 lmer 和方程的斜率的各个水平
Individual levels of slopes per lmer and equations
我想寻求一些帮助来描述 lmer() 模型生成的斜率。
我掌握的数据是不同老鼠在不同日子的质量体积。每只老鼠测量该体积的时间点不同。
对于大鼠 1，我在第 c(89,110,117,124,131,138) 天的体积 c(78,304,352,690,952,1250) 属于智利国家
对于大鼠 2，我在第 c(75,89,96,103,110) 天有卷 c(202,440,520,870,1380)，属于智利国家。
对于老鼠 3，我在第 c(75,89,96,103,110) 天有卷 c(186,370,620,850,1150)，属于智利国家。
对于大鼠 4，我在第 c(47,61,75,82,89,97,103,110) 天的卷 c(92,250,430,450,510,850,1000,1200) 属于英格兰国家。
对于大鼠 5，我在第 c(47,61,75,82) 天的卷 c(110,510,710,1200) 属于英格兰国家。
对于大鼠 6，我在第 c(47,61,75,82,89,97,103,110) 天的卷 c(115,380,480,540,560,850,1150,1350) 属于英格兰国家。
lmer 模型是：
m1 <- lmer(lVolume ~ Country*Day + (1|Rat))
我使用以下方法绘制了模型的曲线：
m1%>%
augment() %>%
clean_names() %>%
ggplot(data = .,
mapping = aes(x = day,
y = exp(l_volume),
group = rat)) +
geom_point(alpha = 0.5) +
geom_line(alpha = 0.5) +
geom_point(aes(y = exp(fitted)),
color = "red") +
geom_line(aes(y = exp(fitted)),
color = "red") +
expand_limits(x = 0 , y = 0)
这个模型根据模型 m1 为我预测了全国每只老鼠的新数据点。
从这个 lmer() 我在整个测量中有一个斜率，这是：
并通过 exp（预测）：
但是，我想用不同的方式来绘制它。我想绘制由我拥有的国家/地区的每个级别生成的斜率。
红线是智利和英格兰生成的 exp(斜率)，但也描绘了包含两个水平的整个模型的 exp(斜率)。
所以，最初我认为创建三个 lmer() 模型：
m1 <- lmer(lVolume ~ Country*Day + (1|Rat))
m2 <- lmer(lVolume ~ Day + (1|Rat)) (Rats in Chile)
m3 <- lmer(lVolume ~ Day + (1|Rat)) (Rats in England)
但我注意到 m2 和 m3 是完全不同的模型，因为它们没有来自 Country 的交互，这是我想检查的。所以，我不知道该怎么办。
更新
我尝试了这个并且有点奏效：
Final.Fixed<-effect(c("Country*Day"), m1,
xlevels=list(Day=seq(0,168,14)))
Final.Fixed<-as.data.frame(Final.Fixed)
Final.Fixed.Plot <-ggplot(data = Final.Fixed, aes(x = Day, y =exp(fit), group=Country))+
coord_cartesian(xlim=c(0,170),ylim = c(0,8000))+
geom_line(aes(color=Country), size=2)+
geom_ribbon(aes(ymin=exp(fit-se), ymax=exp(fit+se),fill=Country),alpha=.2)+
xlab("Day")+
ylab("Volume")+
scale_color_manual(values=c("blue", "red"))+
scale_fill_manual(values=c("blue", "red"))+
theme_bw()+
theme(text=element_text(face="bold", size=12),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(fill = NA, colour = "NA"),
axis.line = element_line(size = 1, colour = "grey80"),
legend.title=element_blank(),
legend.position = c(.2, .92))
Final.Fixed.Plot
这样可以吗？我认为我仍在考虑 m1 与 country*Day 的互动。如果我有问题，请纠正我！另外，我不知道如何为整个模型和该图中的原始数据点添加 exp(fit) 曲线。
可以给我一些 hint/help 吗？
清理顶部摘要
第一个代码块包含一个清理版本，它使用评论中的一些输入解决了问题的所有要点。我把原来的答案留在下面，一步步构建到最终情节。
library(tidyverse)
library(lme4)
library(broom.mixed)
library(ggeffects)
m1 <- lme4::lmer(lVolume ~ Country*Day + (1|Rat), data = df_rats %>%
dplyr::mutate(lVolume = log(Volume)))
# predictions for each country
syn_df <- tidyr::expand_grid(
Day = 1:170,
Country = c("Chile", "England")
) %>%
dplyr::mutate(lVolume = predict(m1, ., re.form = ~0))
# marginal effects for variable "Day"
df_day_marginal <- ggeffect(model = m1, terms = "Day", type = "fe") %>%
as.data.frame() %>%
dplyr::rename(Day = x, lVolume = predicted) %>%
dplyr::mutate(Country = "overall")
#combine prediction curves
df_preds <- bind_rows(syn_df, df_day_marginal)
# manually assemble formulas [units missing]
y0 <- round(fixef(m1)[["(Intercept)"]], 2)
beta_day <- round(fixef(m1)[["Day"]], 3)
beta_englday <- round(fixef(m1)[["CountryEngland:Day"]], 3)
beta_engl <- round(fixef(m1)[["CountryEngland"]], 2)
f_chile <- paste0("volume = exp(", y0, " + ", beta_day, " * days)")
f_england <- paste0("volume = exp(", y0 + beta_engl , " + ", beta_day + beta_englday, " * days)")
df_labels <- data.frame(
x = c(50, 50),
y = c(1300, 1400),
form = c(f_chile, f_england),
country = c("Chile", "England")
)
m1 %>%
broom.mixed::augment()%>%
ggplot(aes(x = Day, y = exp(lVolume), color = Country)) +
geom_ribbon(data = df_preds, aes(ymin = exp(conf.low), ymax = exp(conf.high), color = NULL, fill = Country), alpha = 0.3) +
geom_line(data = df_preds, size = 1.5) +
geom_line(aes(group = Rat)) +
geom_point() +
coord_cartesian(ylim = c(0, 1500), xlim = c(0, 150)) +
geom_text(data = df_labels, aes(x = x, y = y, label = form, color = country)) +
labs(x = "days", y = "volume")
原回答
对于问题的第一部分，我已尝试尽可能接近您的初始代码。
第一个块训练模型并population-level 预测智利和英格兰在指定日期内的表现。 （使用 re.form = ~0 参数，如 所述）
library(tidyverse)
library(lme4)
library(broom.mixed)
#helpful to specify in that `lVolume` is the log of the data you provid in the question
m1 <- lme4::lmer(lVolume ~ Country*Day + (1|Rat), data = df_rats %>%
dplyr::mutate(lVolume = log(Volume)))
days <- seq(0,168,14)
syn_df <- tidyr::expand_grid(
Day = 1:170,
Country = c("Chile", "England")
)
syn_df <- syn_df %>%
dplyr::mutate(l_volume = predict(m1, syn_df, re.form = ~0)) %>%
janitor::clean_names()
然后可以将其添加到您的原始情节中，稍作修改：
m1 %>%
broom.mixed::augment() %>%
janitor::clean_names() %>%
ggplot(data = .,
mapping = aes(x = day,
y = exp(l_volume),
color = country)) +
geom_point(alpha = 0.7) +
geom_line(aes(group = rat), alpha = 0.7) +
expand_limits(x = 0 , y = 0) +
geom_line(data = syn_df, alpha = 1, size = 1.5) +
coord_cartesian(ylim = c(NA, 1500), xlim = c(NA, 150))
已添加
此外，我们可以在图中添加几天的边际效应。
df_day_marginal <- ggeffect(model = m1, terms = "Day", type = "fe")
m1 %>%
broom.mixed::augment() %>%
janitor::clean_names() %>%
ggplot() +
geom_ribbon(data = df_day_marginal, aes(x = x, ymin = exp(conf.low), ymax = exp(conf.high)), alpha = 0.3) +
geom_line(data = syn_df, aes(x = day, y = exp(l_volume), color = country), size = 1.5) +
geom_line(data = df_day_marginal, aes(x = x, y = exp(predicted)), size = 1.5) +
geom_point(aes(x = day, y = exp(l_volume), color = country), alpha = 0.7) +
geom_line(aes(x = day, y = exp(l_volume), color = country, group = rat), alpha = 0.7) +
expand_limits(x = 0 , y = 0) +
coord_cartesian(ylim = c(NA, 1500), xlim = c(NA, 150)) +
labs(x = "days", y = "volume")
我想寻求一些帮助来描述 lmer() 模型生成的斜率。
我掌握的数据是不同老鼠在不同日子的质量体积。每只老鼠测量该体积的时间点不同。
对于大鼠 1，我在第 c(89,110,117,124,131,138) 天的体积 c(78,304,352,690,952,1250) 属于智利国家
对于大鼠 2，我在第 c(75,89,96,103,110) 天有卷 c(202,440,520,870,1380)，属于智利国家。
对于老鼠 3，我在第 c(75,89,96,103,110) 天有卷 c(186,370,620,850,1150)，属于智利国家。
对于大鼠 4，我在第 c(47,61,75,82,89,97,103,110) 天的卷 c(92,250,430,450,510,850,1000,1200) 属于英格兰国家。
对于大鼠 5，我在第 c(47,61,75,82) 天的卷 c(110,510,710,1200) 属于英格兰国家。
对于大鼠 6，我在第 c(47,61,75,82,89,97,103,110) 天的卷 c(115,380,480,540,560,850,1150,1350) 属于英格兰国家。
lmer 模型是：
m1 <- lmer(lVolume ~ Country*Day + (1|Rat))
我使用以下方法绘制了模型的曲线：
m1%>%
augment() %>%
clean_names() %>%
ggplot(data = .,
mapping = aes(x = day,
y = exp(l_volume),
group = rat)) +
geom_point(alpha = 0.5) +
geom_line(alpha = 0.5) +
geom_point(aes(y = exp(fitted)),
color = "red") +
geom_line(aes(y = exp(fitted)),
color = "red") +
expand_limits(x = 0 , y = 0)
这个模型根据模型 m1 为我预测了全国每只老鼠的新数据点。
从这个 lmer() 我在整个测量中有一个斜率，这是：
并通过 exp（预测）：
但是，我想用不同的方式来绘制它。我想绘制由我拥有的国家/地区的每个级别生成的斜率。
红线是智利和英格兰生成的 exp(斜率)，但也描绘了包含两个水平的整个模型的 exp(斜率)。
所以，最初我认为创建三个 lmer() 模型：
m1 <- lmer(lVolume ~ Country*Day + (1|Rat))
m2 <- lmer(lVolume ~ Day + (1|Rat)) (Rats in Chile)
m3 <- lmer(lVolume ~ Day + (1|Rat)) (Rats in England)
但我注意到 m2 和 m3 是完全不同的模型，因为它们没有来自 Country 的交互，这是我想检查的。所以，我不知道该怎么办。
更新
我尝试了这个并且有点奏效：
Final.Fixed<-effect(c("Country*Day"), m1,
xlevels=list(Day=seq(0,168,14)))
Final.Fixed<-as.data.frame(Final.Fixed)
Final.Fixed.Plot <-ggplot(data = Final.Fixed, aes(x = Day, y =exp(fit), group=Country))+
coord_cartesian(xlim=c(0,170),ylim = c(0,8000))+
geom_line(aes(color=Country), size=2)+
geom_ribbon(aes(ymin=exp(fit-se), ymax=exp(fit+se),fill=Country),alpha=.2)+
xlab("Day")+
ylab("Volume")+
scale_color_manual(values=c("blue", "red"))+
scale_fill_manual(values=c("blue", "red"))+
theme_bw()+
theme(text=element_text(face="bold", size=12),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(fill = NA, colour = "NA"),
axis.line = element_line(size = 1, colour = "grey80"),
legend.title=element_blank(),
legend.position = c(.2, .92))
Final.Fixed.Plot
这样可以吗？我认为我仍在考虑 m1 与 country*Day 的互动。如果我有问题，请纠正我！另外，我不知道如何为整个模型和该图中的原始数据点添加 exp(fit) 曲线。
可以给我一些 hint/help 吗？
清理顶部摘要
第一个代码块包含一个清理版本，它使用评论中的一些输入解决了问题的所有要点。我把原来的答案留在下面，一步步构建到最终情节。
library(tidyverse)
library(lme4)
library(broom.mixed)
library(ggeffects)
m1 <- lme4::lmer(lVolume ~ Country*Day + (1|Rat), data = df_rats %>%
dplyr::mutate(lVolume = log(Volume)))
# predictions for each country
syn_df <- tidyr::expand_grid(
Day = 1:170,
Country = c("Chile", "England")
) %>%
dplyr::mutate(lVolume = predict(m1, ., re.form = ~0))
# marginal effects for variable "Day"
df_day_marginal <- ggeffect(model = m1, terms = "Day", type = "fe") %>%
as.data.frame() %>%
dplyr::rename(Day = x, lVolume = predicted) %>%
dplyr::mutate(Country = "overall")
#combine prediction curves
df_preds <- bind_rows(syn_df, df_day_marginal)
# manually assemble formulas [units missing]
y0 <- round(fixef(m1)[["(Intercept)"]], 2)
beta_day <- round(fixef(m1)[["Day"]], 3)
beta_englday <- round(fixef(m1)[["CountryEngland:Day"]], 3)
beta_engl <- round(fixef(m1)[["CountryEngland"]], 2)
f_chile <- paste0("volume = exp(", y0, " + ", beta_day, " * days)")
f_england <- paste0("volume = exp(", y0 + beta_engl , " + ", beta_day + beta_englday, " * days)")
df_labels <- data.frame(
x = c(50, 50),
y = c(1300, 1400),
form = c(f_chile, f_england),
country = c("Chile", "England")
)
m1 %>%
broom.mixed::augment()%>%
ggplot(aes(x = Day, y = exp(lVolume), color = Country)) +
geom_ribbon(data = df_preds, aes(ymin = exp(conf.low), ymax = exp(conf.high), color = NULL, fill = Country), alpha = 0.3) +
geom_line(data = df_preds, size = 1.5) +
geom_line(aes(group = Rat)) +
geom_point() +
coord_cartesian(ylim = c(0, 1500), xlim = c(0, 150)) +
geom_text(data = df_labels, aes(x = x, y = y, label = form, color = country)) +
labs(x = "days", y = "volume")
原回答
对于问题的第一部分，我已尝试尽可能接近您的初始代码。
第一个块训练模型并population-level 预测智利和英格兰在指定日期内的表现。 （使用 re.form = ~0 参数，如
library(tidyverse)
library(lme4)
library(broom.mixed)
#helpful to specify in that `lVolume` is the log of the data you provid in the question
m1 <- lme4::lmer(lVolume ~ Country*Day + (1|Rat), data = df_rats %>%
dplyr::mutate(lVolume = log(Volume)))
days <- seq(0,168,14)
syn_df <- tidyr::expand_grid(
Day = 1:170,
Country = c("Chile", "England")
)
syn_df <- syn_df %>%
dplyr::mutate(l_volume = predict(m1, syn_df, re.form = ~0)) %>%
janitor::clean_names()
然后可以将其添加到您的原始情节中，稍作修改：
m1 %>%
broom.mixed::augment() %>%
janitor::clean_names() %>%
ggplot(data = .,
mapping = aes(x = day,
y = exp(l_volume),
color = country)) +
geom_point(alpha = 0.7) +
geom_line(aes(group = rat), alpha = 0.7) +
expand_limits(x = 0 , y = 0) +
geom_line(data = syn_df, alpha = 1, size = 1.5) +
coord_cartesian(ylim = c(NA, 1500), xlim = c(NA, 150))
已添加
此外，我们可以在图中添加几天的边际效应。
df_day_marginal <- ggeffect(model = m1, terms = "Day", type = "fe")
m1 %>%
broom.mixed::augment() %>%
janitor::clean_names() %>%
ggplot() +
geom_ribbon(data = df_day_marginal, aes(x = x, ymin = exp(conf.low), ymax = exp(conf.high)), alpha = 0.3) +
geom_line(data = syn_df, aes(x = day, y = exp(l_volume), color = country), size = 1.5) +
geom_line(data = df_day_marginal, aes(x = x, y = exp(predicted)), size = 1.5) +
geom_point(aes(x = day, y = exp(l_volume), color = country), alpha = 0.7) +
geom_line(aes(x = day, y = exp(l_volume), color = country, group = rat), alpha = 0.7) +
expand_limits(x = 0 , y = 0) +
coord_cartesian(ylim = c(NA, 1500), xlim = c(NA, 150)) +
labs(x = "days", y = "volume")