R/coxph.r

In mgcv: Mixed GAM Computation Vehicle with Automatic Smoothness Estimation

## (c) Simon N. Wood (2013, 2014) coxph model general family. 
## Released under GPL2 ...

cox.ph <- function (link = "identity") { 
## Extended family object for Cox PH.
  linktemp <- substitute(link)
  if (!is.character(linktemp)) linktemp <- deparse(linktemp)
  if (linktemp %in% c("identity")) stats <- make.link(linktemp)
  else if (is.character(link)) {
    stats <- make.link(link)
    linktemp <- link
  } else {
    if (inherits(link, "link-glm")) {
       stats <- link
            if (!is.null(stats$name))
                linktemp <- stats$name
        }
        else stop(linktemp, " link not available for coxph family; available link is \"identity\" ")
  }
  env <- new.env(parent = .GlobalEnv)
  validmu <- function(mu) all(is.finite(mu))

    

    ## initialization is tough here... need data frame in reverse time order,
    ## and intercept removed from X...
  

    preinitialize <- function(G) {
      ## G is a gam pre-fit object. Pre-initialize can manipulate some of its
      ## elements, returning a named list of the modified ones.
      ## sort y (time) into decending order, and
      ## re-order weights and rows of X accordingly
      ## matrix y has strat as second column
      G$family$data <- list()
      y.order <- if (is.matrix(G$y)) order(G$y[,2],G$y[,1],decreasing=TRUE) else
                 order(G$y,decreasing=TRUE)
      G$family$data$y.order <- y.order
      G$y <- if (is.matrix(G$y)) G$y[y.order,] else G$y[y.order]
      attrX <- attributes(G$X)
      G$X <- G$X[y.order,,drop=FALSE]
      attributes(G$X) <- attrX
      G$w <- G$w[y.order]
      G$offset <- G$offset[y.order]
      list(family=G$family,y=G$y,X=G$X,w=G$w,offset=G$offset)
    } ##  preinitialize
    
    postproc <- expression({
    ## code to evaluate in estimate.gam, to do with data ordering and 
    ## baseline hazard estimation...
      ## first get the estimated hazard and prediction information...
      G$X <- Sl.initial.repara(G$Sl,G$X,inverse=TRUE,cov=FALSE,both.sides=FALSE)
      object$family$data <- G$family$hazard(G$y,G$X,object$coefficients,G$w,G$offset)
      rumblefish <- G$family$hazard(G$y,matrix(0,nrow(G$X),0),object$coefficients,G$w)
      s0.base <- exp(-rumblefish$h[rumblefish$r]) ## no model baseline survival 
      s0.base[s0.base >= 1] <- 1 - 2*.Machine$double.eps ## avoid NA later
      ## now put the survivor function in object$fitted
      object$fitted.values <- exp(-object$family$data$h[object$family$data$r]*exp(object$linear.predictors))
      ## compute the null deviance...
      s.base <- exp(-object$family$data$h[object$family$data$r]) ## baseline survival
      s.base[s.base >= 1] <- 1 - 2*.Machine$double.eps ## avoid NA later
      object$null.deviance <- ## sum of squares of null deviance residuals
      2*sum(abs((object$prior.weights + log(s0.base) + object$prior.weights*(log(-log(s0.base)))))) 
      ## and undo the re-ordering...
      y.order <- G$family$data$y.order
      object$linear.predictors[y.order] <- object$linear.predictors
      object$fitted.values[y.order] <- object$fitted.values
      if (is.matrix(object$y)) object$y[y.order,] <- object$y else object$y[y.order] <- object$y  
      object$prior.weights[y.order] <- object$prior.weights
      object$family$data$Rs[y.order,] <- object$family$data$Rs
      object$family$data$Rc[y.order,] <- object$family$data$Rc
    })
    
    initialize <- expression({
        n <- rep(1, nobs)
        if (is.null(start)) start <- rep(0,ncol(x))
    })

    hazard <- function(y, X,beta,wt,offset=0) {
    ## get the baseline hazard function information, given times in descending order in y
    ## model matrix (same ordering) in X, coefs in beta and censoring in wt (1 = death, 0
    ## = censoring)
    
      if (is.matrix(y)) { 
        ## first column is time, second is *numeric* code indicating strata
        strat <- y[,2] ## stratification variable
        y <- y[,1] ## event times 
	strat.lev <- unique(strat) 
	ns <- length(strat.lev) ## number of strata
	nt <- 0;for (i in 1:ns) nt <- nt + length(unique(y[strat==strat.lev[i]]))
        tr.strat <- tr <- rep(0,nt)
        k <- 1
        for (i in 1:ns) {
	  tr0 <- unique(y[strat==strat.lev[i]])
          ind <- k:(k+length(tr0)-1);k <- k + length(tr0)
          tr[ind] <- tr0 ## unique times at this stratification level
	  tr.strat[ind] <- strat.lev[i] ## stratication index for tr,h,q,km and a
        }
      }	else {
        ns <- 1
	tr <- unique(y)
	nt <- length(tr)
      }
      p <- ncol(X)
      eta <- as.double(X%*%beta) + offset
      gamma <- exp(eta)
      if (ns==1) {
        r <- match(y,tr)
        oo <- .C("coxpp",eta,A=as.double(X),as.integer(r),d=as.integer(wt),
               h=as.double(rep(0,nt)),q=as.double(rep(0,nt)),km=as.double(rep(0,nt)),
               n=as.integer(nrow(X)),p=as.integer(ncol(X)),
               nt=as.integer(nt),PACKAGE="mgcv")
	## Now compute the Shoenfeld residuals, set to zero for censored...
        if (p) {
          X <- (X - apply(gamma*X,2,cumsum)/cumsum(gamma))*wt ## Schoenfeld
	  n <- nrow(X)
	  R <- apply(X[n:1,,drop=FALSE],2,cumsum)[n:1,,drop=FALSE] ## score residuals
	  ## now remove the penalization induced drift...
	  x <- (1:n)[wt!=0]; dx <- seq(1,0,length=sum(wt!=0))
	  drift <- approx(x,dx,xout=1:n,method="constant",rule=2)$y
	  R <- R - t(t(matrix(drift,n,ncol(R)))*R[1,])
	  X[wt==0,] <- NA
	} else R <- X  
	return(list(tr=tr,h=oo$h,q=oo$q,a=matrix(oo$A[1:(p*nt)],p,nt),nt=nt,r=r,km=oo$km,Rs = X,Rc = R))
      } else {
        R <- X 
        r <- y*0;a <- matrix(0,p,nt)
	h <- q <- km <- rep(0,nt)
        for (i in 1:ns) { ## loop over strata
          ind <- which(strat==strat.lev[i])
	  trind <- which(tr.strat==strat.lev[i])
	  r0 <- match(y[ind],tr[trind])
	  nti <- length(trind)
	  etai <- if (p>0) eta[ind] else eta
          oo <- .C("coxpp",etai,A=as.double(X[ind,]),as.integer(r0),d=as.integer(wt[ind]),
               h=as.double(rep(0,nti)),q=as.double(rep(0,nti)),km=as.double(rep(0,nti)),
               n=as.integer(length(ind)),p=as.integer(p),
               nt=as.integer(nti),PACKAGE="mgcv")
	  ## now paste all this into return fields
	  h[trind] <- oo$h
	  q[trind] <- oo$q
	  km[trind] <- oo$km
	  r[ind] <- r0 ## note that indexing is to subsetted tr
	  a[,trind] <- matrix(oo$A[1:(p*nti)],p,nti)
	  ## compute Schoenfeld resiudals, within stratum
          if (p) {
            Xs <- X[ind,] <- (X[ind,,drop=FALSE] - apply(gamma[ind]*X[ind,,drop=FALSE],2,cumsum)/cumsum(gamma[ind]))*wt[ind]
	    n <- nrow(Xs)
	    Rs <- apply(Xs[n:1,,drop=FALSE],2,cumsum)[n:1,,drop=FALSE] ## score residuals
	    ## now remove the penalization induced drift...
	    x <- (1:n)[wt[ind]!=0];
	    if (length(x)) {
	      dx <- seq(1,0,length=sum(wt[ind]!=0))
	      drift <- approx(x,dx,xout=1:n,method="constant",rule=2)$y
	      Rs <- Rs - t(t(matrix(drift,n,ncol(R)))*Rs[1,])
	    }  
	    R[ind,] <- Rs
	  }  
        }
	X[wt==0,] <- NA
	return(list(tr=tr,h=h,q=q,a=a,nt=nt,r=r,km=km,strat=strat,tr.strat=tr.strat,Rs=X,Rc=R))
      }
    } ## hazard

    residuals <- function(object,type=c("deviance","martingale","score","schoenfeld")) {
      type <- match.arg(type)
      w <- object$prior.weights;log.s <- log(object$fitted.values)
      res <- w + log.s ## martingale residuals
      if (type=="deviance") { 
        log.s[log.s>-1e-50] <- -1e-50
        res <- sign(res)*sqrt(-2*(res + w * log(-log.s)))
      } else if (type=="score") {
        res <- object$family$data$Rc
	term <- list()
        if (object$nsdf) {
	  ii <- 1:object$nsdf; sd <- sqrt(diag(object$Vp)[ii])
	  res[,ii] <- t(sd*t(res[,ii]))
	  term$parametric <- ii
        }
	m <- length(object$smooth)
	if (m) for (i in 1:m) {
          ii <- object$smooth[[i]]$first.para:object$smooth[[i]]$last.para
	  R <- chol(object$Vp[ii,ii]);
	  ## so R^T R = V_b for this term. If b' = R^{-T} b, V_b' = I
	  ## and Xb = XR^T b'...
	  res[,ii] <- res[,ii] %*% t(R)
	  term[[object$smooth[[i]]$label]] <- ii
        }
        attr(res,"term") <- term 
	colnames(res) <- names(coef(object))
      } else if (type=="schoenfeld") {
        res <- object$family$data$Rs
	colnames(res) <- names(coef(object))
      }
      res 
    } ## residuals


    predict <- function(family,se=FALSE,eta=NULL,y=NULL,
               X=NULL,beta=NULL,off=NULL,Vb=NULL) {
      ## prediction function.
      if (is.matrix(y)) { 
	strat <- y[,2];y <- y[,1]
	if (is.null(family$data$strat)) stop("something wrong with stratified prediction")
	ii <- order(strat,y,decreasing=TRUE) ## C code expects non-increasing
	strat <- strat[ii]
      } else {
        ii <- order(y,decreasing=TRUE) ## C code expects non-increasing
        strat <- NULL
      }
      if (sum(is.na(y))>0) stop("NA times supplied for cox.ph prediction")
      X <- X[ii,,drop=FALSE];y <- y[ii];
      n <- nrow(X)
      if (is.null(off)) off <- rep(0,n)
      if (is.null(strat)) {
        oo <- .C("coxpred",as.double(X),t=as.double(y),as.double(beta),as.double(off),as.double(Vb),
                a=as.double(family$data$a),h=as.double(family$data$h),q=as.double(family$data$q),
                tr = as.double(family$data$tr),
                n=as.integer(n),p=as.integer(ncol(X)),nt = as.integer(family$data$nt),
                s=as.double(rep(0,n)),se=as.double(rep(0,n)),PACKAGE="mgcv")
        s <- sef <- oo$s
        s[ii] <- oo$s
        sef[ii] <- oo$se 
      } else { ## stratified fit, need to unravel everything by strata
        pstrata <- unique(strat)
	ns <- length(pstrata)
	p <- ncol(X)
	s <- sef <- rep(0,length(y)) 
	for (i in 1:ns) {
          ind <- which(strat==pstrata[i]) ## prediction data index
	  trind <- which(family$data$tr.strat == pstrata[i])
	  n <- length(ind)
	  oo <- .C("coxpred",as.double(X[ind,]),t=as.double(y[ind]),as.double(beta),as.double(off),as.double(Vb),
                a=as.double(family$data$a[,trind]),h=as.double(family$data$h[trind]),q=as.double(family$data$q[trind]),
                tr = as.double(family$data$tr[trind]),
                n=as.integer(n),p=as.integer(p),nt = as.integer(length(trind)),
                s=as.double(rep(0,n)),se=as.double(rep(0,n)),PACKAGE="mgcv")
	  s[ind] <- oo$s
	  sef[ind] <- oo$se
        } ## strata loop
	s[ii] <- s
        sef[ii] <- sef   
      }    
      if (se) return(list(fit=s,se.fit=sef)) else return(list(fit=s))
    } ## predict

    rd <- qf <- NULL ## these functions currently undefined for Cox PH

    ll <- function(y,X,coef,wt,family,offset=NULL,deriv=0,d1b=0,d2b=0,Hp=NULL,rank=0,fh=NULL,D=NULL) {
    ## function defining the Cox model log lik.
    ## Calls C code "coxlpl"
    ## deriv codes: 0   - evaluate the log likelihood
    ##              1   - evaluate the grad and Hessian, H, of log lik w.r.t. coefs (beta)
    ##              2   - evaluate tr(Hp^{-1}dH/drho) - Hp^{-1} must be supplied in fh
    ##                    and db/drho in d1b. Not clear if a more efficient algorithm exists for this.
    ##              3   - evaluate d1H =dH/drho given db/drho in d1b 
    ##                    (2 is for evaluation of diagonal only)
    ##              4 -  given d1b and d2b evaluate trHid2H= tr(Hp^{-1}d2H/drhodrho')
    ## Hp is the preconditioned penalized Hessian of the log lik
    ##    which is of rank 'rank'.
    ## fh is a factorization of Hp - either its eigen decomposition 
    ##    or its Choleski factor
    ## D is the diagonal pre-conditioning matrix used to obtain Hp
    ##   if Hr is the raw Hp then Hp = D*t(D*Hr)
      
      if (is.matrix(y)) {
        ## first column is time, second is *numeric* code indicating strata
        strat <- y[,2] ## stratification variable
        y <- y[,1] ## event times 
	strat.lev <- unique(strat) 
	ns <- length(strat.lev) ## number of strata
      }	else ns <- 1
     
      p <- ncol(X)
      deriv <- deriv - 1
      mu <- X%*%coef + offset
      g <- rep(0,p);H <- rep(0,p*p)
      if (deriv > 0) {
        M <- ncol(d1b)
        d1H <- #if (deriv==1) rep(0,p*M) else
	       rep(0,p*p*M)
      } else M <- d1Ho <- d1H <- 0
      if (deriv > 2) {
        d2H <- rep(0,p*M*(M+1)/2)
        if (is.list(fh)) {
          ev <- fh
        } else  { ## need to compute eigen here
          ev <- eigen(Hp,symmetric=TRUE)
          if (rank < p) ev$values[(rank+1):p] <- 0
        } 
        X <- X%*%(ev$vectors*D)
        d1b <- t(ev$vectors)%*%(d1b/D); d2b <- t(ev$vectors)%*%(d2b/D)
      } else d2Ho <- trHid2H <- d2H <- 0

      for (j in 1:ns) { ## loop over strata
        ind <- if (ns==1) 1:length(y) else which(strat==strat.lev[j]) ## index for points in this strata
        tr <- unique(y[ind])
        r <- match(y[ind],tr)
        ## note that the following call can not use .C(C_coxlpl,...) since the ll
        ## function is not in the mgcv namespace.
	cderiv <- if (deriv==1) 2 else deriv ## reset 1 to 2 to get whole d1H returned
        oo <- .C("coxlpl",as.double(mu[ind]),as.double(X[ind,]),as.integer(r),as.integer(wt[ind]),
            as.double(tr),n=as.integer(length(y[ind])),p=as.integer(p),nt=as.integer(length(tr)),
            lp=as.double(0),g=as.double(g),H=as.double(H),
            d1b=as.double(d1b),d1H=as.double(d1H),d2b=as.double(d2b),d2H=as.double(d2H),
            n.sp=as.integer(M),deriv=as.integer(cderiv),PACKAGE="mgcv");
	if (j==1) {
          lp <- oo$lp
	  lb <- oo$g
	  lbb <- matrix(oo$H,p,p)
        } else { ## accumulating over strata
          lp <- oo$lp + lp
	  lb <- oo$g + lb
	  lbb <- matrix(oo$H,p,p) + lbb
	}
        if (deriv==1) {
	  #d1Ho <- matrix(oo$d1H,p,M) + if (j==1) 0 else d1Ho
	  ind <- 1:(p^2)
          if (j==1) d1Ho <- rep(0,M)
          for (i in 1:M) { 
            d1Ho[i] <- d1Ho[i] + sum(fh*matrix(oo$d1H[ind],p,p)) ## tr(Hp^{-1}dH/drho_i)
            ind <- ind + p^2
          }
	} else if (deriv>1) {
          ind <- 1:(p^2)
          if (j==1) d1Ho <- list()
          for (i in 1:M) { 
            d1Ho[[i]] <- if (j==1) matrix(oo$d1H[ind],p,p) else matrix(oo$d1H[ind],p,p) + d1Ho[[i]]
            ind <- ind + p^2
          }
        } 
        if (deriv>2) { 
          d2Ho <- if (j==1) matrix(oo$d2H,p,M*(M+1)/2) else d2Ho + matrix(oo$d2H,p,M*(M+1)/2)
        }
      } ## strata loop
      if (deriv>2) {
        d <- ev$values
        d[d>0] <- 1/d[d>0];d[d<=0] <- 0
        trHid2H <- colSums(d2Ho*d)
      }
      assign(".log.partial.likelihood", oo$lp, envir=environment(sys.function()))
      list(l=lp,lb=lb,lbb=lbb,d1H=d1Ho,d2H=d2Ho,trHid2H=trHid2H)
    } ## ll

    # environment(dev.resids) <- environment(aic) <- environment(getTheta) <- 
    # environment(rd)<- environment(qf)<- environment(variance) <- environment(putTheta) 
    #environment(aic) <- 
    environment(ll) <- env
    structure(list(family = "Cox PH", link = linktemp, linkfun = stats$linkfun,
        linkinv = stats$linkinv, ll=ll, ## aic = aic, 
        mu.eta = stats$mu.eta, 
        initialize = initialize,preinitialize=preinitialize,postproc=postproc,
        hazard=hazard,predict=predict,residuals=residuals,
        validmu = validmu, valideta = stats$valideta, 
        rd=rd,qf=qf,drop.intercept = TRUE,
        ls=1, ## signal ls not needed
        available.derivs = 2 ## can use full Newton here
        ),
        class = c("general.family","extended.family","family"))
} ## cox.ph