# Xpose 4
# An R-based population pharmacokinetic/
# pharmacodynamic model building aid for NONMEM.
# Copyright (C) 1998-2004 E. Niclas Jonsson and Mats Karlsson.
# Copyright (C) 2005-2008 Andrew C. Hooker, Justin J. Wilkins,
# Mats O. Karlsson and E. Niclas Jonsson.
# Copyright (C) 2009-2010 Andrew C. Hooker, Mats O. Karlsson and
# E. Niclas Jonsson.
# This file is a part of Xpose 4.
# Xpose 4 is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public License
# as published by the Free Software Foundation, either version 3
# of the License, or (at your option) any later version.
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
# You should have received a copy of the GNU Lesser General Public License
# along with this program. A copy can be cound in the R installation
# directory under \share\licenses. If not, see http://www.gnu.org/licenses/.
#' Kaplan-Meier plots of (repeated) time-to-event data
#'
#' Kaplan-Meier plots of (repeated) time-to-event data. Includes VPCs.
#'
#'
#' @param x The independent variable.
#' @param y The dependent variable. event (>0) or no event (0).
#' @param id The ID variable in the dataset.
#' @param data A dataset can be used instead of the data in an Xpose object.
#' Must have the same form as an xpose data object \code{xpdb@Data}.
#' @param evid The EVID data item. If not present then all rows are considered
#' events (can be censored or an event). Otherwise, EVID!=0 are dropped from
#' the data set.
#' @param by A vector of conditioning variables.
#' @param xlab X-axis label
#' @param ylab Y-axis label
#' @param object An Xpose object. Needed if no \code{data} is supplied.
#' @param events.to.plot Vector of events to be plotted. "All" means that all
#' events are plotted.
#' @param sim.data The simulated data file. Should be a table file with one
#' header row and have, at least, columns with headers corresponding to
#' \code{x}, \code{y}, \code{id}, \code{by} (if used), \code{nsim.lab} and
#' \code{sim.evct.lab}.
#' @param sim.zip.file The \code{sim.data} can be in \.zip format and xpose
#' will unzip the file before reading in the data. Must have the same
#' structure as described above in \code{sim.data}.
#' @param VPC \code{TRUE} or \code{FALSE}. If \code{TRUE} then Xpose will
#' search for a zipped file with name
#' \code{paste("simtab",object@Runno,".zip",sep="")}, for example
#' "simtab42.zip".
#' @param nsim.lab The column header for \code{sim.data} that contains the
#' simulation number for that row in the data.
#' @param sim.evct.lab The column header for \code{sim.data} that contains the
#' individual event counter information. For each individual the event counter
#' should increase by one for each event (or censored event) that occurs.
#' @param probs The probabilities (non-parametric percentiles) to use in
#' computation of the prediction intervals for the simulated data.
#' @param add.baseline Should a (x=0,y=1) baseline measurement be added to each
#' individual in the dataset. Otherwise each plot will begin at the first event
#' in the dataset.
#' @param add.last.area Should an area be added to the VPC extending the last
#' PI?
#' @param subset The subset of the data and sim.data to use.
#' @param main The title of the plot. Can also be \code{NULL} or
#' \code{"Default"}.
#' @param main.sub The title of the subplots. Must be a list, the same length
#' as the number of subplots (actual graphs), or \code{NULL} or
#' \code{"Default"}.
#' @param main.sub.cex The size of the title of the subplots.
#' @param nbins The number of bins to use in the VPC. If \code{NULL}, the the
#' number of unique \code{x} values in \code{sim.data} is used.
#' @param real.se Should the standard errors of the real (non simulated) data
#' be plotted? Calculated using \code{\link[survival]{survfit}}.
#' @param real.se.type Type for the standard errors.
#' @param real.type Type for the real data.
#' @param real.lwd Line width (lwd) for the real data.
#' @param real.lty Line type (lty) for the curve of the original (or real)
#' data.
#' @param real.col Color for the curve of the original (or real) data.
#' @param real.se.lty Line type (lty) for the standard error lines.
#' @param real.se.lwd Line width (lwd) for the standard error lines.
#' @param real.se.col Color for the standard error lines.
#' @param cens.type Type for the censored lines.
#' @param cens.lty Line type (lty) for the censored lines.
#' @param cens.col Color for the censored lines.
#' @param cens.lwd Line width for the censored lines.
#' @param cens.rll The relative line length of the censored line compared to the limits of the y-axis.
#' @param cov The covariate in the dataset to plot instead of the survival
#' curve.
#' @param cov.fun The summary function for the covariate in the dataset to plot
#' instead of the survival curve.
#' @param inclZeroWRES Include WRES=0 rows from the real data set in the plots?
#' @param onlyfirst Include only the first measurement for the real data in the
#' plots?
#' @param samp Simulated data in the xpose data object can be used as the
#' "real" data. \code{samp} is a number selecting which simulated data set to
#' use.
#' @param poly.alpha The transparency of the VPC shaded region.
#' @param poly.fill The fill color of the VPC shaded region.
#' @param poly.line.col The line colors for the VPC region.
#' @param poly.lty The line type for the VPC region.
#' @param censor.lines Should censored observations be marked on the plot?
#' @param ylim Limits for the y-axes
#' @param \dots Additional arguments passed to the function.
#' @return returns an object of class "xpose.multiple.plot".
#' @author Andrew C. Hooker
#' @seealso \code{\link[survival]{survfit}}, \code{\link[survival]{Surv}},
#' \code{\link{xpose.multiple.plot}}.
#' @examples
#'
#' \dontrun{
#' library(xpose4)
#'
#' ## Read in the data
#' runno <- "57"
#' xpdb <- xpose.data(runno)
#'
#' ####################################
#' # here are the real data plots
#' ####################################
#'
#' kaplan.plot(x="TIME",y="DV",object=xpdb)
#' kaplan.plot(x="TIME",y="DV",object=xpdb,
#' events.to.plot=c(1,2),
#' by=c("DOSE==0","DOSE!=0"))
#' kaplan.plot(x="TIME",y="DV",object=xpdb,
#' events.to.plot=c(1,2),
#' by=c("DOSE==0","DOSE==10",
#' "DOSE==50","DOSE==200"))
#'
#' ## make a PDF of the plots
#' pdf(file=paste("run",runno,"_kaplan.pdf",sep=""))
#' kaplan.plot(x="TIME",y="DV",object=xpdb,
#' by=c("DOSE==0","DOSE==10",
#' "DOSE==50","DOSE==200"))
#' dev.off()
#'
#' ####################################
#' ## VPC plots
#' ####################################
#'
#' kaplan.plot(x="TIME",y="DV",object=xpdb,VPC=T,events.to.plot=c(1))
#' kaplan.plot(x="TIME",y="DV",object=xpdb,VPC=T,
#' events.to.plot=c(1,2,3),
#' by=c("DOSE==0","DOSE!=0"))
#' kaplan.plot(x="TIME",y="DV",object=xpdb,VPC=T,
#' events.to.plot=c(1),
#' by=c("DOSE==0","DOSE==10","DOSE==50","DOSE==200"))
#'
#' ## make a PDF of all plots
#' pdf(file=paste("run",runno,"_kaplan.pdf",sep=""))
#' kaplan.plot(x="TIME",y="DV",object=xpdb,VPC=T,
#' by=c("DOSE==0","DOSE==10","DOSE==50","DOSE==200"))
#' dev.off()
#' }
#'
#' @export
#' @family specific functions
kaplan.plot <-
function(x="TIME",y="DV",id="ID",
data= NULL,
evid="EVID",
by=NULL,
xlab="Time",
#ylab="Survival (%)",
ylab="Default",
object=NULL,
events.to.plot="All",
sim.data=NULL,
sim.zip.file=NULL,
VPC = FALSE,
nsim.lab="simNumber",
sim.evct.lab="counter",
probs=c(0.025,0.975),
add.baseline=T,
add.last.area=T,
subset=NULL,
##subset.real=NULL,
main="Default",
main.sub="Default",
main.sub.cex = 0.8,
nbins=NULL,
real.type="l",
real.lty=1,
real.lwd=1,
real.col="blue",
real.se= if(!is.null(sim.data)) F else T,
real.se.type="l",
real.se.lty=2,
real.se.lwd=0.5,
real.se.col="red",
cens.type="l",
cens.lty=1,
cens.col="black",
cens.lwd=1,
cens.rll=0.02,
inclZeroWRES=TRUE,
onlyfirst=FALSE,
#RTTE=FALSE,
samp=NULL,
poly.alpha=1,
poly.fill="lightgreen",
poly.line.col="darkgreen",
poly.lty=2,
censor.lines=TRUE,
ylim=c(-5,105),
cov = NULL,
cov.fun = "mean",
...)
{
## example call
## kaplan.plot(x="TIME",y="DV",object=xpdb,VPC=T,cov="DOSE",cov.fun="mean")
## make max/min of VPC and read data the limits of the x and y-axis
## use latticeExtra C command
##Get data for xpose object
if(is.null(object) && is.null(data)) cat("one of data or object should be defined in function input.")
if(!is.null(data)){
if(!is.null(subset)){
#on.exit(detach(data))
#attach(data,warn.conflicts=F)
#data <- data[eval(parse(text=paste("data$", subset))),]
data<-with(data,data[eval(parse(text=subset)),])
if(dim(data)[1]==0) return(NULL)
}
}
if(!is.null(object)){
if(!is.null(samp)) {
data <- SData(object,inclZeroWRES,onlyfirst=onlyfirst,
subset=subset,samp=samp)
} else {
data <- Data(object,inclZeroWRES,onlyfirst=onlyfirst,subset=subset)
}
}
if (VPC && !is.null(object)) sim.zip.file <- paste("simtab",object@Runno,".zip",sep="")
if (!is.null(sim.zip.file)){
sim.data <- read.table(unz(sim.zip.file,
sub(".zip", "",sim.zip.file)),
skip=0,header=T)
}
if(!is.null(sim.data)){
if(!is.null(subset)){
#on.exit(detach(sim.data))
#attach(sim.data,warn.conflicts=F)
#sim.data <- sim.data[eval(parse(text=paste("sim.data$", subset))),]
sim.data<-with(sim.data,sim.data[eval(parse(text=subset)),])
if(dim(sim.data)[1]==0) return(NULL)
}
}
## Check if this is RTTE or TTE data
data$counter <- 0
counter <- 0
old.id <- 0
for(i in 1:length(data$counter)){
new.id <- data[[id]][i]
if(new.id != old.id){counter <- 0}
old.id <- new.id
if(!is.null(data[[evid]])){ # filter on evid
if (data[[evid]][i]==0) counter=counter+1
} else { # all rows are events
counter=counter+1
}
data$counter[i]=counter
}
events <- max(data$counter)
RTTE = FALSE
if(events>1) RTTE = TRUE
max.events <- events
if(!is.null(sim.data)){
sim.events <- max(sim.data[[sim.evct.lab]])
#if(sim.events>1) RTTE = TRUE
#max.events <- max(sim.events,events)
}
full.data <- data
y.true <- y
if(!is.null(sim.data)){
full.sim.data <- sim.data
## add sim.ID
## add a unique ID identifier to simulated data.
rle.result <- rle(full.sim.data[[id]])
rle.result$values <- 1:length(rle.result$values)
new.id <- inverse.rle(rle.result)
full.sim.data$sim.ID <- new.id
}
#browser()
#by <- c("DOSE==0","DOSE!=0")
by.no <- length(by)
#by.no
if(!all(is.na(match(events.to.plot,"All")))){
num.of.plots <- max.events
if(by.no>0) num.of.plots <- num.of.plots*by.no
plotList <- vector("list",num.of.plots)
event.list <- 1:max.events
} else {
num.of.plots <- length(events.to.plot)
if(by.no>0) num.of.plots <- num.of.plots*by.no
plotList <- vector("list",num.of.plots)
event.list <- events.to.plot
}
plot.num <- 0
for(event.no in event.list){
if(is.null(by)) by=c(1)
for(by.val in by){
#i=1 # for testing
if(by.val==1){
by.val=NULL
by=NULL
}
tmp.cex <- main.sub.cex
if (is.null(main.sub) | num.of.plots==1){
tmp.name <- NULL
} else {
if(!is.na(match(main.sub,"Default"))) {
tmp.name <- paste("Event",event.no,sep=" ")
if(!is.null(by.val)) tmp.name <- paste(tmp.name,by.val,sep=", ")
} else {
tmp.name <- main.sub[plot.num+1]
}
}
data <- full.data
if(!is.null(sim.data)) sim.data <- full.sim.data
data <- subset(data,counter<=event.no)
if(!is.null(data[[evid]])){
data <- data[data[evid]==0,]
}
if(!is.null(sim.data)) sim.data <- subset(sim.data,eval(parse(text=sim.evct.lab))<=event.no)
data <- data[!duplicated(data[[id]], fromLast = TRUE),]
if(!is.null(sim.data)) sim.data <- sim.data[!duplicated(sim.data$sim.ID, fromLast = TRUE),]
## account for censoring and values larger than 1
data$tmp.event <- 0
data$tmp.event[data[y.true]!=0] = 1 # values larger than 1 set to 1
data$tmp.event[data$counter < event.no] = 0 # censored events set to zero
y="tmp.event"
## need to account for simulated data as well
if(!is.null(sim.data)){
sim.data$tmp.event <- 0
sim.data$tmp.event[sim.data[y.true]!=0] = 1 # values larger than 1 set to 1
sim.data$tmp.event[sim.data[[sim.evct.lab]] < event.no] = 0 # censored events set to zero
}
if(!is.null(by.val)){
data <- subset(data,eval(parse(text=by.val)))
#data <- full.data[eval(parse(text=paste("full.data$", by.val))),]
if(dim(data)[1]==0) return(NULL)
if(!is.null(sim.data)){
sim.data <- subset(sim.data,eval(parse(text=by.val)))
#sim.data <- full.sim.data[eval(parse(text=paste("full.sim.data$", by.val))),]
if(dim(sim.data)[1]==0) return(NULL)
}
}
#sim.data$tmp.event <- sim.data[[y]]
##d.sub[c("ID","new.DV","TIME")]
##d.sub$ETMP = 0
#d.sub$ETMP[d.sub$TYPE!=3] = 1
#names(d.sub)
#kaplan.plot(x="TIME",y="DV",data=d.sub)
#browser()
## make new column for when censoring (no event) occurs
## if(!RTTE){
## data$tmp.event <- 0
## data$tmp.event[data[y]!=0] = 1
## sim.data$tmp.event <- sim.data[[y]]
## y="tmp.event"
## }
#browser()
S <- survival::Surv(data[,x],data[,y])
##f.1 <- survfit(S)
f.1 <- survival::survfit(S~1)
a.1 <- summary(f.1)
# plot(f.1)
## f.2 <- survfit(S~data[,"DOSE"])
## plot(f.2,lty=1:4,col=1:4,yscale=100,ylab="Survival (%)", xlab="Time")
## legend(350, 1.0, c("Placebo", "10 mg","50 mg", "100 mg"), lty = 1:4,col=1:4)
##browser()
## ## covariate plots
cov.plot <- FALSE
censored.ids <- FALSE
if(!is.null(cov)) cov.plot <- TRUE
if(cov.plot){
if(is.factor(data[,cov])){
cat("\n Transforming",cov,"from levels to numeric\n",sep=" ")
data[,cov] <- as.numeric(levels(data[,cov]))[data[,cov]]
}
tmp.y.cov <- c()
## the full covariate set
tmp.cov.prev <- data[,cov]
for(i in 1:length(f.1$time)){
if(f.1$n.censor[i]>0 && f.1$n.event[i]==0){ # we have only censored events at this time point
##tmp.cov <- data[,cov][data[,x]>=f.1$time[i]] # Include the censored ID in covariate list
#print("hi\n")
tmp.cov <- tmp.cov.prev
} else { # Include IDs that survive beyond that time in cov list
tmp.cov <- data[,cov][data[,x]>f.1$time[i]]
}
tmp.y.cov <- c(tmp.y.cov,do.call(cov.fun,list(tmp.cov)))
#print(mean(tmp))
tmp.cov.prev <- tmp.cov
}
}
##fit a Kaplan-Meier and plot it
##fit <- survfit(Surv(time, status) ~ x, data = aml)
##plot(fit, lty = 2:3)
##legend(100, .8, c("Maintained", "Nonmaintained"), lty = 2:3)
##fit a Cox proportional hazards model and plot the
##predicted survival for a 60 year old
##fit <- coxph(Surv(futime, fustat) ~ age, data = ovarian)
##plot(survfit(fit, newdata=data.frame(age=60)),
## xscale=365.25, xlab = "Years", ylab="Survival")
if(!is.null(sim.data)){
times <- sort(unique(sim.data[,x]))
m1 <- matrix(nrow=length(unique(sim.data[,nsim.lab])),#number of simulaitons
ncol=length(times)) # number of times
if(cov.plot){
m2 <- m1
}
ii=0
for(i in unique(sim.data[,nsim.lab])){
ii=ii+1
##tmp <- subset(sim.data[,nsim==i],nsim==i)
tmp <- sim.data[eval(parse(text=paste("sim.data$", nsim.lab,"==",i))),]
#if(dim(tmp)[1]==0) browser()
S.sim <- survival::Surv(tmp[,x],tmp[,y])
##f.sim <- survfit(S.sim~tmp$TRT)
f.1.sim <- survival::survfit(S.sim~1)
##a.sim <- summary(f.sim)
a.1.sim <- summary(f.1.sim)
tmp.times <- f.1.sim$time
col.index <- match(tmp.times,times)
## plot the simulated k-m curves
#browser()
#plot(f.1.sim)
if(cov.plot){
if(is.factor(tmp[,cov])){
cat("\n Transforming",cov,"from levels to numeric\n",sep=" ")
tmp[,cov] <- as.numeric(levels(tmp[,cov]))[tmp[,cov]]
}
tmp.y.cov.sim.n <- c()
for(j in 1:length(f.1.sim$time)){
if(f.1.sim$n.censor[j]>0 && f.1.sim$n.event[j]==0){
tmp1.cov <- tmp[,cov][tmp[,x]>=f.1.sim$time[j]]
#print("hi\n")
} else {
tmp1.cov <- tmp[,cov][tmp[,x]>f.1.sim$time[j]]
}
tmp.y.cov.sim.n <- c(tmp.y.cov.sim.n,do.call(cov.fun,list(tmp1.cov)))
#print(mean(tmp1))
}
}
##print(dim(m1))
##print(col.index)
##print(i)
##print(by.val)
## add the new information
m1[ii,col.index] <- f.1.sim$surv
## make sure the total prop is carried through in the times
if(is.na(m1[ii,1])){
m1[ii,1] <- 1
}
for(j in 2:length(times)){
if(is.na(m1[ii,j])){
m1[ii,j] <- m1[ii,j-1]
}
}
if(cov.plot){
## add the new information
m2[ii,col.index] <- tmp.y.cov.sim.n
## make sure the total is carried through in the times
if(is.na(m2[ii,1])){
m2[ii,1] <- do.call(cov.fun,list(tmp[,cov]))
}
for(j in 2:length(times)){
if(is.na(m2[ii,j])){
m2[ii,j] <- m2[ii,j-1]
}
}
}
}
if(cov.plot){
m1 <- m2
}
#browser()
#m1[,1] # this list should actally not contain any NA
## should be 1 or a value.
#dim(m1)
#m1[,129]
if(is.null(nbins)){
quants <- matrix(nrow=2,ncol=dim(m1)[2])
time.bin <- times
for(j in 1:dim(quants)[2]){
quants[,j] <- quantile(m1[,j],probs=probs,na.rm=F)
}
} else {
remainder.times <- 0
ncomb <- floor(dim(m1)[2]/nbins)
if(ncomb==0){
ncomb <- 1
quants <- matrix(nrow=2,ncol=dim(m1)[2])
time.bin <- c()
} else {
remainder.times <- dim(m1)[2]-ncomb*nbins
quants <- matrix(nrow=2,ncol=nbins)
time.bin <- c()
}
k.add.this <- 0
for(j in 1:dim(quants)[2]){
##j <- 1
tmp.bin <- c()
if(remainder.times>0){
ncomb.tmp <- ncomb+1
remainder.times <- remainder.times - 1
} else {
ncomb.tmp <- ncomb
}
for(k in 1:ncomb.tmp){
##k.add.this <- (j-1)*ncomb
if(k==1) time.bin <- c(time.bin,times[k+k.add.this])
if(k+k.add.this>dim(m1)[2]) next
tmp.bin <- c(tmp.bin,m1[,k+k.add.this])
}
k.add.this <- k.add.this+ncomb.tmp
quants[,j] <- quantile(tmp.bin,probs=probs,na.rm=F)
}
}
if(add.baseline){
tmp.mat <- matrix(nrow=dim(quants)[1],ncol=dim(quants)[2]+1)
tmp.mat[,1] <- 1
if(cov.plot) tmp.mat[,1] <- do.call(cov.fun,list(sim.data[,cov]))
tmp.mat[,-1] <- quants
tmp.x <- c(0,time.bin)
} else {
tmp.mat <- quants
tmp.x <- time.bin
}
PI.up <- c()
PI.down <- c()
PI.times <- c()
n.times <- length(tmp.x)
for(i in 1:n.times){
PI.reps=2
time.reps=2
if(i==1) time.reps=1
if(i==n.times) PI.reps=1
PI.up <- c(PI.up,rep(tmp.mat[2,i],PI.reps))
PI.down <- c(PI.down,rep(tmp.mat[1,i],PI.reps))
PI.times <- c(PI.times,rep(tmp.x[i],time.reps))
}
if(add.last.area){
if(tail(tail(PI.times,n=1))==tail(times,n=1)){
time.pt <- (tail(PI.times,n=1)-PI.times[1])*1.01
} else {
time.pt <- (tail(times,n=1)-PI.times[1])*1.01
}
PI.times=c(PI.times,time.pt)
PI.up <- c(PI.up,tail(PI.up,n=1))
PI.down <- c(PI.down,tail(PI.down,n=1))
}
if(!cov.plot){
PI.up <- PI.up*100
PI.down <- PI.down*100
}
} else{
PI.up=NULL
PI.down=NULL
PI.times=NULL
} # end simulation stuff
if(add.baseline){
tmp.y <- c(1,f.1$surv)
tmp.x <- c(0,f.1$time)
tmp.y.upper <- c(1,f.1$upper)
tmp.y.lower <- c(1,f.1$lower)
if(cov.plot) tmp.y.cov <- c(do.call(cov.fun,list(data[,cov])),tmp.y.cov)
} else {
tmp.y <- c(f.1$surv)
tmp.x <- c(f.1$time)
tmp.y.upper <- c(f.1$upper)
tmp.y.lower <- c(f.1$lower)
}
if(cov.plot){
tmp.y <- tmp.y.cov
}
real.data <- c()
real.times <- c()
real.data.upper <- c()
real.data.lower <- c()
for(i in 1:length(tmp.y)){
y.reps=2
x.reps=2
if(i==1) x.reps=1
if(i==length(tmp.y)) y.reps=1
real.data <- c(real.data,rep(tmp.y[i],y.reps))
real.times <- c(real.times,rep(tmp.x[i],x.reps))
real.data.upper <- c(real.data.upper,rep(tmp.y.upper[i],y.reps))
real.data.lower <- c(real.data.lower,rep(tmp.y.lower[i],y.reps))
}
if(!cov.plot){
real.data <- real.data*100
real.data.upper <- real.data.upper*100
real.data.lower <- real.data.lower*100
}
if(cov.plot){
if(!is.null(sim.data)) {
if(!any(grepl("ylim",names(match.call())))){
ylim <- c(min(range(real.data)[1],range(PI.down)[1])*0.9,max(range(real.data)[2],range(PI.up)[2])*1.1)
}
} else {
if(!any(grepl("ylim",names(match.call())))){
ylim <- c(min(range(real.data)[1])*0.9,max(range(real.data)[2])*1.1)
}
}
real.se <- FALSE
#censor.lines <- FALSE
if(!is.na(match(ylab,"Default"))) {
ylab <- paste(cov.fun,cov)
}
cen.x0 <- f.1$time[f.1$n.censor>0]
cen.x1 <- f.1$time[f.1$n.censor>0]
cen.y0 <- tmp.y[match(f.1$time[f.1$n.censor>0],tmp.x)] -(ylim[2] - ylim[1])*cens.rll
cen.y1 <- tmp.y[match(f.1$time[f.1$n.censor>0],tmp.x)] +(ylim[2] - ylim[1])*cens.rll
}
#print(xyplot(real.data~real.times))
xplot <-
xyplot(real.data~real.times,
main=list(tmp.name,cex=tmp.cex),
ylim = ylim,
xlab=xlab,
ylab=if(!is.na(match(ylab,"Default"))){"Survival (%)"}else{ylab},
real.type=real.type,
PI.up=PI.up,PI.down=PI.down,
PI.times=PI.times,
real.data.upper=real.data.upper,
real.data.lower=real.data.lower,
real.se.type=real.se.type,
real.se.lty=real.se.lty,
real.se.col=real.se.col,
f.1=f.1,
...,
panel=function(x,y,PI.up,PI.down,PI.times,
real.data.upper,real.data.lower,
real.se.type=real.se.type,
real.type=real.type,
real.se.lty=real.se.lty,
real.se.col=real.se.col,
f.1=f.1,
...){
if(!is.null(sim.data)){
grid.polygon(c(PI.times,rev(PI.times)),c(PI.up,rev(PI.down)),
default.units="native",
gp=gpar(fill=poly.fill,alpha=poly.alpha,col=poly.line.col,lty=poly.lty)
)
}
panel.xyplot(x,y,type=real.type,lwd=real.lwd,col=real.col,lty=real.lty,...)
if(real.se) panel.xyplot(x,real.data.upper,type=real.se.type,
lty=real.se.lty,
col=real.se.col,
lwd=real.se.lwd,
...)
if(real.se) panel.xyplot(x,real.data.lower,
type=real.se.type,
lty=real.se.lty,
col=real.se.col,
lwd=real.se.lwd,
...)
if(censor.lines){
if(any(f.1$n.censor>0)){
if(cov.plot){
panel.segments(x0=cen.x0,
y0=cen.y0,
x1=cen.x1,
y1=cen.y1,
...)
} else {
# cen.y0 <- tmp.y[match(f.1$time[f.1$n.censor>0],tmp.x)] -(ylim[2] - ylim[1])*0.01
# cen.y1 <- tmp.y[match(f.1$time[f.1$n.censor>0],tmp.x)] +(ylim[2] - ylim[1])*0.01
# browser()
panel.segments(x0=c(f.1$time[f.1$n.censor>0]),
y0=c(f.1$surv[f.1$n.censor>0]*100-(ylim[2] - ylim[1])*cens.rll),
x1=c(f.1$time[f.1$n.censor>0]),
y1=c(f.1$surv[f.1$n.censor>0]*100+(ylim[2] - ylim[1])*cens.rll),
type=cens.type,
lty=cens.lty,
col=cens.col,
lwd=cens.lwd,
...)
}
}
}
}
)
plot.num <- plot.num+1
#if(plot.num>152) browser()
plotList[[plot.num]] <- xplot
}
}
default.plot.title <- "Kaplan-Meier plots"
if(cov.plot) default.plot.title <- paste(cov.fun,cov)
if (num.of.plots==1) default.plot.title <- paste("Kaplan-Meier plot for event",event.list[1])
if (num.of.plots==1 && cov.plot) default.plot.title <- paste(cov.fun, cov, "for event",event.list[1])
if(!is.null(object)){
plotTitle <- xpose.multiple.plot.title(object=object,
plot.text = default.plot.title,
main=main,
...)
} else {
if (is.null(main)){
plotTitle <- NULL
} else {
if(!is.na(match(main,"Default"))) {
plotTitle <- default.plot.title
if (!is.null(subset)){
plotTitle <- paste(plotTitle,"\n[",subset,"]",sep="")
}
} else {
plotTitle <- main
}
}
}
obj <- xpose.multiple.plot(plotList,plotTitle,...)
return(obj)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.