R/fitObject.R

In multiRec: Analysis of Multi-Type Recurrent Events

#' Construct a object returned by multiRec
#'
#' @param fit the output from optim
#' @param modelObj the model object, as created by multiRec
#' @param robust logical, TRUE if robust standard errors are required
#' @param init numeric, the vector of initial parameter estimates
#' @param fitDetails logical, if TRUE the returned fit object contains
#'   additional information used to create fit diagnostics.
#' @param hessian character, use the Hessian calculated in optim or calculated
#'   separately in numDeriv::hessian().
#' @param omit omit, information about dropped observations
#'
#' @return an object of class multiRec
#' @importFrom stats pnorm
#' @importFrom MASS ginv
#' @keywords internal
makeFitObject = function(fit, modelObj, robust=TRUE, init=init, fitDetails=FALSE,
                         hessian=c('optim','numDeriv'), omit=NULL) {
  param = fit$par
  log.lik=fit$value
  deviance = -2*log.lik

  hessian = match.arg(hessian)

  obj = list()

  # make names for the rows
  shortNames = sprintf('%s[%s] %s',
                       pkgEnv$names[,1],
                       gsub('(Intrinsic)','*',pkgEnv$names[,2], fixed=TRUE),
                       gsub('(Intercept)','(Int)',pkgEnv$names[,3], fixed=TRUE)
  )

  # Get the hessian
  if(hessian == 'optim') H = fit$hessian
  else H  = numDeriv::hessian(likelihood, param, modelObj=modelObj)

  J.int  = numDeriv::jacobian(likelihood, param, modelObj=modelObj, type='likelihood.i')
  J = matrix(NA,nrow = length(unique(modelObj$id)),ncol = ncol(J.int))
  for (j in 1:ncol(J)){
    temp = ave(J.int[,j],modelObj$id,FUN=sum)
    temp = temp[!duplicated(modelObj$id)]
    J[,j]=temp
  }
  naive.var = -ginv(H)

  if (robust) {
    # Get the Jacobian and calculate the sandwich estimate of the variance
    # covariance matrix
    n     = nrow(J)
    bread = naive.var # H
    meat  = crossprod(J) # crossprod(J)/n
    obj$var = (bread %*% meat %*% bread)
    obj$naive.var = naive.var
    rownames(obj$naive.var) = colnames(obj$naive.var) = shortNames
  }
  else obj$var=naive.var
  rownames(obj$var) = colnames(obj$var) = shortNames

  # Get the diagonal elements. It's possible to get negative values,
  # particularly with the identity link.
  var = diag(obj$var)
  if (any(var<0)) {
    fit$convergence = 99
    fit$message = 'Negative variance estimates'
    var[var<0] = NA
  }
  se = sqrt(var)
  coef = data.frame(estimate=param,
                    se=se,
                    z=param/se,
                    'P(>|Z|)'=round(2*pnorm(abs(param/se),lower.tail=FALSE),4),
                    check.names=FALSE)
  obj$coefficients = cbind(pkgEnv$names, coef)
  rownames(obj$coefficients) = shortNames

  # Calculate additional quantities
  obj$calls    = sapply(modelObj$hazards,attr,'call')
  obj$loglik   = log.lik
  obj$deviance = deviance
  obj$data     = modelObj$data
  obj$n        = nrow(modelObj$data)
  obj$nEvent   = table(modelObj$data[[modelObj$eventVar]])
  obj$first    = J
  obj$convergence = fit$convergence
  obj$iterations = fit$counts['function']
  obj$link = modelObj$link[c('link', 'linkfun','linkinv')]
  obj$na.action = omit

  # Calculate something like martingale residuals
  residuals = likelihood(param, modelObj, type='resid')
  obj$residuals = residuals$residuals
  obj$deviance.components = residuals$dev.i
  obj$score.components = J

  if (fit$convergence!=0) warning('Model did not converge')

  if (fitDetails) {
    obj$fitDetails = list(param=obj$coefficients[,c('hazard', 'component',
                                                    'parameter', 'estimate', 'se')],
                          init=init,
                          modelObj=modelObj)
  }

  class(obj) = 'multiRec'
  return(obj)
}