R/anova.coxph.R

In therneau/survival: Survival Analysis

# The anova function for a coxph object
anova.coxph <- function (object, ...,  test = 'Chisq') {
    if (!inherits(object, "coxph"))
        stop ("argument must be a cox model")

    # All the ... args need to be coxph or coxme fits.  If any of them
    #  have a name attached, e.g., 'charlie=T' we assume a priori
    #  that they are illegal
    #
    dotargs <- list(...)
    named <- if (is.null(names(dotargs))) 
	           rep(FALSE, length(dotargs))
             else (names(dotargs) != "")
    if (any(named)) 
        warning(paste("The following arguments to anova.coxph(..)", 
            "are invalid and dropped:", paste(deparse(dotargs[named]), 
                collapse = ", ")))
    dotargs <- dotargs[!named]
    
    single <- (inherits(object, "coxph") || inherits(object, "coxme"))
    if (length(dotargs) >0  || !single) {
        # there are multiple arguments, either the object itself is a list
        #  of models, or there were multiple arguments.
        # paste them all together into a single list
        if (single) object <- list(object)
        if (length(dotargs)>0) object <- c(object, dotargs)

        # coxme and coxphms models get sent elsewhere
        is.coxme <-  sapply(object, function(x) inherits(x, "coxme"))
        is.multi <-  sapply(object, function(x) inherits(x, "coxphms"))
        if (any(is.multi)) return(anova.coxphms(object, test=test))
        
        if (any(is.coxme)) {
            # We need the anova.coxmelist function from coxme
            # If coxme is not loaded the line below returns NULL
            temp <- getS3method("anova", "coxmelist", optional=TRUE)
            if (is.null(temp)) 
                stop("a coxme model was found and library coxme is not loaded")
            else return(temp(object, test = test))
        }
        else return(anova.coxphlist(object, test = test))
    }

    #
    # The argument is a single Cox model 
    # Show the nested list of models generated by this model.
    # By tradition the sequence is main effects (in the order found in
    #  the model statement), then 2 way interactions, then 3, etc.
    #  One does this by using the "assign" attribute of the model matrix.
    #  (This does not work for penalized terms.)
    # Remember to propogate any the method argument
    mtie <- object$method
    if (inherits(object, "coxphms")) return(anova.coxphms(object, test=test))
    if (length(object$rscore)>0)
        stop("Can't do anova tables with robust variances")
    
    has.strata <- !is.null(attr(terms(object), "specials")$strata)
    if (is.null(object[['y']]) || (has.strata && is.null(object$strata))) {
        # We need the model frame
        mf <- stats::model.frame(object)
        Y <- model.response(mf)
        X <- model.matrix(object, mf)
        if (has.strata) {
            stemp <- untangle.specials(terms(object), "strata")
            if (length(stemp$vars)==1) strata.keep <- mf[[stemp$vars]]
            else strata.keep <- strata(mf[,stemp$vars], shortlabel=TRUE)
            strats <- as.numeric(strata.keep)
        }
    }
    else {
        Y <- object[['y']]
        X <- model.matrix(object)
        if (has.strata) strats <- object$strata
    }

    assign <- attr(X, 'assign') 
    alevels <- sort(unique(assign)) #if there are strata the sequence has holes
    nmodel <- length(alevels)

    df <- integer(nmodel+1)  #this will hold the df vector
    loglik <- double(nmodel+1) #and the loglike vector
    df[nmodel+1] <- if (is.null(object$df)) sum(!is.na(object$coefficients))
                    else sum(object$df)
    loglik[nmodel+1] <- object$loglik[2]
    df[1] <- 0
    loglik[1] <- object$loglik[1]
    
    # Now refit intermediate models
    for (i in seq.int(1, length.out =nmodel-1)){
        if (length(object$offset)) {
            if (has.strata) 
                tfit <- coxph(Y ~ X[,assign <= alevels[i]] + strata(strats) +
                              offset(object$offset), ties=mtie)
            else tfit <- coxph(Y ~ X[, assign<= alevels[i]] +
                               offset(object$offset), ties=mtie)
        }
        else {
            if (has.strata) 
                tfit <- coxph(Y ~ X[,assign <= alevels[i]] + strata(strats),
                              ties=mtie)
            else tfit <- coxph(Y ~ X[,assign <= alevels[i]], ties=mtie)
        }
        df[i+1] <- sum(!is.na(tfit$coefficients))
        loglik[i+1] <- tfit$loglik[2]
    }
                              
    table <- data.frame(loglik=loglik, Chisq=c(NA, 2*diff(loglik)), 
                        Df=c(NA, diff(df))) 

    if (nmodel == 0) #failsafe for a NULL model
             table <- table[1, , drop = FALSE]

    if (length(test) >0 && test[1]=='Chisq') {
        table[['Pr(>|Chi|)']] <- pchisq(table$Chisq, table$Df, lower.tail=FALSE)
        }
  
    temp <- attr(terms(object), "term.labels")
    if (has.strata && length(stemp$terms)>0 ) temp <- temp[-stemp$terms]
    row.names(table) <- c('NULL', temp)

    title <- paste("Analysis of Deviance Table\n Cox model: response is ",
		   deparse(object$terms[[2]]),
		   "\nTerms added sequentially (first to last)\n", 
		   sep = "")
    structure(table, heading = title, class = c("anova", "data.frame"))
}

anova.coxphms <- function(object, ..., test="chisq") {
    stop("anova method not yet available for multi-state coxph fits")
}