R/likelihood.R

In multiRec: Analysis of Multi-Type Recurrent Events

#' Calculate the likelihood)
#'
#' @param modelObj the model object constructed by mtre
#' @param params the parameter vector
#' @param type character,if "likelihood" = returns log-likelihood , "likelihood.i" = returns 
#' log-likelihood components , if "deviance" = returns deviance , if "deviance.i" = returns deviance 
#' components , if "resid" = returns score residuals, if 'robust' = returns components for the calculation of robust standard errors
#' @importFrom Rfast rowsums
#' @keywords internal
likelihood = function(params, modelObj, 
                      type=c('likelihood','likelihood.i','deviance','deviance.i', 'resid','robust')) {
 
   type = match.arg(type)

  # Calculate the hazards for all possible events. Returns a matrix with
  # columns for each event type, and rows for each participant and time
  hazards = lapply(modelObj$hazards, function(hazard) {
    # Each component is a consisting of an intrinsic pieces and 0 or more
    # prior pieces.
    pieces = sapply(hazard, function(component) {
      eta   = component$X %*% params[component$beta.ix]
      evName = component$event

      # Multiply by f(counts), except if it's an intrinsic component, in which
      # case the countName will be NULL
      countName = component$countName
      if (!is.null(countName)) {
        N = modelObj$data[[countName]]

        if (!is.null(component$alpha.ix)) alpha = params[component$alpha.ix]
        else alpha=component$alpha

        eta = eta * (N)^alpha
      }

      return(eta)
    })

    hazard = modelObj$link$linkinv(rowSums(pieces),params[modelObj$link$param.ix])
    return(hazard)
  })
  
  hazards = do.call(cbind, hazards)

  # The cumulative hazard of no event (assuming mutually exclusive events), i.e.
  # log(S)
  noEvHaz  = -rowSums(hazards)*modelObj$w


  # Use the selector to pick out the hazard of the event that actually happened
  # at the end of the interval (or 1 if no event)
  hazard = hazards[modelObj$selector]
  hazard[is.na(hazard)] = 1
  
  lik.i = log(hazard) + noEvHaz
  dev.i = -2*lik.i
  
  if (type=='likelihood') return(sum(lik.i))
  else if (type=='likelihood.i') return(lik.i)
  else if (type=='deviance.i') return(dev.i)
  else if (type=='deviance') return(sum(dev.i))
  else if (type=='robust') stop("NOT YET IMPLEMENTED")
  else if (type=='resid') {

    # Calculate analogs of martingale residuals
    cumHazards = modelObj$w*hazards
    d          = matrix(0, nrow=nrow(cumHazards), ncol=ncol(cumHazards))
    d[modelObj$selector] = 1
    residuals  = d-cumHazards

    return(list(residuals=residuals,
                dev.i=dev.i))
  }
}