# TODO: Add comment # # Author: E.Korsching 4.8.2011 ############################################################################### multiSurv <- function(x, t.surv, d.a, eff, xlab="time [month]", ylab="survival [%]", t.text=NULL, l.title=NULL, l.text=NULL, bty="n", cox=F, pv=NULL, xl=10, yl=0.4, adj=c(0,0), ci="none", subfolder.filename="/results/a1", add=0, lty=1, col=c("blue","turquoise1","limegreen","yellow4","brown1","orangered4"), text.col="black" ){ # try to find master effectors by performing all single dependency survival plots # and one cox combined effect model # get: data.frame with all information: x, survival period : t.surv, dead or alive : d.a (at the end of the observartion period) 0,1 F,T 1,2 (alive,dead) # vector of effector col(s): eff (will be decomposed into levels) # na.action: is: options()$na.action: default is normally: na.omit() function , can be explicitly applied by na.omit(x) # col identifier: number or name (name : bug in legend!) # t.text: title # l.title: legend titel # l.text: legend items - needs to be the same number of entries as in levels # bty: "o": box, "n": no box see par() # pv: p value at yl+pv position or NULL # text.col: color of the legend text # ci: confidence intervals: "none", "plain", "log", "log-log" # subfolder.filename="/results/a1" : w/o subfolder only "/a1" or NULL # cox: additionally cox ph ? : T:yes or F:no # add: 0,1,2 : 0: full plot, 1: less full plot, 2: add line plots to existing (survival) plot # xl: horizontal legend position # yl: vertical legend position: 0.4 (standard start), and e.g. 0.32, 0.24, 0.16, ... # lty: one type (1: solid) -or- a sufficient number of types (integers) # col: vector of colors or intergers (for colors) - sufficient number #eval(parse(text=paste(c(1,2,3), collapse="+"))) #ini nrx <- nrow(x) ncy <- ncol(x) names.x <- names(x) if(!is.null(subfolder.filename)){ res.split <- unlist(strsplit(x=subfolder.filename, split="/", fixed=T)) len.split <- length(res.split) } eff.len <- length(eff) #count number of survival plots cex <- par("cex") # require(survival) require(plotrix) #addtable2plot - existing plot will be used, also: require(gplots):textplot - new plot will be created #work #stat1 stat.n <- nrow(x) # number of values # Kaplan-Meier if(!is.null(subfolder.filename)){ pdf(file=paste(getwd(),subfolder.filename,".","Kaplan-Meier",".",format(Sys.time(), "%Y%m%d%H%M"),".pdf", sep=""), width=11, height=7, onefile=T, title=res.split[len.split], pointsize=12) } for(i in 1:eff.len){ #trigger p values off in case pvT <- 1 # on #stat2 eff.i <- eff[i] stat.na <- sum(is.na(x[,eff.i])) # how many NA # if(is.null(stat.na)){ stat.na <- 0 } stat.nlevel <- nlevels(factor(x[,eff.i])) # number of levels stat.npl <- as.data.frame(table(factor(x[,eff.i])))[,2] # count per level if(stat.nlevel==1){ fit1 <- survfit( Surv(x[,t.surv], x[,d.a]) ~ 1, type="kaplan-meier", conf.type=ci) pvT <- 0 # off : no log-rank test etc. - only one curve }else{ fit1 <- survfit( Surv(x[,t.surv], x[,d.a]) ~ x[,eff.i], type="kaplan-meier", conf.type=ci) #default normally: "na.omit" test.diff1 <- survdiff( Surv(x[,t.surv], x[,d.a]) ~ x[,eff.i], rho=0) #log-rank test on difference between curves cat("\n in test : n=",test.diff1$n) p.val1 <- 1 - pchisq(test.diff1$chisq, length(test.diff1$n) - 1) #p value(s) test.diff2 <- survdiff( Surv(x[,t.surv], x[,d.a]) ~ x[,eff.i], rho=1) #Peto & Peto modification of the Gehan-Wilcoxon test on difference between curves p.val2 <- 1 - pchisq(test.diff2$chisq, length(test.diff2$n) - 1) #p value(s) } # png(filename=paste(getwd(),subfolder.filename,".",format(Sys.time(), "%Y%m%d%H%M."),i,".png", sep=""), # width=800, height=800, units="px", pointsize=12, bg="white") cat("\n col : ",names.x[eff.i]," levels : ",stat.nlevel) if(add==0){ # full plot plot(fit1, xlab=xlab, ylab=ylab, lty=lty, col=col, mark.time=T) if(is.null(l.text)){ ll.text <- as.character(levels(factor(x[,eff.i]))) }else{ if(pvT==0){ ll.text <- l.text[1] }else{ ll.text <- l.text } } legend(x=xl,y=yl, legend=ll.text, lty=lty, col=col, title=if(is.null(l.title)){names.x[eff.i]}else{l.title}, bty=bty, text.col=text.col ) #plot coordinates if(!is.null(t.text)){ title(main=paste(t.text,"n=",stat.n,", NA:",stat.na,", level(s):",stat.nlevel,", count per level:",paste(stat.npl,collapse=","),sep=" ")) } } if(add==2){ # only curves lines(fit1, type="s", lty=lty, col=col, mark.time=T) # and legend if(is.null(l.text)){ ll.text <- as.character(levels(factor(x[,eff.i]))) }else{ if(pvT==0){ ll.text <- l.text[1] }else{ ll.text <- l.text } } legend(x=xl,y=yl, legend=ll.text, lty=lty, col=col, title=l.title, bty=bty, text.col=text.col ) } if(!is.null(pv) & pvT==1){ # p value text(x=xl,y=yl+pv, labels=paste("log-rank test - p :", round(p.val1,digits=4), ", ", "Peto&Peto modification of Gehan-Wilcoxon test - p :", round(p.val2,digits=4)), adj=adj, vfont=NULL, font=NULL, cex=cex-0.1, col=text.col) } # dev.off() } if(!is.null(subfolder.filename)){ dev.off() } cat("\n") # Cox proportional hazards if(cox){ if(!is.null(subfolder.filename)){ pdf(file=paste(getwd(),subfolder.filename,".","Cox.ph",".",format(Sys.time(), "%Y%m%d%H%M"),".pdf", sep=""), width=11, height=7, onefile=T, title=res.split[len.split], pointsize=12) } for(i in 1:eff.len){ par(mfrow=c(1,2)) eff.i <- eff[i] # cat(paste("Surv( ",names.x[t.surv], " , ", names.x[d.a], " ) ~ ",names.x[eff.i], sep=""), "\n") # try to keep the symbolic names tmp <- as.formula( paste("Surv( ",names.x[t.surv], " , ", names.x[d.a], " ) ~ ",names.x[eff.i], sep="") ) # do cox fit2 <- coxph( eval(tmp), data=eval(x), method="efron", model=F, x=F) #default normally: "na.omit"; other cox often use method='breslow' #print summary in left graph area plot(x=1,y=1, xlim=c(0,10), ylim=c(0,10), type="n", xlab="", ylab="", axes=F) tmp2 <- summary(fit2) tmp3 <- matrix(c(round(tmp2$logtest[1],3), round(tmp2$logtest[2],3), round(tmp2$logtest[3],3), round(tmp2$waldtest[1],3), round(tmp2$waldtest[2],3), round(tmp2$waldtest[3],3), round(tmp2$sctest[1],3), round(tmp2$sctest[2],3), round(tmp2$sctest[3],3)), 3,3,byrow=T) dimnames(tmp3)[[1]] <- c("Likelihood ratio test","Wald test","Score (logrank) test") dimnames(tmp3)[[2]] <- c("test","df","p") text(x=0.1, y=9.5, labels="Cox proportional hazards model (on the right side per level)", adj=c(0,0), cex=0.8) addtable2plot(x=0.1,y=8.5, round(tmp2$coefficients,3), bty="n", cex=0.8, display.colnames=T, display.rownames=T, hlines=T, vlines=F, xjust=0, yjust=0.5, title="Coefficients") addtable2plot(x=0.1,y=1.5, tmp3, bty="n", cex=0.8, display.colnames=T, display.rownames=T, hlines=T, vlines=F, xjust=0, yjust=0.5, title="Test values") # fit a Cox proportional hazards model and plot the predicted survival for each effector level # plot survival graph in the right graph area new.data <- as.numeric(levels( as.factor(x[,eff.i]) )) # levels from the effector variable cat("\n new.data",new.data) new.data.len <- length(new.data) if(new.data.len==1){ assign(x=names.x[eff.i], value=new.data[1]) plot(survfit(fit2, newdata=data.frame(get(x=names.x[eff.i], pos=-1, inherits=F)), conf.type="log"), xlab=xlab, ylab=ylab, lty=lty, col=col[1]) }else{ assign(x=names.x[eff.i], value=new.data[1]) plot(survfit(fit2, newdata=data.frame(get(x=names.x[eff.i], pos=-1, inherits=F)), conf.type="log"), xlab=xlab, ylab=ylab, lty=lty, col=col[1]) for(j in 2:new.data.len){ assign(x=names.x[eff.i], value=new.data[j]) lines(survfit(fit2, newdata=data.frame(get(x=names.x[eff.i], pos=-1, inherits=F)), conf.type="log"), type="s", lty=lty, col=col[j]) } } } par(mfrow=c(1,1)) # restore if(!is.null(subfolder.filename)){ dev.off() } } return() } #paste(strwrap(x=paste(rep("we ta ",40),collapse=""), width=20), collapse='\n') #capture.output #leitet output in var um #aa <- summary(md) #gibt ein obj zurück, das mit print ausgegeben wird #aa$logtest #oder einzel #aa$waldtest #aa$sctest #fit a Cox proportional hazards model and plot the #predicted survival for a 60 year old etc. #afit <- coxph(Surv(futime, fustat) ~ age, data = ovarian) #plot(survfit(afit, newdata=data.frame(age=60),conf.type="none"), xlab="Months", ylab="Survival") #lines(survfit(afit, newdata=data.frame(age=50),conf.type="none")) #lines(survfit(afit, newdata=data.frame(age=40),conf.type="none")) #lines(survfit(afit, newdata=data.frame(age=30),conf.type="none"))