R/lifetable.R

Defines functions LifeTable lifetable.formula lifetable.matrix

Documented in lifetable.formula lifetable.matrix

##' @export
`lifetable` <- function(x,...) UseMethod("lifetable")

##' Create simple life table 
##' 
##' @title Life table
##' @param x time formula (Surv) or matrix/data.frame with columns time,status or entry,exit,status
##' @param strata Strata
##' @param data data.frame
##' @param breaks Time intervals
##' @param confint If TRUE 95\% confidence limits are calculated
##' @param ... Additional arguments to lower level functions
##' @author Klaus K. Holst
##' @aliases lifetable lifetable.matrix lifetable.formula
##' @usage
##'  \method{lifetable}{matrix}(x, strata = list(), breaks = c(),
##'    confint = FALSE, ...)
##'
##'  \method{lifetable}{formula}(x, data=parent.frame(), breaks = c(),
##'    confint = FALSE, ...)
##' @examples
##' library(timereg)
##' data(TRACE)
##' \donttest{
##'     lifetable(Surv(time,status==9)~sex+I(cut(wmi,c(-Inf,1,1.5,Inf))),
##'               data=TRACE,breaks=c(0.2),confint=TRUE)
##' }
##' 
##' d <- with(TRACE,lifetable(Surv(time,status==9)~sex+vf,breaks=c(0,0.2,0.5,8.5)))
##' summary(glm(events ~ offset(log(atrisk))+factor(int.end)*vf + sex*vf,
##'             data=d,poisson))
##' @export
lifetable.matrix <- function(x,strata=list(),breaks=c(),confint=FALSE,...) {
    if (ncol(x)==3) {
        status <- x[,3]
        entry <- x[,1]
        time <- x[,2]
    } else {
        status <- x[,2]
        time <- x[,1]
        entry <- rep(0,length(time))
    }
    LifeTable(time,status,entry,strata,breaks,confint,...)
}

##' @export
lifetable.formula <- function(x,data=parent.frame(),breaks=c(),confint=FALSE,...) {
    cl <- match.call()
    mf <- model.frame(x,data)
    Y <- model.extract(mf, "response")
    Terms <- terms(x, data = data)
    if (!is.Surv(Y)) stop("Expected a 'Surv'-object")
    if (ncol(Y)==2) {
        exit <- Y[,1]
        entry <- NULL ## rep(0,nrow(Y))
        status <- Y[,2]
    } else {
        entry <- Y[,1]
        exit <- Y[,2]
        status <- Y[,3]
    }
    strata <- list()
    X <- model.matrix(Terms, data)
    if (!is.null(intpos  <- attributes(Terms)$intercept))
        X <- X[,-intpos,drop=FALSE]
    if (ncol(X)>0) {
        strata <- as.list(model.frame(Terms,data)[,-1,drop=FALSE])
    }
    LifeTable(exit,status,entry,strata,breaks,confint,...)       
}


LifeTable <- function(time,status,entry=NULL,strata=list(),breaks=c(),confint=FALSE) {
    if (is.null(entry)) entry <- rep(0,NROW(time))
    if ((is.matrix(time) || is.data.frame(time)) && ncol(time)>1) {
        if (ncol(time)==3) {
            status <- time[,3]
            entry <- time[,1]
            time <- time[,2]
        } else {
            status <- time[,2]
            time <- time[,1]
            entry <- rep(0,length(time))
        }
    }
    if (length(strata)>0) {
        a <- by(cbind(entry,time,status), strata,
                FUN=LifeTable, breaks=breaks, confint=confint)
        cl <- lapply(strata,class)
        nulls <- which(unlist(lapply(a,is.null)))
        nn <- do.call("expand.grid",attributes(a)$dimnames)
        if (length(nulls)>0) nn <- nn[-nulls,,drop=FALSE]
        nam <- nn[rep(seq(NROW(nn)),each=NROW(a[[1]])),,drop=FALSE]
        xx <- list()
        for (i in seq(ncol(nam))) {
            if (cl[i]%in%c("numeric","integer"))
                xx <- c(xx,list(as.numeric(as.character(nam[,i]))))
            else                    
                xx <- c(xx, list(do.call(paste("as.",as.character(cl[i]),sep=""),list(nam[,i]))))
        }
        xx <- as.data.frame(xx); colnames(xx) <- colnames(nam)
        res <- Reduce("rbind",a)
        res <- cbind(res,xx)
        return(res)
    }
    if (length(breaks)==0) breaks <- c(0,max(time,na.rm=TRUE))
    if (length(breaks)==1) breaks <- c(0,breaks)
    breaks <- sort(unique(breaks))
    ## en <- matrix(unlist(lapply(c(-Inf, breaks),function(x) pmax(x,entry))),
    ##               ncol=length(breaks)+1)
    ## ex <- matrix(unlist(lapply(c(breaks, Inf),function(x) pmin(x,time))),
    ##               ncol=length(breaks)+1)
    en <- matrix(unlist(lapply(breaks[-length(breaks)],function(x) pmax(x,entry))),
                  ncol=length(breaks)-1)
    ex <- matrix(unlist(lapply(breaks[-1],function(x) pmin(x,time))),
                  ncol=length(breaks)-1)
    dur <- ex-en
    endur <- rbind(breaks[-1])%x%cbind(rep(1,nrow(en)))-en
    ##endur <- rbind(c(breaks,Inf))%x%cbind(rep(1,nrow(en)))-en
    dur[dur<0] <- NA
    enter <- colSums(!is.na(dur))
    atrisk <- colSums(dur,na.rm=TRUE)
    eventcens <- dur<endur
    eventcens <- rbind(apply(dur<endur,2,function(x) x*(status+1)))
    lost <- colSums(eventcens==1,na.rm=TRUE)
    events <- colSums(eventcens==2,na.rm=TRUE)
    res <- subset(data.frame(enter=enter,
                             atrisk=atrisk,
                             lost=lost,
                             events=events,
                             ## int.start=c(-Inf,breaks),
                             ## int.end=c(breaks,Inf),
                             int.start=breaks[-length(breaks)],
                             int.end=breaks[-1],
                             surv=0,
                             rate=events/atrisk))
    cumsum.na <- function(x,...) { x[is.na(x)] <- 0; cumsum(x) }
    res$surv <- with(res, exp(-cumsum.na(rate*(int.end-int.start))))    
    if (confint) {
        ff <- events ~ offset(log(atrisk))
        if (length(breaks)>2) ff <- update(ff,.~.+factor(int.end)-1)
        g <- glm(ff,data=res,poisson)
        suppressMessages(ci <- rbind(exp(stats::confint(g))))
        res[,"2.5%"] <- ci[,1]
        res[,"97.5%"] <- ci[,2]
    }
    res
}

Try the mets package in your browser

Any scripts or data that you put into this service are public.

mets documentation built on May 31, 2017, 1:52 a.m.