The R Project SVN R-packages

#SCCS  @(#)coxph.s  5.12 06/12/00
# Version with general penalized likelihoods

coxph <- function(formula=formula(data), data=parent.frame(),
    weights, subset, na.action, init,
    control, method= c("efron", "breslow", "exact"),
    singular.ok =TRUE, robust=FALSE,
    model=FALSE, x=FALSE, y=TRUE, ...) {

    method <- match.arg(method)
    call <- match.call()
    m <- match.call(expand.dots=FALSE)
    temp <- c("", "formula", "data", "weights", "subset", "na.action")
    m <- m[ match(temp, names(m), nomatch=0)]
    special <- c("strata", "cluster")
    Terms <- if(missing(data)) terms(formula, special)
         else              terms(formula, special, data=data)
    m$formula <- Terms
    m[[1]] <- as.name("model.frame")
    m <- eval(m, parent.frame())
    if (NROW(m)==0)
        stop("No (non-missing) observations")
    
    if (missing(control)) control <- coxph.control(...)
    Y <- model.extract(m, "response")
    if (!inherits(Y, "Surv")) stop("Response must be a survival object")
    weights <- model.extract(m, 'weights')
    offset<- attr(Terms, "offset")
    tt <- length(offset)
    offset <- if(tt == 0)
            rep(0, nrow(Y))
          else if(tt == 1)
              m[[offset]]
          else {
            ff <- m[[offset[1]]]
            for(i in 2:tt)
                ff <- ff + m[[offset[i]]]
            ff
            }

    attr(Terms,"intercept")<- 1  #Cox model always has \Lambda_0
    strats <- attr(Terms, "specials")$strata
    cluster<- attr(Terms, "specials")$cluster
    dropx <- NULL
    if (length(cluster)) {
    if (missing(robust)) robust <- TRUE
    tempc <- untangle.specials(Terms, 'cluster', 1:10)
    ord <- attr(Terms, 'order')[tempc$terms]
    if (any(ord>1)) stop ("Cluster can not be used in an interaction")
    cluster <- strata(m[,tempc$vars], shortlabel=TRUE)  #allow multiples
    dropx <- tempc$terms
    }
    if (length(strats)) {
    temp <- untangle.specials(Terms, 'strata', 1)
    dropx <- c(dropx, temp$terms)
    if (length(temp$vars)==1) strata.keep <- m[[temp$vars]]
    else strata.keep <- strata(m[,temp$vars], shortlabel=TRUE)
    strats <- as.numeric(strata.keep)
    }

    ##if (length(dropx)) X <- model.matrix(Terms[-dropx], m)[,-1,drop=F]
    ##else               X <- model.matrix(Terms, m)[,-1,drop=F]
    ### this is inefficient, but subscripting loses the assign attribute
    if (length(dropx))
        newTerms<-Terms[-dropx]
    else
        newTerms<-Terms
    X<-model.matrix(newTerms,m)
    assign<-lapply(attrassign(X,newTerms)[-1],function(x) x-1)
    X<-X[,-1,drop=FALSE]
    
    
    type <- attr(Y, "type")
    if (type!='right' && type!='counting')
    stop(paste("Cox model doesn't support \"", type,
              "\" survival data", sep=''))
    if (missing(init)) init <- NULL

    # Check for penalized terms
    pterms <- sapply(m, inherits, 'coxph.penalty')
    if (any(pterms)) {
    pattr <- lapply(m[pterms], attributes)
    # 
    # the 'order' attribute has the same components as 'term.labels'
    #   pterms always has 1 more (response), sometimes 2 (offset)
    # drop the extra parts from pterms
    temp <- c(attr(Terms, 'response'), attr(Terms, 'offset'))
    if (length(dropx)) temp <- c(temp, dropx+1)
    pterms <- pterms[-temp]
    temp <- match((names(pterms))[pterms], attr(Terms, 'term.labels'))
    ord <- attr(Terms, 'order')[temp]
    if (any(ord>1)) stop ('Penalty terms cannot be in an interaction')
    ##pcols <- (attr(X, 'assign')[-1])[pterms]  
        pcols<-assign[pterms]
        
    #penalized are hard sometimes   
    if (control$eps.miss)   control$eps <- 1e-7
    if (control$iter.miss)  control$iter.max <- 20  

        fit <- coxpenal.fit(X, Y, strats, offset, init=init,
                control,
                weights=weights, method=method,
                row.names(m), pcols, pattr, assign)
    }
    else {
    if( method=="breslow" || method =="efron") {
        if (type== 'right')  fitter <- get("coxph.fit")
        else                 fitter <- get("agreg.fit")
        }
    else if (method=='exact') fitter <- get("agexact.fit")
    else stop(paste ("Unknown method", method))

    fit <- fitter(X, Y, strats, offset, init, control, weights=weights,
                method=method, row.names(m))
    }

    if (is.character(fit)) {
    fit <- list(fail=fit)
    class(fit) <- 'coxph'
    }
    else {
    if (!is.null(fit$coefficients) && any(is.na(fit$coefficients))) {
       vars <- (1:length(fit$coefficients))[is.na(fit$coefficients)]
       msg <-paste("X matrix deemed to be singular; variable",
               paste(vars, collapse=" "))
       if (singular.ok) warning(msg)
       else             stop(msg)
       }
    fit$n <- nrow(Y)
    class(fit) <- fit$method
    fit$terms <- Terms
    ##fit$assign <- attr(X, 'assign')
        fit$assign<-assign
    if (robust) {
        fit$naive.var <- fit$var
        fit$method    <- method
        # a little sneaky here: by calling resid before adding the
        #   na.action method, I avoid having missings re-inserted
        # I also make sure that it doesn't have to reconstruct X and Y
        fit2 <- c(fit, list(x=X, y=Y, weights=weights))
        if (length(strats)) fit2$strata <- strata.keep
        if (length(cluster)) {
        temp <- residuals.coxph(fit2, type='dfbeta', collapse=cluster,
                      weighted=TRUE)
        # get score for null model
        if (is.null(init))
            fit2$linear.predictors <- 0*fit$linear.predictors
        else fit2$linear.predictors <- c(X %*% init)
        temp0 <- residuals.coxph(fit2, type='score', collapse=cluster,
                     weighted=TRUE)
        }
        else {
        temp <- residuals.coxph(fit2, type='dfbeta', weighted=TRUE)
        fit2$linear.predictors <- 0*fit$linear.predictors
        temp0 <- residuals.coxph(fit2, type='score', weighted=TRUE)
            }
        fit$var <- t(temp) %*% temp
        u <- apply(as.matrix(temp0), 2, sum)
        fit$rscore <- coxph.wtest(t(temp0)%*%temp0, u, control$toler.chol)$test
        }
    #Wald test
    if (length(fit$coefficients) && is.null(fit$wald.test)) {
        #not for intercept only models, or if test is already done
        nabeta <- !is.na(fit$coefficients)
        if (is.null(init)) temp <- fit$coefficients[nabeta]
        else temp <- (fit$coefficients - init)[nabeta]
        fit$wald.test <-  coxph.wtest(fit$var[nabeta,nabeta], temp,
                      control$toler.chol)$test
        }
    na.action <- attr(m, "na.action")
    if (length(na.action)) fit$na.action <- na.action
    if (model) fit$model <- m
        ## else { ## we might want model=T and X=T
        if (x)  {
            fit$x <- X
            if (length(strats)) fit$strata <- strata.keep
        }
        if (y)     fit$y <- Y
        ##    }
    }
    if (!is.null(weights) && any(weights!=1)) fit$weights <- weights
    
    ##fit$formula <- as.vector(attr(Terms, "formula"))
    fit$formula<-formula(Terms) ## get the environments right
    fit$call <- call
    fit$method <- method
    fit
    }