emmip - 有没有办法绘制 'merged' 治疗 x 时间因子？
emmip - is there a way to plot a 'merged' treatment x time factor?
我附上了一些来自 link 的已发布数据（最后）：
https://www.middleprofessor.com/files/applied-biostatistics_bookdown/_book/pre-post.html
我可以 运行 一个相当标准的线性混合模型，包含治疗、时间和治疗 x 时间交互，并得到相同的结果：
# LDA (mixed model)----
> m3 <- lme4::lmer(y ~ treat*time + (1|id), data = dat)
> summary(m3)
Linear mixed model fit by REML ['lmerMod']
Formula: y ~ treat * time + (1 | id)
Data: dat
REML criterion at convergence: 263.3
Scaled residuals:
Min 1Q Median 3Q Max
-1.7804 -0.5552 0.1779 0.6036 1.2683
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 26.57 5.155
Residual 46.49 6.819
Number of obs: 40, groups: id, 20
Fixed effects:
Estimate Std. Error t value
(Intercept) 54.656 2.849 19.182
treatil_17 2.472 3.842 0.644
timepost -1.155 3.214 -0.359
treatil_17:timepost 11.712 4.334 2.702
Correlation of Fixed Effects:
(Intr) trt_17 timpst
treatil_17 -0.742
timepost -0.564 0.418
trtl_17:tmp 0.418 -0.564 -0.742
并使用 emmip 很好地绘制结果：
> emmeans(m3, revpairwise ~ treat | time)
$emmeans
time = pre:
treat emmean SE df lower.CL upper.CL
veh 54.7 2.85 31.8 48.9 60.5
il_17 57.1 2.58 31.8 51.9 62.4
time = post:
treat emmean SE df lower.CL upper.CL
veh 53.5 2.85 31.8 47.7 59.3
il_17 67.7 2.58 31.8 62.4 72.9
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
time = pre:
contrast estimate SE df t.ratio p.value
il_17 - veh 2.47 3.84 31.8 0.644 0.5245
time = post:
contrast estimate SE df t.ratio p.value
il_17 - veh 14.18 3.84 31.8 3.692 0.0008
Degrees-of-freedom method: kenward-roger
> emmip(m3, treat ~ time) +
+ geom_point() +
+ geom_line(size = 1) +
+ scale_y_continuous(limits = c(50, 70), breaks = seq(50, 70, by = 2)) +
+ theme_bw(base_size = 15)
geom_path: Each group consists of only one observation. Do you need to adjust the group aesthetic?
我现在想要 运行 模型的约束版本，基本上强制组的主效应等于 0。为此，我按照 'merge' 处理和时间因素与共同的基线水平相结合，然后是每个 post 治疗的新水平（第 12 页，以下 link）：
https://bendixcarstensen.com/Notes/rm.pdf
# CONSTRAINED LDA (mixed model) ----
# Create a 'merged' treat*time variable with 3 levels:
# - a common pre level for both treatments
# - a new post level for each treatment
dat <- transform(dat, treat_time = Epi::Relevel(interaction(treat, time), list(baseline = 1:2)))
# Inspect levels
with(dat, ftable(treat, time))
with(dat, ftable(treat_time, treat, time))
> m4 <- lme4::lmer(y ~ treat_time + (1|id), data = dat)
> summary(m4)
Linear mixed model fit by REML ['lmerMod']
Formula: y ~ treat_time + (1 | id)
Data: dat
REML criterion at convergence: 268.2
Scaled residuals:
Min 1Q Median 3Q Max
-1.7986 -0.5202 0.1392 0.6089 1.2810
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 25.74 5.074
Residual 46.02 6.784
Number of obs: 40, groups: id, 20
Fixed effects:
Estimate Std. Error t value
(Intercept) 56.016 1.894 29.571
treat_timeveh.post -2.027 2.902 -0.698
treat_timeil_17.post 11.270 2.676 4.212
Correlation of Fixed Effects:
(Intr) trt_t.
trt_tmvh.ps -0.419
trt_tml_17. -0.454 0.190
> emmeans(m4, revpairwise ~ treat_time)
$emmeans
treat_time emmean SE df lower.CL upper.CL
baseline 56.0 1.89 32.9 52.2 59.9
veh.post 54.0 2.79 37.0 48.3 59.6
il_17.post 67.3 2.53 36.9 62.2 72.4
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
veh.post - baseline -2.03 2.97 24.2 -0.684 0.7752
il_17.post - baseline 11.27 2.72 22.8 4.141 0.0011
il_17.post - veh.post 13.30 3.72 32.5 3.573 0.0031
Degrees-of-freedom method: kenward-roger
P value adjustment: tukey method for comparing a family of 3 estimates
> emmip(m4, ~ treat_time) +
+ geom_point() +
+ geom_line(size = 1) +
+ scale_y_continuous(limits = c(50, 70), breaks = seq(50, 70, by = 2)) +
+ theme_bw(base_size = 15)
geom_path: Each group consists of only one observation. Do you need to adjust the group aesthetic?
>
但我不知道如何（或者是否有办法）使用 emmip 生成类似的图，考虑到因子水平是组合的。我希望能够从本质上重现上面的第一个图，显示每个组从公共基线的变化。我想我可以在 ggplot 中手动执行此操作 - 只是想知道是否有某种方法可以直接使用 emmip 进行操作？
dat <- structure(list(id = c(1L, 1L, 2L, 2L, 3L, 3L, 4L, 4L, 5L, 5L,
6L, 6L, 7L, 7L, 8L, 8L, 9L, 9L, 10L, 10L, 11L, 11L, 12L, 12L,
13L, 13L, 14L, 14L, 15L, 15L, 16L, 16L, 17L, 17L, 18L, 18L, 19L,
19L, 20L, 20L), treat = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L), .Label = c("veh", "il_17"), class = "factor"), time = structure(c(1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L), .Label = c("pre", "post"), class = "factor"),
y = c(63.9871851, 62.7837203, 58.4825871, 59.4795929, 59.7735716,
42.4923719, 51.0754017, 55.5395164, 55.8808155, 42.8811745,
38.686015, 40.6969206, 46.8213372, 61.6, 62.6084263, 62.987013,
54.5901379, 53.0518311, 60.99353, 65.8177226, 66.2765671,
77.6594363, 60.1260258, 80.5687204, 49.444707, 63.7658228,
55.1138485, 77.6594363, 61.4070066, 80.5687204, 52.4932976,
49.444707, 55.1138485, 57.5908965, 61.4070066, 54.4305875,
52.4932976, 68.6955215, 53.5461576, 68.3395943), treat_time = structure(c(1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L,
1L, 2L, 1L, 3L, 1L, 3L, 1L, 3L, 1L, 3L, 1L, 3L, 1L, 3L, 1L,
3L, 1L, 3L, 1L, 3L, 1L, 3L, 1L, 3L), .Label = c("baseline",
"veh.post", "il_17.post"), class = "factor")), class = "data.frame", row.names = c(NA,
-40L))
我想如果你对每个治疗使用相同的基线值，你可以这样做：
foo = emmeans(m4, "treat_time")
bar = contrast(foo, list(b1 = c(1,0,0), p1 = c(0,1,0),
b2 = c(1,0,0), p2 = c(0,0,1)) )
levels(bar) = list(time = c("baseline", "post"),
treat = c("veh", "il_17") )
emmip(bar, treat ~ time)
更简单的方法
我发现你可以大大简化第二步：
foo = emmeans(m4, "treat_time")
bar = foo[c(1,2,1,3)]
levels(bar) = list(time = c("baseline", "post"),
treat = c("veh", "il_17") )
emmip(bar, treat ~ time)
我附上了一些来自 link 的已发布数据（最后）：
https://www.middleprofessor.com/files/applied-biostatistics_bookdown/_book/pre-post.html
我可以 运行 一个相当标准的线性混合模型，包含治疗、时间和治疗 x 时间交互，并得到相同的结果：
# LDA (mixed model)----
> m3 <- lme4::lmer(y ~ treat*time + (1|id), data = dat)
> summary(m3)
Linear mixed model fit by REML ['lmerMod']
Formula: y ~ treat * time + (1 | id)
Data: dat
REML criterion at convergence: 263.3
Scaled residuals:
Min 1Q Median 3Q Max
-1.7804 -0.5552 0.1779 0.6036 1.2683
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 26.57 5.155
Residual 46.49 6.819
Number of obs: 40, groups: id, 20
Fixed effects:
Estimate Std. Error t value
(Intercept) 54.656 2.849 19.182
treatil_17 2.472 3.842 0.644
timepost -1.155 3.214 -0.359
treatil_17:timepost 11.712 4.334 2.702
Correlation of Fixed Effects:
(Intr) trt_17 timpst
treatil_17 -0.742
timepost -0.564 0.418
trtl_17:tmp 0.418 -0.564 -0.742
并使用 emmip 很好地绘制结果：
> emmeans(m3, revpairwise ~ treat | time)
$emmeans
time = pre:
treat emmean SE df lower.CL upper.CL
veh 54.7 2.85 31.8 48.9 60.5
il_17 57.1 2.58 31.8 51.9 62.4
time = post:
treat emmean SE df lower.CL upper.CL
veh 53.5 2.85 31.8 47.7 59.3
il_17 67.7 2.58 31.8 62.4 72.9
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
time = pre:
contrast estimate SE df t.ratio p.value
il_17 - veh 2.47 3.84 31.8 0.644 0.5245
time = post:
contrast estimate SE df t.ratio p.value
il_17 - veh 14.18 3.84 31.8 3.692 0.0008
Degrees-of-freedom method: kenward-roger
> emmip(m3, treat ~ time) +
+ geom_point() +
+ geom_line(size = 1) +
+ scale_y_continuous(limits = c(50, 70), breaks = seq(50, 70, by = 2)) +
+ theme_bw(base_size = 15)
geom_path: Each group consists of only one observation. Do you need to adjust the group aesthetic?
我现在想要 运行 模型的约束版本，基本上强制组的主效应等于 0。为此，我按照 'merge' 处理和时间因素与共同的基线水平相结合，然后是每个 post 治疗的新水平（第 12 页，以下 link）：
https://bendixcarstensen.com/Notes/rm.pdf
# CONSTRAINED LDA (mixed model) ----
# Create a 'merged' treat*time variable with 3 levels:
# - a common pre level for both treatments
# - a new post level for each treatment
dat <- transform(dat, treat_time = Epi::Relevel(interaction(treat, time), list(baseline = 1:2)))
# Inspect levels
with(dat, ftable(treat, time))
with(dat, ftable(treat_time, treat, time))
> m4 <- lme4::lmer(y ~ treat_time + (1|id), data = dat)
> summary(m4)
Linear mixed model fit by REML ['lmerMod']
Formula: y ~ treat_time + (1 | id)
Data: dat
REML criterion at convergence: 268.2
Scaled residuals:
Min 1Q Median 3Q Max
-1.7986 -0.5202 0.1392 0.6089 1.2810
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 25.74 5.074
Residual 46.02 6.784
Number of obs: 40, groups: id, 20
Fixed effects:
Estimate Std. Error t value
(Intercept) 56.016 1.894 29.571
treat_timeveh.post -2.027 2.902 -0.698
treat_timeil_17.post 11.270 2.676 4.212
Correlation of Fixed Effects:
(Intr) trt_t.
trt_tmvh.ps -0.419
trt_tml_17. -0.454 0.190
> emmeans(m4, revpairwise ~ treat_time)
$emmeans
treat_time emmean SE df lower.CL upper.CL
baseline 56.0 1.89 32.9 52.2 59.9
veh.post 54.0 2.79 37.0 48.3 59.6
il_17.post 67.3 2.53 36.9 62.2 72.4
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
veh.post - baseline -2.03 2.97 24.2 -0.684 0.7752
il_17.post - baseline 11.27 2.72 22.8 4.141 0.0011
il_17.post - veh.post 13.30 3.72 32.5 3.573 0.0031
Degrees-of-freedom method: kenward-roger
P value adjustment: tukey method for comparing a family of 3 estimates
> emmip(m4, ~ treat_time) +
+ geom_point() +
+ geom_line(size = 1) +
+ scale_y_continuous(limits = c(50, 70), breaks = seq(50, 70, by = 2)) +
+ theme_bw(base_size = 15)
geom_path: Each group consists of only one observation. Do you need to adjust the group aesthetic?
>
但我不知道如何（或者是否有办法）使用 emmip 生成类似的图，考虑到因子水平是组合的。我希望能够从本质上重现上面的第一个图，显示每个组从公共基线的变化。我想我可以在 ggplot 中手动执行此操作 - 只是想知道是否有某种方法可以直接使用 emmip 进行操作？
dat <- structure(list(id = c(1L, 1L, 2L, 2L, 3L, 3L, 4L, 4L, 5L, 5L,
6L, 6L, 7L, 7L, 8L, 8L, 9L, 9L, 10L, 10L, 11L, 11L, 12L, 12L,
13L, 13L, 14L, 14L, 15L, 15L, 16L, 16L, 17L, 17L, 18L, 18L, 19L,
19L, 20L, 20L), treat = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L), .Label = c("veh", "il_17"), class = "factor"), time = structure(c(1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L), .Label = c("pre", "post"), class = "factor"),
y = c(63.9871851, 62.7837203, 58.4825871, 59.4795929, 59.7735716,
42.4923719, 51.0754017, 55.5395164, 55.8808155, 42.8811745,
38.686015, 40.6969206, 46.8213372, 61.6, 62.6084263, 62.987013,
54.5901379, 53.0518311, 60.99353, 65.8177226, 66.2765671,
77.6594363, 60.1260258, 80.5687204, 49.444707, 63.7658228,
55.1138485, 77.6594363, 61.4070066, 80.5687204, 52.4932976,
49.444707, 55.1138485, 57.5908965, 61.4070066, 54.4305875,
52.4932976, 68.6955215, 53.5461576, 68.3395943), treat_time = structure(c(1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L,
1L, 2L, 1L, 3L, 1L, 3L, 1L, 3L, 1L, 3L, 1L, 3L, 1L, 3L, 1L,
3L, 1L, 3L, 1L, 3L, 1L, 3L, 1L, 3L), .Label = c("baseline",
"veh.post", "il_17.post"), class = "factor")), class = "data.frame", row.names = c(NA,
-40L))
我想如果你对每个治疗使用相同的基线值，你可以这样做：
foo = emmeans(m4, "treat_time")
bar = contrast(foo, list(b1 = c(1,0,0), p1 = c(0,1,0),
b2 = c(1,0,0), p2 = c(0,0,1)) )
levels(bar) = list(time = c("baseline", "post"),
treat = c("veh", "il_17") )
emmip(bar, treat ~ time)
更简单的方法
我发现你可以大大简化第二步：
foo = emmeans(m4, "treat_time")
bar = foo[c(1,2,1,3)]
levels(bar) = list(time = c("baseline", "post"),
treat = c("veh", "il_17") )
emmip(bar, treat ~ time)