自行编写的survival curve的function，旨在画出生存曲线和曲线下的no at risk表。
结果运行后报错，求帮助。
报错trackback如下：
Error in parse(text = x, keep.source = FALSE) :
<text>:2:0: unexpected end of input
1: srv~
^
7.
parse(text = x, keep.source = FALSE)
6.
eval(parse(text = x, keep.source = FALSE)[[1L]])
5.
formula(eval(parse(text = x, keep.source = FALSE)[[1L]]))
4.
formula.character(object, env = baseenv())
3.
formula(object, env = baseenv())
2.
as.formula(paste("srv~", marker))
1.
ggsrv(time = "RFS", event = "resultR", data = data)
原代码如下：
ggsrv<-function(
time="os",
event="death",
event.marker=1,
marker="",
tabletitle="",
xlab="Months since surgery",
ylab="Overall Survival rate",
ystratalabs=c("High", "Low"),
pval=TRUE,
n.risk=FALSE,
timeby=2,
xmax=16,
fun="survival",
data = dat,
conf.int=FALSE,
returns = F,
...)
{
#Load the related packages
require(ggplot2)
require(survival)
require(gridExtra)
require(grid)
require(plyr)
require(cowplot)
#create survival fit and data
srv <- Surv(data[,time], data[,event]==event.marker)
form <- as.formula(paste("srv~",marker))
sfit=survfit(form,data = data)
sdiff=survdiff(form,data = data)
if(is.null(ystratalabs)) {
ystratalabs <- as.character(levels(summary(sfit)$strata))
}
m <- max(nchar(ystratalabs))
# if(xmax<=max(sdata$time)){
# xlims=c(0, round(xmax/timeby, digits=0)*timeby)
# }else{
# xlims=c(0, round((max(sdata$time))/timeby, digits=0)*timeby)
# }
times <- seq(0, xmax, by = timeby)
sdata <- data.frame(time = sfit$time,
n.risk = sfit$n.risk,
n.event = sfit$n.event,
surv = sfit$surv,
strata = summary(sfit, censored = T)$strata,
upper = sfit$upper,
lower = sfit$lower,
n.censor = sfit$n.censor)
levels(sdata$strata) <- ystratalabs
zeros <- data.frame(time = 0,
surv = 1,
strata = factor(ystratalabs, levels=levels( sdata$strata)),
upper = 1,
lower = 1,
n.censor=0)
sdata <- rbind.fill(zeros, sdata)
d <- length(levels(sdata$strata))
dat.cens <- subset(sdata, n.censor != 0)
if (fun=="event"){
p <- ggplot(sdata, aes(x = time, y = 1-surv,group = strata, colour = strata)) +
geom_step(size = 1) +
scale_x_continuous(xlab, breaks = times, limits = c(0, max(times)), expand = c(0, 0)) +
scale_y_continuous(ylab, breaks=seq(0,1,by=0.2), limits=c(0,1),
labels=c(as.character(seq(0,0.8,by=0.2)),"1.0"),expand = c(0, 0)) +
theme_classic() +
theme(axis.title.x = element_text(vjust = 0.5)) +
theme(panel.grid.minor = element_blank()) +
theme(legend.position = c(ifelse(m < 10, .28, .15), ifelse(d < 4, .75, .95))) +
theme(legend.key = element_rect(colour = NA)) +
ggtitle(tabletitle)+
geom_point(data = subset(sdata, n.censor != 0), aes(x = time, y = 1-surv), shape = 3)
}else{
p <- ggplot(sdata, aes(x = time, y = surv,group = strata, colour = strata)) +
geom_step(size = 1) +
scale_x_continuous(xlab, breaks = times, limits = c(0, max(times)), expand = c(0, 0)) +
scale_y_continuous(ylab, breaks=seq(0,1,by=0.2), limits=c(0,1),
labels=c(as.character(seq(0,0.8,by=0.2)),"1.0"),expand = c(0, 0)) +
theme_classic() +
theme(axis.title.x = element_text(vjust = 0.5)) +
theme(panel.grid.minor = element_blank()) +
theme(legend.position = c(ifelse(m < 10, .28, .15), ifelse(d < 4, .35, .45))) +
theme(legend.key = element_rect(colour = NA)) +
ggtitle(tabletitle)+
geom_point(data = subset(sdata, n.censor != 0), aes(x = time, y = surv), shape = 3)
}
#Add P value to plot
if(pval) {
sdiff <- survdiff(eval(sfit$call$formula), data = eval(sfit$call$data))
pval <- pchisq(sdiff$chisq, length(sdiff$n)-1, lower.tail = FALSE)
pvaltxt <- ifelse(pval < 0.0001, "p < 0.0001", paste("p =", signif(pval, 3)))
p <- p + annotate("text", x = 0.75 * max(times), y = 0.15, label = pvaltxt)
}
if(n.risk) {
## Create table graphic to include at-risk numbers
risk.data <- data.frame(strata = summary(sfit, times = times, extend = TRUE)$strata,
time = as.numeric(as.factor(summary(sfit, times = times, extend = TRUE)$time)),
n.risk = summary(sfit, times = times, extend = TRUE)$n.risk)
data.table <- ggplot(risk.data, aes(x = time, y = strata,
label = format(n.risk, nsmall = 0))) +
#, color = strata)) +
geom_text(size=4) +
scale_y_discrete(breaks = as.character(levels(risk.data$strata)),
labels = ystratalabs) +
scale_x_discrete(expand = c(0, 0))+
theme(
panel.grid.major = element_blank(),
legend.position = "none",
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
legend.position="none",
axis.line = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_text(size=10, face="bold", color = 'black'),
axis.ticks=element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
plot.title = element_blank()
)
both = rbind(ggplotGrob(p), ggplotGrob(data.table), size="last")
panels <- both$layout$t[grep("panel", both$layout$name)]
both$heights[panels] <- list(unit(1,"null"), unit(2, "lines"))
both <- gtable_add_rows(both, heights = unit(1,"line"), 8)
both <- gtable_add_grob(both,
textGrob("Number at risk", hjust=0, x=0),
t=9, l=2, r=4)
plot_grid(both)
} else {
return(p)
}
}
