我尝试创建一个生存预测图表
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
但这是一次不成功的尝试。感谢任何想法
最佳答案
通过 reshape 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)
这给出了以下情节:
关于r - 如何将误差范围的阴影添加到图表中,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/59824646/