在同一ggplot上绘制多个模型规格

debugcn 发表于 Dev

Beddotcom

我正在使用线性模型，我想使用条件之间均值差异的单点范围图来显示效果对不同规格的鲁棒性。据我所知（下文MWE）。

我有三个重要的虚拟治疗指标，以及五个协变量。

我现在想做的是用五个其他模型的图叠加该图中三种处理的估计，每个模型都包含不同的协变量，然后添加图例/形状/颜色以区分它们。我假设我可以group_by()和do()五个单独的模型一起使用，但是置信区间名称被替换了，而且我不确定如何让ggplot读取多个模型（尤其是在tidyverse中，这对我来说是陌生的）。

我无法弄清楚或找不到任何可以解决此类问题的现有线程。能做到吗？感谢您的提示！

具有示例数据的MWE：

treatment1 = rep(seq(0, 1, 1), 300)
treatment2 = sample(seq(from = 0, to = 1, by = 1), size = 300, replace = TRUE)
treatment3 = rep(seq(0, 1, 1), each=300)
response = rnorm(n = 300, mean = 3, sd = 1)
cov1 = rnorm(n = 300, mean = 0, sd = 1)*response
cov2 = rnorm(n = 300, mean = 0, sd = 1)/response
cov3 = rnorm(n = 300, mean = 0, sd = 1)-response
cov4 = rnorm(n = 300, mean = 0, sd = 1)+response
cov5 = rnorm(n = 300, mean = 0, sd = 1)*log(response)
df <- data.frame(treatment1,treatment2,treatment3,
           response,cov1,cov2,cov3,cov4,cov5)

model <- df %>% group_by(treatment1, treatment2,treatment3) %>% 
  do(data.frame(tidy(lm(response ~ treatment1*treatment2*treatment3, data = .),conf.int=T, conf.level = 0.95 )))

facet.labs <- c("T1=0", "T1=1")
names(facet.labs) <- c("0", "1")


model$treatment3 <-factor(model$treatment3, labels = c("T3=1","T3=0"))
model$treatment3 <-factor(model$treatment3, levels = c("T3=1","T3=0"))

ggplot(model, aes(x=estimate, y=treatment2, shape = treatment3)) +
  geom_pointrange(position = position_dodge(width = 1), aes(xmin=conf.low, xmax=conf.high), size=.75) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(),
        axis.line = element_line(colour = "black"),
        panel.border = element_rect(colour = "black", fill=NA, size=.8)) +
  scale_y_continuous(name ="", breaks = c(.1,1.22), labels=c("T2=0","T2=1")) +
geom_vline(xintercept=0, linetype="dotted") +
  labs(title="") + xlab("")  +
  labs(shape="")+ theme(axis.ticks = element_blank()) +
  theme(axis.text.y = element_text(angle = 90, vjust = 0.5, hjust=1)) +#,col="Treatment 1")+
  guides(fill = guide_legend(override.aes = list(linetype = 0,fill=NA)))+
  facet_wrap(~treatment1, labeller = labeller(treatment1 = facet.labs)) +
  guides(shape = guide_legend(reverse=T))

YBS

也许这符合您的需求。请注意，我只包括3个型号。您可以添加更多。

library(broom)

treatment1 = rep(seq(0, 1, 1), 300)
treatment2 = sample(seq(from = 0, to = 1, by = 1), size = 300, replace = TRUE)
treatment3 = rep(seq(0, 1, 1), each=300)
response = rnorm(n = 300, mean = 3, sd = 1)
cov1 = rnorm(n = 300, mean = 0, sd = 1)*response
cov2 = rnorm(n = 300, mean = 0, sd = 1)/response
cov3 = rnorm(n = 300, mean = 0, sd = 1)-response
cov4 = rnorm(n = 300, mean = 0, sd = 1)+response
cov5 = rnorm(n = 300, mean = 0, sd = 1)*log(response)
df1 <- data.frame(treatment1,treatment2,treatment3,
                 response,cov1,cov2,cov3,cov4,cov5)

facet.labs <- c("T1=0", "T1=1")
names(facet.labs) <- c("0", "1")

