如何通过针对协变量等调整的治疗来绘制预测因子？

Question

我想绘制具有立方预测变量和调整的许多协变量和交互作用的拟合的治疗效果。使用 ggplot 我可以轻松地按治疗对数据进行分组并添加 geom_smooth() 来获得它，但无需调整。我对我的问题应用了 this answer，但是当你在你的模型中进行大量调整并且当你有长格式数据时几乎不适用时这会很痛苦。所以我的问题是是否有更简单的方法来获得我想要的东西。

一些数据

set.seed(42)
n <- 1e4
D <- rbinom(n, 1, .5)  # treatment indicator
X <- .5 + rnorm(n)  # bunch of covariates and other adjustemnts
P <- 5.54 + 0.35*D -.24*X + rnorm(n)  # predictor
Y <- 1.49 - 1.35*P + .5*P^2 - 0.04*P^3 - 0.83*D + 0.43*X + rnorm(n, 0, 6)
df1 <- data.frame(D, X, P, Y)

模型，指定完整和不完整

true <- lm(Y ~ P + I(P^2) + I(P^3) + D + X , df1)  # true model
bias <- lm(Y ~ P + I(P^2) + I(P^3) + D, df1)

> round(rbind(true=coef(true), bias=c(coef(bias), NA)),
+       3)
     (Intercept)     P I(P^2) I(P^3)      D    X
true      -4.023 1.803 -0.088 -0.005 -0.728 0.42
bias      -3.426 1.753 -0.091 -0.005 -0.702   NA

所以 ggplot 将向我展示的内容与常规情节中的真实模型相比有很大差异。

绘制 w/o 协变量

library(ggplot2)
p1 <- ggplot(df1, aes(P, Y, color=as.factor(D), group=D)) +
  geom_smooth(se=FALSE) +
  theme_bw()
p2 <- ggplot(df1, aes(P, Y, color=as.factor(D), group=D)) +
  stat_smooth(method="lm", formula=y ~ poly(x, 3, raw=TRUE), se=FALSE) +
  theme_bw()
egg::ggarrange(p1, p2)

将 mentioned solution 应用于我的问题得到以下结果。

预测

n.data <- data.frame(D=rep(range(D), each=n/2),
                   P=rep(seq(range(df1$P)[1], range(df1$P)[2],
                              length.out=n/2), times=2),
                   X=rep(seq(range(df1$X)[1], range(df1$X)[2],  # assume this dozens of!
                             length.out=n/2), times=2))
df1.2 <- data.frame(n.data, pred=predict(true, n.data))

带协变量的绘图

p1a <- ggplot(df1, aes(x=P, y=Y, color=as.factor(D)))  +
  geom_smooth(data=df1.2, aes(x=P, y=pred, color=as.factor(D))) +
  theme_bw()
p2a <- ggplot(df1, aes(x=P, y=Y, color=D)) +
  stat_smooth(method="lm", formula=y ~ poly(x, 3, raw=TRUE),
              data=df1.2, aes(x=P, y=pred, color=as.factor(D))) +
  theme_bw()
egg::ggarrange(p1a, p2a)

看起来是我想要的，不过我不是很信任。无论如何，是否有更简单、更可靠的方法来获得这样的图？

Answer 1

我知道这个问题很老了，所以仅供参考！

我会选择 effects 套餐：

set.seed(42)
n <- 1e4
D <- rbinom(n, 1, .5)  # treatment indicator
X <- .5 + rnorm(n)  # bunch of covariates and other adjustemnts
P <- 5.54 + 0.35*D -.24*X + rnorm(n)  # predictor
Y <- 1.49 - 1.35*P + .5*P^2 - 0.04*P^3 - 0.83*D + 0.43*X + rnorm(n, 0, 6)
df1 <- data.frame(D = factor(D, labels = c("Control", "Treatment")), X, P, Y)
true <- lm(Y ~ poly(P, 3, raw = TRUE):D + X , df1)  # true model
library(effects)
plot(predictorEffect("P", true), lines=list(multiline=TRUE))

如果你想要 ggplot，可以使用 ggeffects 包，它的功能基本相同，但使用的是 ggplot2 系统。

如何通过针对协变量等调整的治疗来绘制预测因子？

How to plot predictor by treatment adjusted for covariates etc.?

r

prediction

ggplot2