Блог

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

  

Возможно, это соответствует вашим потребностям. Пожалуйста, обратите внимание, что я включил только 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()