model1 <- df1 %>% group_by(treatment1, treatment2,treatment3) %>% 
  do(data.frame(tidy(lm(response ~ treatment1*treatment2*treatment3, data = .),conf.int=T, conf.level = 0.95 )))

#model1$treatment3 <-factor(model1$treatment3, labels= c("T3=1","T3=0"), levels = c("T3=1","T3=0"))

model11 <- data.frame(model1,model=1)

treatment1 = rep(seq(0, 1, 1), 300)
treatment2 = sample(seq(from = 0, to = 1, by = 1), size = 300, replace = TRUE)
treatment3 = rep(seq(0, 1, 1), each=300)
response = rnorm(n = 300, mean = 4, sd = 1)
cov1 = rnorm(n = 300, mean = 0, sd = 1)*response*2
cov2 = rnorm(n = 300, mean = 0, sd = 1)/response
cov3 = rnorm(n = 300, mean = 0, sd = 1)-response
cov4 = rnorm(n = 300, mean = 0, sd = 1)+response
cov5 = rnorm(n = 300, mean = 0, sd = 1)*log2(response)
df2 <- data.frame(treatment1,treatment2,treatment3,
                  response,cov1,cov2,cov3,cov4,cov5)

model2 <- df2 %>% group_by(treatment1, treatment2, treatment3) %>% 
  do(data.frame(tidy(lm(response ~ treatment1*treatment2*treatment3, data = .),conf.int=T, conf.level = 0.95 )))

#model2$treatment3 <-factor(model2$treatment3, labels = c("T3=1","T3=0"), levels = c("T3=1","T3=0"))

model22 <- data.frame(model2,model=2)

cov1 = rnorm(n = 300, mean = 0, sd = 1)*response*0.5
cov5 = rnorm(n = 300, mean = 0, sd = 1)*log10(response)
df3 <- data.frame(treatment1,treatment2,treatment3,
                  response,cov1,cov2,cov3,cov4,cov5)

model3 <- df3 %>% group_by(treatment1, treatment2, treatment3) %>% 
  do(data.frame(tidy(lm(response ~ treatment1*treatment2*treatment3, data = .),conf.int=T, conf.level = 0.95 )))

#model3$treatment3 <-factor(model2$treatment3, labels = c("T3=1","T3=0"), levels = c("T3=1","T3=0"))

model33 <- data.frame(model3,model=3)

model <- rbind(model11,model22,model33)

myshapes <- c(15, 17)
mycolors <- c("blue","orange")
mygroup <- c("T3=1","T3=0")

modelb <- transform(model,trt2_model = paste0("model ",model, " - trt2 ", treatment2))

ggplot(modelb, aes(x=estimate, y=trt2_model, xmin=conf.low, xmax=conf.high, 
                  shape = factor(treatment3), color=factor(treatment3) )) + 
  geom_pointrange(position = position_dodge(width = 1), aes(xmin=conf.low, xmax=conf.high), size=.75) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(),
        legend.key = element_rect(fill = "white", colour = "white"),
        axis.ticks = element_blank(),
        axis.line = element_line(colour = "black"),
        #axis.text.y = element_text(angle = 90, vjust = 1, hjust=0.5),
        panel.border = element_rect(colour = "blue", fill=NA, size=.8)) +
  #scale_y_continuous(name ="", breaks = c(.1,1.22), labels=c("T2=0","T2=1")) +
  geom_vline(xintercept=0, linetype="dotted", lwd=1, color="red") +
  labs(shape="", x="", y="", title="")+ 
  scale_shape_manual(name = " ",
                     labels = mygroup,
                     values = myshapes) +  ## choice of shapes
  scale_color_manual(name = " ",
                     labels = mygroup,
                     values = mycolors ) +  ## colors of your choice
  
  guides(color='none', fill = guide_legend(override.aes = list(linetype = 0,fill=NA)))+
  facet_wrap(~treatment1, labeller = labeller(treatment1 = facet.labs)) +
  guides(shape = guide_legend(override.aes=list(col=mycolors, lty=0, pt.cex=1.5, reverse=T)) ) +
  theme_bw()