Я хочу нанести на график эффект лечения подгонки с кубическими предикторами и множеством ковариат и взаимодействий, откорректированных для. С помощью 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)
Похоже, это то, чего я хочу, хотя я не очень доверяю. В любом случае, может быть есть более простой и надежный способ получения таких графиков?