r - 如何在图表中添加误差范围的阴影

Question

我尝试创建一个生存预测'图表

library("survival")
# fit regression
res.cox <- coxph(Surv(time, status) ~ age + sex + wt.loss, data =  lung)
res.cox

拟合新数据

sex_df <- with(lung,
               data.frame(sex = c(1, 2), 
                          age = rep(mean(age, na.rm = TRUE), 2),
                          wt.loss = rep(mean(wt.loss, na.rm = TRUE), 2)  ))

图表

library("ggplot2")
fit <- survfit(res.cox, newdata = sex_df)
library(reshape2)
dat = data.frame(surv = fit$surv,lower= fit$lower, upper = fit$upper,time= fit$time)
head(dat)
head(melt(dat, id="time"))
data = melt(dat, id="time")

obj = strsplit(as.character(data$variable), "[.]") # делим текст на объекты по запятой

data$line = sapply(obj, '[', 1)
data$number = sapply(obj, '[', 2)

ggplot(data, aes(x=time, y=value, group=variable)) +
  geom_line(aes(linetype=line, color=as.factor(number), size=line)) +
  # geom_point(aes(color=number)) +
  theme(legend.position="top", axis.text = element_text(size = 20), 
        axis.title = element_text(size = 20), 
        legend.text=element_text(size=40),
        legend.key.size = unit(3,"line"))+
  scale_linetype_manual(values=c( 2,1,2))+ # "dotted", "twodash","dotted"
  scale_color_manual(values=c("#E7B800", "#2E9FDF", 'red'))+
  scale_size_manual(values=c(2, 3.5, 2)) +
  scale_x_continuous(limits=c(0, 840),
                     breaks=seq(0, 840, 120)) + ylab("Surv prob") + 
  guides(linetype = FALSE, size = FALSE, color = guide_legend(override.aes = list(size=5))) + labs(color='') + 
  geom_ribbon(aes(ymin = rep(data$value[data$line == 'lower' & 
                                          data$number == "1"],6), 
       ymax = rep(data$value[data$line == 'upper' & data$number == "1"],6)), 
       fill = "#E7B800",alpha=0.1) +
  geom_ribbon(aes(ymin = rep(data$value[data$line == 'lower' & data$number == "2"],6), 
                  ymax = rep(data$value[data$line == 'upper' & data$number == "2"],6)), 
              fill = "#2E9FDF",alpha=0.1) 

问题 图表没问题，但我必须用手加上这个

geom_ribbon(aes(ymin = rep(data$value[data$line == 'lower' & data$number == "2"],6), 
                      ymax = rep(data$value[data$line == 'upper' & data$number == "2"],6)), 
                  fill = "#2E9FDF",alpha=0.1) 

如果新数据中有三个但不是两个元素，您将不得不重写代码。是否可以重写代码，使其不依赖于新数据的元素数量？我尝试使用循环

temp = list()
uniq <- unique(unlist(data$number))
for (i in 1:length(levels(as.factor(data$number)))) {
  n = geom_ribbon(aes(ymin = rep(data$value[data$line == 'lower' & data$number == uniq[i]],6), 
                        ymax = rep(data$value[data$line == 'upper' & data$number == uniq[i]],6)), 
                  fill = "#2E9FDF", alpha=0.1) # 
  temp = append(n, temp)
  }
temp

但这是一次不成功的尝试。感谢您的任何想法

标签： rggplot2rstudioreshapesurvival

Answer 1

通过重塑 data.frame 使surv,lower和upper是单独的向量，您可以geom_ribbon按元素而不是行的“含义”对 进行分组。

下面是使用tidyr包的代码；第一部分只是您生成数据的代码。

library(survival)
library(reshape2)
library(ggplot2)

# fit regression
res.cox <- coxph(Surv(time, status) ~ age + sex + wt.loss, data =  lung)
res.cox

sex_df <- with(lung,
               data.frame(sex = c(1, 2), 
                          age = rep(mean(age, na.rm = TRUE), 2),
                          wt.loss = rep(mean(wt.loss, na.rm = TRUE), 2)  ))


fit <- survfit(res.cox, newdata = sex_df)

dat = data.frame(surv = fit$surv,lower= fit$lower, upper = fit$upper,time= fit$time)

head(dat)
head(melt(dat, id="time"))
data = melt(dat, id="time")

# Reformats the data into format with the survival curve and the confidence intervals in their own columns
library(tidyr)

data_wide <- data %>%
  separate(col = variable, into = c("type", "sex"), sep = "\\.") %>%
  spread(key = type, value = value)


ggplot(data = data_wide) +
  geom_line(aes(x = time, y = surv, group = sex, colour = sex),
            size = 3.5,
            linetype = 1) +
  geom_line(aes(x = time, y = lower, group = sex, colour = sex),
            size = 2,
            linetype = 2) +
  geom_line(aes(x = time, y = upper, group = sex, colour = sex),
            size = 2,
            linetype = 2) +
  # Geom_ribbom now grouped by sex
  geom_ribbon(aes(x = time, ymin = lower, ymax = upper, group = sex, fill = sex),
              alpha = 0.1) +
  scale_colour_manual(values = c("#E7B800", "#2E9FDF")) +
  scale_fill_manual(values = c("#E7B800", "#2E9FDF")) +
  scale_x_continuous(limits = c(0, 840),
                     breaks = seq(0, 840, 120)) +
  theme(legend.position = "top",
        axis.text = element_text(size = 20),
        axis.title = element_text(size = 20),
        legend.text = element_text(size = 40),
        legend.key.size = unit(3, "line")) +
  ylab("Surv prob")

这是绘图输出：

我们添加了另一个元素来测试这是否有效，您将不得不向scale_colour_manual和添加更多颜色scale_fill_manual。

library(dplyr)
data_wide2 <- filter(data_wide, sex == "1") %>%
  mutate(sex = "3",
         surv = surv - 0.2,
         upper = upper - 0.2,
         lower = lower - 0.2) %>%
  rbind(data_wide)

这给出了以下情节：