Nothing
#' @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)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.