R/fromtimeline.R

In survival: Survival Analysis

#
# Routine to turn a "timeline" model frame into a "counting process" model
#  frame. 
# The mf arg is most often the result of a call to coxph or survfit (or a
#  survfit.coxphms or residuals.coxphms call that had to recreate the data), and
#  that fcn is using this function to transform to CP style. As such, there will 
#  normally be a column named "(id)", and perhaps "(istate)" and/or "(cluster)".
# The task is fairly simple:
#   1. An id with 10 rows will have 9 rows in the new data set, 2-10 
#     contribute the (time2, outcome) part of the survival response, and
#     original rows 1-9 everything else.
#   2. Missing covariates are filled in using last-value-carried forward.
#   3. If there is no '(istate)', and normally there won't be because a
#     coxph call doesn't need that argument for Surv2 data, then create it.
#
# The repeat attribute of the Surv2 object determines whether events can
#   "stutter", i.e., two of the same type that are adjacent, or only have
#    censored times between them, count as a second event.
#
surv2counting <- function(mf, repeated=FALSE, lvcf=TRUE) {
    Terms <- terms(mf)
    y <- model.response(mf)
    n <- nrow(y)
    states <- attr(y, "states")
    # the next line supports Surv2 objects
    if (missing(repeated) && !is.null(attr(y, "repeated"))) 
        repeated <- attr(y, "repeated")

    id <- model.extract(mf, "id")
    if (length(id) != n) stop("id statement is required")

    # relax this some later day (or not?)
    if (any(is.na(id)) || any(is.na(y[,1])))
        stop("id and time cannot be missing")
    
    # The data isn't necessarity sorted, don't reorder it
    isort <- order(id, y[,1])
    id2 <- id[isort]
    y2  <- y[isort,]
    first <- !duplicated(id2)
    last  <- !duplicated(id2, fromLast=TRUE)
    # check for duplcate times
    temp <- data.frame(id2, y2[,1]) 
    if (any(duplicated(temp)))
        stop("duplicated time for an id")
        
    mf2 <- mf[isort,][!last,]  # starting point for the output

    # The LVCF operation replaces any missing values with the most recent non
    #  missing.  This is easiest to do with the tmerge3 C routine, which expects
    #  the data to be in time within id order
    # Use LVCF on all colums except: doesn't have any missings (not needed),
    #  column 1= response, and those with () names, normally (id), (cluster)
    id3 <- id2[!last]  # I need idi again later
    if (is.factor(id3)) idi <- as.integer(id3)
    else    idi <- match(id3, unique(id3))  # tmerge3 wants an integer id

    if (lvcf) {
        skip <- c(1, which(substring(names(mf), 1, 1) == "("))
        for (i in (1:ncol(mf2))[-skip]) {
            miss <- is.na(mf2[[i]])
            if (any(miss)) {
                k <- .Call(Ctmerge3, idi, miss)
                # k is the row number of the value to be carried forward
                update <- (miss & k>0)  # some values will be updated
                if (any(update))  { # some successful replacements
                    if (is.matrix(mf2[[i]])) 
                        mf2[[i]][update,] <-mf2[[i]][k[update],]
                    else mf2[[i]][update] <- (mf2[[i]])[k[update]]
                }
            }
        }
    }
    
    # If there are missing status values, assume that NA is actually
    #  the code for censored.  This is not uncommon when merge() is
    #  used to create a data set
    if (any(is.na(y2[,2]))) {
        if (is.null(states)) {
            # Really? The user had NA, 1, 2 and didn't use a factor?
            # I'm going to guess not and leave the NA rows to be removed
            # later by na.action
        } else { 
            # input was a factor
            # it's the user's job to make sure that the first level of that
            #   factor is "censored"
            y2[,2] <- ifelse(is.na(y2[,2]), 0, y2[,2])
        }
    }

    if (is.character(repeated) && casefold(repeated) == "first") {
        # only allow the first instance of any outcome. We need to do
        # this before computing istate2, since it may change the definition
        # of 'current state'
        temp <- unlist(tapply(y2[,2], id2, function(x) {
            ifelse(x==0, 0, ifelse(duplicated(x),0,x))
            }))
        y2[,2] <- temp
    }

    # Create a current state vector, which is just LVCF using the state.
    # NA or censored  treated as "not present" for the tmerge3 call.  
    # Y[,2]=0 is censored for either single state or multistate timeline
    y3 <- y2[!last]
    censored <- ifelse(is.na(y2[,2]), TRUE, y2[,2] ==0) # sorted version
    if (all(censored[first])) {
        # special case -- no initial state for anyone, don't create istate
        # this occurs for ordinary 0/1 status, or perhaps competing risks
        istate2 <- NULL
    } 
    else if (any(censored[first]))
        stop("everyone or no one should have an initial state")
    else {
        istate2 <- y2[!last,2]
        censored <- censored[!last]
        if (any(censored)) { # 0 censored rows is possible, but rare
            k <- .Call(Ctmerge3, idi, censored) # replace censored rows
            istate2[censored] <- istate2[k[censored]]
        }
        mf2[, "(istate)"] <- factor(istate2, 1:length(states), states)

        istate <- model.extract(mf, "istate")
        if (!is.null(istate)) {
            # The user had an istate= argument in their multistate (coxph) or 
            #  Aalen-Johansen (survfit) call.  It should agree, perfectly.  
            # Well, they could legally put the states in a different order
            xx <-"istate argument does not agree with initial Surv2 values"
            if (is.numeric(istate) && any(istate2 != istate)) stop(xx)
            else {
                istate <- as.factor(istate)  # in case is is type character
                itemp <- match(levels(istate), states, nomatch=0)
                if (any(itemp==0)) stop(xx)
                if (any(istate2[itemp] != as.numeric(istate))) stop(xx)
            }
        }            
    }
    # Now replace the Surv2 response with one of the standard forms
    #    "right"   :  (time, status) 
    #    "counting":  (time1, time2, status)
    #    "mright"  :  (time, state)
    #    "mcounting": (time1, time2, state)
    status <- y2[!first, 2] # integer version
    if (is.logical(repeated)) {
        if (!repeated  && !is.null(istate2)) 
            status[status== istate2] <- 0 # no repeat of current state
    } else {
        if (!is.character(repeated) || casefold(repeated) != "first") 
            stop("invalid value for repeated option")
    }
        
    if (!is.null(states)) 
        status <- factor(status, 0:length(states), c("censor", states))

    tstart <- y2[!last, 1]
    tstop  <- y2[!first,1]
    # right or mright
    if (!any(duplicated(id3)))  # one obs per id
        mf2[[1]] <- Surv(tstop, status)
    else
        mf2[[1]] <- Surv(tstart, tstop, status)

    #put the data back into the original order
    jj <- order(isort[!last])  # this line is not obvious, but it works!
    mf2[jj,]
}
   

