在 R 中创建“marginsplot”

我认为您正在寻找的东西可以用 sjPlot：:p lot_model（） 来完成：

library(ggplot2)
library(sjPlot)

plot_model(
    model,
    type = "pred",
    terms = c("age", "smoke"),
    ci.lvl = NA
) +
    theme_bw()

您可以使用软件包完成所有这些操作（免责声明：我是维护者）。在网站上，您会发现超过 25 个小插曲，包括一个完整的情节小插曲：marginaleffects

• 主网站： https://vincentarelbundock.github.io/marginaleffects/
• 入门：https://vincentarelbundock.github.io/marginaleffects/articles/marginaleffects.html
• 剧情简介： https://vincentarelbundock.github.io/marginaleffects/articles/plot.html
• 与 Stata 的比较：https://vincentarelbundock.github.io/marginaleffects/articles/alternative_software.html#stata

请注意，我添加了一个交互，使情节看起来更有趣：

library(marginaleffects)
library(ggplot2)
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)

plot_predictions(model, condition = c("age", "smoke")) + theme_minimal()

或者您可以绘制斜率：

plot_slopes(model, variables = "treat", condition = c("age", "smoke")) + theme_minimal()

通过复制 Stata 的 ，您想要的是沿着所需的 x 和 y 变量可视化模型中所有可能的离子，例如在本例中为年龄和烟雾。marginsplotmeanpredict

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

因此，首先，我们使用expand.grid

.newdata <- expand.grid(
  treat=unique(patient$treat),
  age=with(patient, min(age):max(age)),
  sex=unique(patient$sex),
  smoke=unique(patient$smoke)
  )

为了喂食，这结果我们.predictcbind

.newdata <- cbind(.newdata, predict(lm1, newdata=.newdata, interval='conf'))

接下来，我们计算烟雾和年龄变量的 ed 值的平均值，以及置信区间的各自和边界。aggregatefitlwrupr

agg <- aggregate(cbind(fit, lwr, upr) ~ smoke + age, .newdata, mean)

至此，我们已经完成了预处理，并准备好了。plot

par(mar=c(4, 4, 3, 2) + .1)
plot.new();plot.window(range(agg$age) + c(0, 2), range(agg[3:5]) + c(0, 2))
by(agg, agg$smoke, \(x) with(x, lines(age + as.integer(smoke) - 2, fit, col=smoke)))
dec <- agg$age %% 10 == 0
by(agg[dec, ], agg[dec, ]$smoke, \(x) 
   with(x, points(age + as.integer(smoke) - 2, fit, col=smoke, pch=20)))
by(agg[dec, ], agg[dec, ]$smoke, \(x) 
   with(x, arrows(age + as.integer(smoke) - 2, lwr, age + as.integer(smoke) - 2, upr, 
                  col=smoke, code=3, angle=90, length=.05)))
axis(1, axTicks(1)); axis(2, axTicks(2))
mtext('age', 1, 2.5); mtext('pred. outcome', 2, 2.5)
legend('topleft', pch=20, col=1:3, legend=unique(agg$smoke),
       title='smoke', horiz=TRUE, cex=.9)
box()

数据：

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