如何在R中复刻Stata的marginsplot并按smoke分组展示age预测值？

受Chuck Huber的视频启发，我们可以通过调整margins包中cplot()函数的参数，实现类似Stata marginsplot的效果——按smoke的不同水平，展示age递增时的预测结果。

完整代码（含修改后的cplot()语句） #

library(margins)
set.seed(42)
n <- 1000
patient <- data.frame(id=1:n,
                      treat = factor(sample(c('Treat','Control'), n, rep=TRUE, prob=c(.5, .5))),
                      age=sample(18:80, n, replace=TRUE),
                      sex = factor(sample(c('Male','Female'), n, rep=TRUE, prob=c(.6, .4))),
                      smoke=factor(sample(c("Never", 'Former', 'Current'), n, rep=TRUE, prob=c(.25, .6, .15))),
                      outcome=runif(n, min=16, max=45))

model <- lm(outcome ~ treat*age + smoke, data = patient)

# 修改后的cplot语句
cplot(model, 
      x = "age", 
      by = "smoke", 
      overlay = TRUE,
      xlim = range(patient$age),  # 确保覆盖数据中所有年龄范围
      lwd = 2,  # 加粗线条提升可读性
      main = "Outcome预测值随age的变化（按smoke分组）",
      xlab = "年龄",
      ylab = "Outcome预测值")

关键修改说明 #

• xlim = range(patient$age)：强制x轴覆盖数据中age的全部取值范围，避免曲线截断
• lwd = 2：加粗分组线条，让不同smoke水平的曲线更易区分
• 添加main、xlab、ylab：自定义图表标题和坐标轴标签，让结果更直观

运行上述代码后，会生成一张叠加图：三条曲线分别对应smoke的三个水平，横轴为age，纵轴为模型预测的outcome值，完美复刻Stata marginsplot按分组展示连续变量预测趋势的效果。

内容的提问来源于stack exchange，提问作者Sandro