fromtimeline <- function(formula, data, subset, id, repeated= FALSE,
                              lvcf = TRUE, yname=c("tstart", "tstop", "status")){
    # Get the formula for the response
    Call <- match.call()
    indx <- match(c("formula", "data", "subset", "id"),
                  names(Call), nomatch=0)
    if (indx[1] ==0) stop("a formula argument is required")
    if (indx[2] ==0) stop("the data argument is required")
    if (indx[4] ==0) stop("the id argument is required")
    tform <- Call[c(1, indx)]  # only keep arguments we wanted
    tform$na.action <- stats::na.pass
    tform[[1L]] <- quote(stats::model.frame)
    mf <- eval(tform, parent.frame())
    
    Y <- model.response(mf)
    id <- model.extract(mf, "id")
    if (!(is.Surv(Y) || inherits(Y, "Surv2")))
        stop("response must be a survival object")
    if (is.Surv(Y) && attr(Y, "type") %in% c("counting", "mcounting"))
        stop("response cannot be of counting process type")
    if (is.null(id) || !any(duplicated(id)))
        stop("data does not appear to be timeline data")
    new <- surv2counting(mf, repeated=repeated, lvcf=lvcf)

    # Remove (id) and rename (istate) to istate
    ii <- match("(istate)", names(new)) # added by surv2counting
    if (!is.na(ii)) { # not added in some cases
        if (!any(names(mf)== "istate")) names(new)[ii] <- "istate"
    }
    new["(id)"] <- NULL

    #
    # Give new names to response, so that the data set is more useable
    #
    yy <- new[[1]] # the response is always variable 1
    states <- attr(yy, "states")
    ny <- ncol(yy)  # ny=2 for competing risks, for instance
    if (is.null(states)) status <- yy[,ny]     
    else status <- factor(yy[,ny], 0:length(states), c("censor", states))
    if (ncol(yy) ==3) 
        tdata <- data.frame(tstart= yy[,1], tstop = yy[,2], status)
    else tdata <- data.frame(tstart=yy[,1], status)

    if (!missing(yname)) {
        if (any(!is.na(match(yname, names(new)))))
            stop("element of yname conflicts with an existing name in the data")
        if (ny==2) {
            if (!(length(yname) %in% 2:3)) stop("wrong length for yname")
            names(tdata) <- c(yname[1], yname[length(yname)])
        } else {
            if (length(yname) != 3) stop("wrong length for yname") 
            names(tdata) <- yname
        }
    } else {
        # Use a best guess, if possible
        resp <- formula[[2]] # "Surv(age, dead)", say
        if (is.name(resp[[2]])) {
            if (ny==2) yname <- c(as.character(resp[2]), yname[3])
            else yname[1:2] <- paste0(as.character(resp[2]), 1:2)
        }
        if (is.name(resp[[3]])) yname[length(yname)] <- as.character(resp[[3]])
        
        conflict <- match(yname, names(new), nomatch=0)
        if (any(conflict >0))
            yname[conflict>0] <- paste0("_", yname[conflict>0], "_")
        names(tdata) <- yname
    }
    cbind(new[,-1,drop=FALSE], tdata)
}

