R/init.R

In multiRec: Analysis of Multi-Type Recurrent Events

#' Get starting values for proportional hazards style likelihood
#'
#' The estimate of the intercept is the mean of the log baseline hazard.
#'
#' @param modelObj a model object built by mtre
#'
#' @importFrom survival coxph basehaz Surv
#' @importFrom stats as.formula coefficients update
#' @keywords internal
init = function(modelObj, shift=FALSE) {
  resetParamVector()

  #linkName = modelObj$link$link

  # Construct a fake dataset by joining X matrices and giving the columns names
  # based on beta.ix.  Trying to match beta values to terms in the design
  # matrices is too messy given R's odd and complicated way of naming terms in a
  # regression output.
  lapply(modelObj$hazards, function(hazard) {
    # Construct the formula
    hazardName = attr(hazard,'hazardName')
    formula = paste0('Surv(',modelObj$timeVar,',',modelObj$time2Var,',',modelObj$eventVar,'=="',hazardName,'")')
    X = lapply(hazard, function(component) {
      # Handle the alpha, if any
      if (!is.null(component$alpha.ix)) {
        setParamVector(component$alpha.ix,1)
        setParamNames(component$alpha.ix,
                      hazard=hazardName,
                      component=component$event,
                      param='(alpha)')
      }

      # construct the X piece
      X = component$X
      refNames    = sprintf('%s~%s~beta~%s',hazardName,component$event,colnames(X))
      setParamNames(component$beta.ix,
                    hazard=hazardName,
                    component=ifelse(component$event=='(Intrinsic)',
                                     component$event,
                                     paste('N of prior',component$event)),
                    param=colnames(X))

      colnames(X) = paste0(modelObj$uniquePrefix, component$beta.ix)

      countName = component$countName
      if (!is.null(countName)) {
        N = modelObj$data[[countName]]
        for (j in 1:ncol(X)) X[,j] = X[,j]* N
      }
      return(X)
    })

    X = as.data.frame(do.call(cbind, X))
    formula = as.formula(paste0(formula,' ~ ', paste(names(X),collapse=' + ')))

    # Add the left hand side variables
    X[[modelObj$timeVar]]  = modelObj$data[[modelObj$timeVar]]
    X[[modelObj$time2Var]] = modelObj$data[[modelObj$time2Var]]
    X[[modelObj$eventVar]] = modelObj$data[[modelObj$eventVar]]

    # fit a Cox PH model. Note that we're using the PH model regardless of the
    # link function. We just need reasonable starting values and this should be
    # good enough
    fit  = suppressWarnings(coxph(formula=formula,data=X))
    beta = coefficients(fit)

    # estimate the intercept (it'll come back with a value of NA)
    aveHaz = max(1e-6,mean(basehaz(fit, centered=FALSE)$hazard))
    beta[is.na(beta)] =log(aveHaz)

    # Now assign values to the parameter vector
    indices = substring(names(beta),nchar(modelObj$uniquePrefix)+1,nchar(names(beta)))
    indices = as.integer(indices)
    setParamVector(indices, beta)
  })

  # If we're using the Yeo-Johnson link, set its parameter to 0 (i.e. start
  # with the Cox model)
  if (!is.null(modelObj$link$param.ix)) {
    setParamVector(modelObj$link$param.ix, 0)
    setParamNames(modelObj$link$param.ix,
                  hazard='(All)',
                  component='(Link)',
                  param='(Power)')
  }

  return(paramVector())
}