R/residuals.R

In EstimationTools: Maximum Likelihood Estimation for Probability Functions from Data Sets

#' @title Extract Residuals from \code{maxlogL} model.
#'
#' @author Jaime Mosquera Gutiérrez, \email{jmosquerag@unal.edu.co}
#'
#' @description
#' `r lifecycle::badge("experimental")`
#'
#' \code{residuals.,axlogL} is the \code{maxlogLreg} specific method for the
#' generic function residuals which extracts the residuals from a fitted model.
#'
#' @aliases residuals.maxlogL
#'
#' @param object an object of \code{\link{maxlogL}} class which summary is desired.
#' @param type a character with the type of residuals to be computed required.
#'             The default value is \code{type = "rqres"}, which is used to
#'             compute the normalized randomized quantile residuals.
#' @param parameter a character which specifies the parameter whose normalized
#'                  quantile residuals will be computed.
#' @param ... for extra arguments.
#'
#' @details
#' For \code{type = "deviance"}, the residuals are computed using the following
#' expression
#'
#' \deqn{r^D_i = \mbox{sign}(y_i - \hat{\mu}_i) d_i^{1/2},}
#'
#' where \eqn{d_i} is the residual deviance of each data point. In this context,
#' \eqn{\hat{\mu}} is the estimated mean, computed as the expected value using
#' the estimated parameters of a fitted \code{maxlogLreg} model.
#'
#' On the other hand, for \code{type = "response"} the computation is simpler
#'
#' \deqn{r^R_i = (y_i - \hat{\mu}_i).}
#'
#' @export
residuals.maxlogL <- function(object, parameter = NULL,
                              type = c("deviance", "response"),
                              ...){
  type <- match.arg(type)
  if ( is.null(parameter) ) parameter <- object$outputs$par_names[1]
  parameter <- match.arg(parameter, choices = object$outputs$par_names)

  y <- object$outputs$response
  support <- object$inputs$support
  dist <- deparse(object$inputs$y_dist[[3]])
  delta <- object$inputs$cens

  if ( is.null(support) & type %in% c("deviance", "pearson",
                                      "response") ){
    stop(paste0(type, " residuals cannot be computed if a support is",
                "not defined. Please, refit your 'maxlogLreg' model",
                "specifying the 'suport' argument."))
  }
  if (type == 'deviance'){
    mean <- moment_integration_maxlogL(object)
    resid <- sign(y - mean)*sqrt(sum(dev.resids(object)))
  }
  if (type == 'response'){
    mean <- moment_integration_maxlogL(object)
    resid <- y - mean
  }
  if (type == 'martingale'){
    if ( is.Surv(object$inputs$y_dist) ){
      cumHaz <- hazard_integration_maxlogL(object)
      delta <-
      resid <- delta[,2] - cumHaz  # Martingale for right censored data
    }
  }
  names(resid) <- 1:length(y)
  return(resid)
}

#==============================================================================
# Deviance for each data point ------------------------------------------------
#==============================================================================
# title Deviance of each data point for \code{maxlogLreg} outputs
# This function computes the deviance for each data point from the
# response variable given a fitted model.
#
# param object an object of \code{maxlogL} class generated by
#       \code{\link{maxlogLreg}} function.
#
# details
# The function requires a fitted model with \code{maxlogLreg}.
#
# return
# This function returns a list with the following elements:
# \enumerate{
#    \item \code{deviance_i}: the contribution of each value of the
#    response variable in the data set to the deviance.
#    \item \code{proposed_deviance_i}: log-likelihood of data given the fitted
#    model in \code{object} argument.
#    \item \code{saturated_deviance_i}. log-likelihood of data given the
#    saturated model.
# }
dev.resids <- function(object){
  if (object$outputs$type != "maxlogLreg")
    stop(paste("'dev.resids()' method only useful for models with ",
               "covariates. \n Use 'maxlogLreg' in order to ",
               "take advantage of this method."))
  # Details of Proposed model
  distr <- object$inputs$distr
  fitted.values0 <- object$outputs$fitted.values
  y <- object$outputs$response
  loglik0_i <- do.call(what = distr,
                       args = c(list(x = y, log = TRUE), fitted.values0))

  # Saturated model
  saturated_model <- saturated_maxlogL(object)
  fitted.valuesS <- saturated_model$outputs$fitted.values
  loglikS_i <- do.call(what = distr,
                       args = c(list(x = y, log = TRUE), fitted.valuesS))
  output <- list(deviance_i = 2 * (loglikS_i - loglik0_i),
                 proposed_deviance_i = loglik0_i,
                 saturated_deviance_i = loglikS_i)
  return(output)
}
#==============================================================================
# Computation of expected value (and high order moments) for a fitted model ---
#==============================================================================
moment_integration_maxlogL <- function(object, ord = 1, routine,
                               ...){
  parameters <- object$outputs$fitted.values
  par_names <- names(parameters)
  parameters <- matrix(unlist(parameters), nrow = object$outputs$n)
  colnames(parameters) <- par_names
  distr <- object$inputs$distr
  support <- object$support$type

  integrand <- function(distr){
    nm_distr <- as.name(distr)
    pars <- formals(args(distr))
    log_par <- pars[names(pars) == 'log']
    pars <- sapply(object$outputs$par_names, as.name)
    distr_call <- as.call(c(nm_distr, as.name('x'), pars, log_par))
    body_fun <- str2expression(paste('x ^', ord, '*', deparse(distr_call)))
    func <- function() 'body'
    formals(func) <- formals(args(distr))
    formals(func)[object$outputs$par_names] <-
      sapply(object$outputs$par_names, function(x) x <- bquote())
    body(func) <- body_fun
    return(func)
  }
  EX <- integrand(distr)

  if ( missing(routine) ){
    if (support$type == 'continuos'){
      routine <- 'integrate'
    } else if (support$type == 'discrete'){
      routine <- 'sumate'
    }
  }
  mean_computation <- function(x)
    do.call(what = 'integration',
            args = c(list(fun = EX, lower = support$interval[1],
                          upper = support$interval[2],
                          routine = routine, ...), x))
  result <- apply(parameters, MARGIN = 1, mean_computation)
  return(result)
}
#==============================================================================
# Computation of cumulative for a fitted model --------------------------------
#==============================================================================
hazard_integration_maxlogL <- function(object, routine, ...){
  parameters <- object$outputs$fitted.values
  par_names <- names(parameters)
  parameters <- matrix(unlist(parameters), nrow = object$outputs$n)
  colnames(parameters) <- par_names
  distr <- object$inputs$distr
  support <- object$support$type

  integrand <- hazard_fun(distr, support)

  if ( missing(routine) ){
    if (support$type == 'continuos'){
      routine <- 'integrate'
    } else if (support$type == 'discrete'){
      routine <- 'summate'
    }
  }
  haz_computation <- function(x)
    do.call(what = 'integration',
            args = c(list(fun = integrand, lower = support$interval[1],
                          upper = support$interval[2],
                          routine = routine, ...), x))
  result <- apply(parameters, MARGIN = 1, haz_computation)
  return(result)
}