Predict in sample hazard rates that could be used for limited information maximum likelihood two-stage selection in Cox proportional hazard models.

coxLIML.R

#' Predict in sample hazard rates that could be used for limited information maximum likelihood two-stage selection in Cox proportional hazard models
#' 
#'  @param obj a fitted model object created with \code{coxph}.
#'  @param data the name of the data frame used to create \code{obj}.
#'  @param idvar character string of the variable in \code{data} that identifies individuals.
#'  @param timeVar character string of the variable in \code{data} that identifies the \code{time} variable used in \code{obj}.
#'  
#'   @return Returns the original \code{data} data frame with an additional variable added at the end called \code{PredictedHRate}. This contains the hazard rate predicted from the model \code{obj} for each observation. 
#'   
#'   @importFrom survival basehaz
#'   @export 




coxLIML <- function(obj, data, idVar, timeVar){
  require(survival)
  if (class(obj) != 'coxph'){
    stop('obj must be a coxph class object.')
  }
  
  # Function to find one row's predicted hazard rate
  OneRow <- function(x, overValue, ncoef = ncoef, Coef = Coef){
    To <- ncoef + (overValue - 1)
    Row1 <- x[, overValue:To]
    PHR <- x$hazard * exp(sum(Row1*Coef))
    PHR
  }
  
  # Find baseline hazards and coefficient estimates
  BH <- basehaz(obj)  
  Coef <- obj$coefficients
  ncoef <- length(Coef)
  NRows <- nrow
  
  # Created data set of only the timevar + model independent variables
  if (!is.null(BH$strata)){
    StrataName <- gsub('=.', '', BH$strata[1])
    CoefTime <- c(idVar, timeVar, StrataName, names(Coef))
    DataSub <- data[, CoefTime]
    BH$strata <- gsub('^.*?=', '', BH$strata)
    names(BH) <- c('hazard', timeVar, StrataName)
    DataSub <- merge(DataSub, BH, by = c(timeVar, StrataName), all = TRUE)
    
    TempVect <- vector()  
    for(i in 1:nrow(DataSub)){
      Temp <- OneRow(x = DataSub[i,], overValue = 4, ncoef = ncoef, Coef = Coef)
      TempVect <- c(TempVect, Temp)
    }
    DataSub$PredictedHRate <- TempVect
    DataSubSmall <- DataSub[, c(idVar, timeVar, 'PredictedHRate')]  
    
  } else if (is.null(BH$strata)){
    CoefTime <- c(idVar, timeVar, names(Coef))
    DataSub <- data[, CoefTime]
    names(BH) <- c('hazard', timeVar)
    DataSub <- merge(DataSub, BH, by = timeVar, all = TRUE)
    
    TempVect <- vector()  
    for(i in 1:nrow(DataSub)){
      Temp <- OneRow(x = DataSub[i,], overValue = 3, ncoef = ncoef, Coef = Coef)
      TempVect <- c(TempVect, Temp)
    }
    DataSub$PredictedHRate <- TempVect
    DataSubSmall <- DataSub[, c(idVar, timeVar, 'PredictedHRate')]  
  }
  Out <- merge(data, DataSubSmall, by = c(idVar, timeVar), all = TRUE)    
}