R/ah.fit.R

In katehu/addhazard: Fit Additive Hazards Models for Survival Analysis

#######################################################
# Author: Kate HU
# Date: July 27th, 2016
# Task: write function when assuming there may be ties
# Date: March 20, 2017
# Task: change censor variable name to "event"
# To do: interval censoring
#######################################################
#' @import ahaz

ah.fit<-function(Y,Z,wts,robust){
  
  n.obs <- dim(Y)[1]
  if(!is.vector(Z)) {
    if(!is.matrix(Z)) warning ("Covariates Z has to be either a vector or a matrix")
    n.par <- dim(Z)[2]
  } else {
    n.par <- 1
  }
  
  t.fail <- Y[,1]
  event <- Y[,2] # 1 means death and 0 means censored
  if (missing(wts) || is.null(wts)) wts <- rep(1, n.obs)
  
  scale(Z, center = rep(1,n.par),scale=F)
  # some statistics about time points----------------------------------------
  # ordering time based on 1st: failure time and 2nd: events
  new.order <- order(Y[,1], Y[,2])
  Y.ord <- Y[new.order, ]
  # Y.ord[,2] is the ordered events
  # Y.ord[,1] is the ordered time
  if (is.vector(Z)){
    Z.ord <- Z[new.order]
    wts.ord <- wts[new.order]
  }
  #add matrix later
  
  # count how many events at each sorted time
  event.table <- tapply(Y.ord[,2], Y.ord[,1], sum)
  # sorted unique time, either censoring or failure time
  t.marks <- as.numeric(names(event.table))
  # number of time points
  n.t.marks <- length(t.marks)
  
  # number of events at each time point
  n.events.at.t.marks <- as.numeric(event.table)
  
  # count total Z for people who are leaving the risk set at each time point
  Z.sum.lost.at.t.marks <- tapply(Z.ord * wts.ord, Y.ord[,1], sum)
  # count total Z for people who have events at each time point
  Z.ord.times.event <- Z.ord * wts.ord * Y.ord[,2]
  Z.sum.at.event.at.t.marks <- tapply(Z.ord.times.event, Y.ord[,1], sum)
  
  # some statistics about time segements --------------------------------------
  # number of people at risk for each segment
  # the first time segment includes all the individuals
  n.in.dur <- c(n.obs, n.obs - tapply(Y.ord[,1] * wts.ord, Y.ord[,1], length))
  # count cumulative total Z lost for each duration
  Z.cumsum.lost.in.dur <- cumsum(c(0, Z.sum.lost.at.t.marks))
  # count cumulative total Z for each duration
  Z.cumsum.in.dur <- sum(Z) - Z.cumsum.lost.in.dur
  
  # divide the timeline to small intervals
  t.diff <- t.marks- c(0,t.marks[-n.t.marks])
  
  # Section 1: calculate eta(s) ------------------------------------------------
  # when Z is a vector
  # n.in.dur for the last duration is 0
  eta <- Z.cumsum.in.dur[-(n.t.marks+1)]/n.in.dur[-(n.t.marks+1)]
  eta.cum <- cumsum(eta * t.diff)
  
  # Section 2: calculate h(s) ------------------------------------------------
  
  h <- n.events.at.t.marks/n.in.dur[-(n.t.marks+1)]
  h.cum <- cumsum(h)
  # Lambda starts from t1 and ends right before T
  Lambda <- h.cum[-n.t.marks] - eta.cum[-1] * theta
  
  
  
  # time when  fail events occur
  tmp <- t.fail * event
  t1.marks <- sort(tmp[tmp!=0])
  
  # mark each person's event time on the timeline
  sub.pos.on.t.marks <- match(t.fail, t.marks)
  
  # mark all the event times on the timeline
  event.pos.on.t.marks <- match(t1.marks, t.marks)
  
  
  eta.num <- eta.den <- matrix(NA, nrow=n.par, ncol = n.t.marks)
  z.death <- matrix(NA, nrow=n.par, ncol = n.t.marks)
  
  n.death<-NULL
  eta.den.vec<-NULL
  
  # Centralized the covariates
  # Z <- scale( Z,center= rep(1,n.par),scale=F)
  
  # Divide the timeline to small intervals according to  events or censoring
  
  t.sorted <- sort(t.fail)
  t.unique <- unique(t.sorted)
  #  The length of eta
  eta.ncol <- length(t.unique)
  t.diff <- t.unique - c(0,t.unique[-eta.ncol])
  
  # Section 1: calculate eta(t) ------------------------------------------------
  # eta(t):  average Z for patients still at risk at time point t
  # eta(t) only changes when a subject fails or censored, i.e., t.unique
  
  
  ### We calculate the numerator and the denomiators of eta separtedly
  ### We create the spaces for the numerator and the denomiators of eta
  ### Z_i is a n.par x 1 vector , therefore
  
  
  ### We calcuate the numerator and the denomiators of eta
  ### At each unique failure time, we calcuate for each covariate(age,sex,...)
  ### Then we calculate for all unique failure times
  z.death <- eta.num <- eta.den <- matrix(data=NA, nrow=n.par, ncol = eta.ncol)
  n.death<-NULL
  eta.den.vec<-NULL
  for ( i in seq_along(eta.ncol)){
    idx1 <- t.fail >= t.unique[i]
    idx2 <- which(t.fail == t.unique[i])
    
    n.death[i] <- sum((event*wts)[idx2])
    eta.den.vec[i]<-sum(wts[idx1])
    
    for (j in 1:n.par){
      eta.num[j,i] <-sum((Z*wts)[idx1,j])
      
      ########################## IF we calculate Lambda ###############
      ### calculate the cumulative eta
      ### The sum of Z of the subjests who have died at the jth dicrete time
      
      z.death[j,i]<-sum((Z[,j]*event*wts)[idx2])
    }
    eta.den[,i]<-eta.den.vec[i]
  }
  
  
  
  ### Calculate eta, which is a  n.par X eta.col matrix
  dLambda1<-n.death/eta.den.vec
  #print(dLambda1)
  eta<-eta.num/eta.den
  ###############################################################################
  
  ###############################################################################
  ###Calculate theta. It doens't dependent on t
  
  ### Calculate the numerator of theta
  
  
  
  ###Generate a new eta matrix(vector) corresponding to each observation such that
  ### eta.each.person_i is the average of Z for patients still at risk at
  ###the failure time of  individual i
  
  ### create the space fo this  n.obs x n.par matrix (vector)
  
  #### Individuals
  
  
  ### find which
  match.eta<-match(t, t.unique)
  eta.i<-matrix(as.vector(eta[,match.eta]),ncol=n.par,byrow=TRUE)
  
  
  ###Calculate theta.numerator
  
  
  # z.cen.1<-(Z-eta.i)*wts*event
  A=fit$D
  
  
  
  
  #theta.denominator<-Reduce("+",theta.denominator.timepoint)
  
  ###Calculate theta.denominator
  
  
  # A<-matrix(0,n.par,n.par)
  
  #for ( i in 1:eta.ncol){
  
  #  center<-Z[failure.time>= time.interval[i]]
  # eta.each.person[failure.time>= time.interval[i]]
  
  #   idx1<-t>=  t.unique[i]
  #  rep.times<-sum(idx1)
  # eta.stack<-matrix(rep(as.vector(eta[,i]),rep.times), byrow=TRUE,ncol=n.par)
  # temp<-(Z[idx1,]-eta.stack)*sqrt(wts[idx1])
  
  #A<-
  #A+t(temp)%*%temp *t.diff[i]
  
  #}
  
  
  
  
  
  # theta.denominator<-Reduce("+",theta.denominator.timepoint)
  
  ### Calculate theta
  iA<-solve(A)
  theta.num<-fit$d
  theta<-iA%*%theta.num
  
  
  
  ###Calculate the variance
  ###B part
  
  
  B<-fit$B
  var<-iA%*% B %*% iA
  
  
  eta.cum<-as.matrix(apply(t(eta)*t.diff,2,cumsum))
  
  
  
  #for (k in 1: n.obs){
  #       int.upper<-match(t[k],t.unique)
  #     for (l in 1: int.upper){
  #        temp<-(Z[k,]-eta[,l])
  #       A.i[,,k] <-A.i[,,k]+temp%*%t(temp)*t.diff[l]*wts[k]
  #      z.cen.2[k,]<-z.cen.2[k,]+temp*dLambda1[l]*wts[k]
  
  # }
  
  # cumhaz<-function(fit){
  #    eta.cum<-matrix(data=NA,nrow=fit$n.par,ncol=fit$eta.ncol)
  #     Lambda1<-rep(0,fit$eta.ncol)
  #     eta.cum[,1]<-fit$eta[,1]*fit$t.diff[1]
  #     Lambda1[1]<-fit$dLambda1[1]
  #     for ( i in 2:fit$eta.ncol){
  #         eta.cum[,i]<-eta.cum[,i-1]+fit$eta[,i]*fit$t.diff[i]
  #         Lambda1[i]<-Lambda1[i-1]+ fit$dLambda1[i]
  #     }
  #     print(Lambda1)
  #     return(Lambda1-t(eta.cum)%*%fit$coef)
  #
  # }
  
  effects<-as.vector(Z%*%theta)
  
  fit<-list(coef=theta,
            var=var,
            A=A,
            iA=iA,
            eta=eta,
            eta.cum=eta.cum,
            match.eta=match.eta,
            eta.i=eta.i,
            eta.den.vec=eta.den.vec,
            dLambda1=dLambda1,
            n.par=n.par,
            n.obs=n.obs,
            z.death=z.death,
            n.death=n.death,
            eta.ncol=eta.ncol,
            event=event,
            time=t,
            t.diff=t.diff,
            t.unique=t.unique,
            effects=effects,
            robust=robust)
  
  
  lambda0<-dLambda1-as.vector(t(eta)%*%theta*t.diff)
  Lambda0<-cumsum(lambda0)
  
  
  
  M<-(event-Lambda0[match.eta]-effects*t)
  
  
  
  z.resid<-Z*M
  eta.resid2<-apply(t(eta)*lambda0,2,cumsum)
  
  
  eta.resid3<-eta.cum[match.eta,]*effects
  
  
  eta.resid<-eta.i*event-eta.resid2[match.eta,]-eta.resid3
  
  ###nobs x n.par
  z.cen.resid<-(z.resid-eta.resid)
  
  #A.i<-array(rep(0,n.par*n.par*n.obs),dim=c(n.par,n.par,n.obs))
  #for (k in 1: n.obs){
  #     int.upper<-match(t[k],t.unique)
  #   for (l in 1: int.upper){
  #      temp<-(Z[k,]-eta[,l])
  #     A.i[,,k] <-A.i[,,k]+temp%*%t(temp)*t.diff[l]*wts[k]
  #    z.cen.2[k,]<-z.cen.2[k,]+temp*dLambda1[l]*wts[k]
  
  #}
  #}
  
  
  fit$Lambda0<-Lambda0
  fit$lambda0<-lambda0
  
  fit$M<-M
  fit$resid<-z.cen.resid
  
  #z.cen.3<-as.matrix(apply(A.i,3,function(s){s%*%fit$coef}))
  #print(A.i[,,1]%*%theta)
  #print(A.i[,,2]%*%theta)
  #print(dim(z.cen.3))
  #print(dim(z.cen.2))
  #fit$z.cen<-(z.cen.1-z.cen.2-z.cen.3)
  
  #print(z.cen.2)
  return(fit)
}


residuals.ah<-function (object, type = c("martingale", "pseudoscore",
                                         "pseudodfbeta","martingale.increment")){
  
  type <- match.arg(type)
  if (type=="martingale") return(object$M)
  if (type=="pseudoscore") return(object$resid)
  if (type=="pseudodfbeta") return(object$iA%*%object$resid)
  #if (type=="martingale.increment"){
  # dM<-matrix(0,nrow=object$n.obs,ncol=object$eta.ncol)
  #for  ( i in 1: object$eta.ncol){
  #    dM[object$match.eta<=i,] <- (object$event-object$lambda0[i]
  #                      -object$effects*object$t.diff[i])[object$match.eta<=i]
  #  }
  #return(dM)
  #}
}