# This function has not been tested (note the browser call), it would
#  be the basis of a Surv(time1, time2, stat) to Surv2 transform function
# We await any definite proof of need before working more
# Until tested, it won't appear in the NAMESPACE file or the man pages

totimeline <- function(formula, data, id, istate) {
    if (missing(formula) || missing(data) || missing(id))
        stop("formula, data, and id arguments are required")
    Call <- match.call()
    tcall <- Call
    tcall[[1L]] <- quote(stats::model.frame)

    mf <- eval(tcall, parent.frame())

    Y <- model.response(mf)
    id <- model.extract(mf, "id")

    if (!inherits(Y, "Surv")) stop("response must be a Surv object")
    type <- attr(Y, "type")
    if (type== "left" || type=="interval")
        stop("not valid for interval censored or left censored data")
    if (ncol(Y) != 3) stop("initial data is not in (time1, time2) form")
    
    # get a name for the resulting time and state variables, by parsing
    #  the formula.  Don't try too hard, though.
    tname <- "(time)"
    sname <- "(state)"  # backup defaults
    lhs <- formula[[2]]
    if ((is.name(lhs[[1]]) && lhs[[1]]== as.name("Surv")) ||
        deparse(lhs[[1]]) == "survival::Surv")  {
        # the lhs if of length 4
        if (is.name(lhs[[3]])) tname <- as.character(lhs[[3]])
        if (is.name(lhs[[4]]))  sname <- as.character(lhs[[4]])
        else if (is.call(lhs[[4]])){
            temp <- lhs[[4]]
            if (deparse(temp[[1]])== "factor" && is.name(temp[[2]]))
                sname <- as.character(temp[[2]])
        }
    }   
    
    if (is.name(Call$id)) {
        idname <- as.character(Call$id)
        i <- match(idname, names(mf))
        if (is.na(i)) names(mf)[match("(id)", names(mf))] <- idname
    } else stop("id must be a simple variable name")

    # Get the list of states and istate
    if (missing(istate)) check <- survcheck2(Y, id)
    else {
        istate <- model.extract(mf, "istate")
        check <- survcheck2(Y, id, istate)
        }
    if (any(check$states == "censor")) states <- c("(censor)", check$states)
    else states <- c("censor", check$states)
    nstate <- length(check$states)
    # In the new data, there is 1 more row per subject.
    # newtime for subject i = c(first time1, time2) 
    # newstate is c(initial, Y[i,3])
    # for covariates, the last row of each subject is a repeat
    first <- !duplicated(id)
    last  <- !duplicated(id, fromLast=TRUE)
    n <- nrow(mf)

    indx1 <- rep(1:n, ifelse(first, 2, 1))
    indx2 <- rep(1:n, ifelse(last,  2, 1))
    newtime <- Y[indx1,2]
    newstat <- c(0L, match(attr(Y, "states"), check$states))[1L+ Y[indx1,3]]
    browser()
    row1 <- duplicated(indx1, fromLast=TRUE) # first row of each subject
    newtime[row1] <- Y[first, 1]
    newstat[row1] <- as.numeric(check$istate[first])
    newdata <- cbind(data.frame("(time)"= newtime, 
                                "(state)"= factor(newstat, 0:nstate, states)),
                     mf[indx2,-1])
    row.names(newdata) <- NULL  # they are annoying and useless
    names(newdata)[1:2] <- c(tname, sname)
    indx <- match(c("(id)", "(istate)"), names(newdata), nomatch=0)
    newdata[, -indx]  # remove the redundant (id) and (istate) columns
}