########################################################### ## A function for making pseudo-values ############### ## for cumulative incidence estimator ##################### ########################################################### pseudoci <- function(time, status, tmax){ library(cmprsk) howmany <- length(time) ## Cumulative incidence function for the whole sample modelfit <- cuminc(time, status) ## The values at time points tmax CI.est <- timepoints(modelfit, times = tmax)\$est ##Pseudoobservations pseudo <- as.data.frame(matrix(nrow=length(time), ncol=2 + 2*length(tmax))) pseudo[,1] <- time pseudo[,2] <- status npoints<-length(tmax) nam<-c("time", "status","rpseudo1","tpseudo1") for(j in 2:npoints){nam1<-c(paste("rpseudo",j,sep=""),paste("tpseudo",j,sep="")) nam<-cbind(nam,nam1)} names(pseudo)<-nam for (j in 1:howmany){ ##CI for sample without the j'th observation modelfitj <- cuminc(time[-j],status[-j]) ## The value of CI at time tmax CI.est.j <- timepoints(modelfitj, times = tmax)\$est ## Pseudoobservation pseudo[j,3:(2+2*length(tmax))] <- howmany*CI.est - (howmany-1)*CI.est.j } detach(package:cmprsk) return(pseudo) } ########################################################### ## A function for making pseudo-values ############### ## for tmax-year survival via the Kaplan-Meier estimator ## ########################################################### pseudosurv <- function(time,cens,tmax){ if (sum(tmax > max(time)) > 0) stop ("tmax greater than largest observation time") else{ library(survival) howmany <- length(time) ##KM for the whole sample kmfit <- survfit(Surv(time,cens)~1) tmax <- sort(tmax) S.of.tmax <- summary(kmfit, times = tmax)\$surv ##Pseudoobservations pseudo <- as.data.frame(matrix(data = NA, ncol = length(tmax)+2, nrow = howmany)) pseudo[,1] <- time pseudo[,2] <- cens names(pseudo) <- c("time", "cens",paste("tmax =", round(tmax,3), sep="")) for (j in 1:howmany){ ##KM for sample without the j'th observation kmfit <- survfit(Surv(time[-j],cens[-j])~1) S.of.tmax.j <- summary(kmfit, times = tmax)\$surv ## Pseudoobservation pseudo[j,-c(1,2)] <- howmany*S.of.tmax-(howmany-1)*S.of.tmax.j } return(pseudo) detach("package:survival") } } ########################################################### ## A function for making pseudo-values ############### ## for restricted mean #################################### ########################################################### pseudomean <- function(time, dead, tmax){ library(survival) rmtime <- ifelse(time >= tmax ,tmax,time) rmdead <- ifelse(time >= tmax ,0, dead) howmany <- length(time) ##KM for rm-sample kmfit <- survfit(Surv(rmtime,rmdead)~1)#km-estimator difft <- diff(c(0, summary(kmfit,censored=FALSE)\$time,tmax)) #dist. between jumps prevsurv <- c(1,summary(kmfit,censored=FALSE)\$surv)#survival values at prev. jump ## Mean value rmean <- sum(difft*prevsurv)#area below KM ##Pseudoobservations psumean <- c(rep(0,howmany)) for (j in 1:howmany){ ##KM for sample without the j'th observation kmfitj <- survfit(Surv(rmtime[-j],rmdead[-j])~1) kmfit <- survfit(Surv(rmtime,rmdead)~1) difftj <- diff(c(0, summary(kmfitj,censored=FALSE)\$time,tmax)) prevsurvj <- c(1,summary(kmfitj,censored=FALSE)\$surv) ## Mean value rmeanj <- sum(difftj*prevsurvj) ## Pseudoobservation psumean[j] <- howmany*rmean-(howmany-1)*rmeanj } return(as.data.frame(cbind(time, dead, psumean))) detach("package:survival") }