r - 如何通过针对协变量等进行调整的治疗来绘制预测因子?
问题描述
我想用三次预测变量和许多协变量和交互调整来绘制拟合的治疗效果。不过,ggplot
我可以很容易地按处理对数据进行分组并添加一个geom_smooth()
来获得它,而无需调整。我将此答案应用于我的问题,但是当您在模型中进行大量调整时会很痛苦,而当您有长格式的数据时几乎不适用。所以我的问题是是否有更简单的方法来获得我想要的东西。
一些数据
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
因此,与习惯情节中的真实模型相比,向我展示的内容存在很大差异。
绘制无协变量
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)
将上述解决方案应用于我的问题会给我以下信息。
预言
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)
看起来好像是我想要的,不过我不太信任它。无论如何,是否有更简单,更可靠的方法来获得这样的情节?
解决方案
我知道这个问题很老,所以仅供参考!
我会选择这个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 系统。
推荐阅读
- java - 插件更新错误“某些插件需要插件 JNA 才能安装 netbeans”
- python - 指定 Windows 路径时找不到文件错误
- javascript - 拖动桌面文档时更改鼠标光标
- jupyter-lab - 如何在 JupyterLab 中拆分和合并单元格
- azure - 用于运行生成或管道的 Azure REST API
- c++ - typename 必须在模板中指定吗?
- c++ - C++ 中的 _Tidy_deallocate_guard 是做什么用的?
- python - 缺少 1 个必需的位置参数:django 模型中的“on_delete”
- python - groupby 应用函数中的缺失值
- jenkins - 我可以将我的 Jenkins 管道编写为 Yaml 而不是 groovy 吗?