我试图创建一个“生存预测”图表。
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
但这是一次失败的尝试。感谢任何想法。