R/Anova.resi.R

In RESI: Robust Effect Size Index (RESI) Estimation

#' Anova method for resi objects
#'
#' After running the \code{\link{resi}} function on a fitted model, this function can be used to print the Anova-style table component. If the resi function was run with the `store.boot = TRUE` option to store the full matrix of bootstrapped estimates, the user can specify a different alpha level for this function's confidence intervals.
#' @param mod an object resulting from resi function
#' @param alpha an optional new specification for the confidence level. Can be vector-valued
#' @param ... ignored
#' @return Returns an `anova` object containing the computed Anova-style table
#' @examples
#' # fit a model
#' mod = lm(charges ~ bmi + sex, data = RESI::insurance)
#'
#' # run resi with the store.boot = TRUE option
#' resi.obj = resi(mod, nboot = 100, store.boot = TRUE, alpha = 0.01)
#'
#' # run Anova, specifying a different alpha level if desired
#' car::Anova(resi.obj, alpha = 0.05)
#' @noRd
#' @export
Anova.resi <- function(mod, alpha = NULL, ...){
  if(is.null(mod$anova)){
    stop('\nresi function was not run with anova = TRUE option')
  }

  if (is.null(alpha)){
    output = mod$anova
  }
  else{
    if (!(all(alpha %in% mod$alpha))){
      if (!is.null(mod$ci.method) && mod$ci.method != "boot" &&
          !is.null(mod$model.full)) {
        # Asymptotic (qf/normal): recompute anova CIs analytically.
        output <- mod$anova[, 1:(which(colnames(mod$anova) == 'RESI')), drop = FALSE]
        for (a in alpha) {
          asym_out <- tryCatch(
            resi_pe_asymptotic(
              model.full   = mod$model.full,
              alpha        = a,
              ci.method    = mod$ci.method,
              vcovfunc     = if (!is.null(mod$vcovfunc)) mod$vcovfunc else sandwich::vcovHC,
              vcov.args    = if (!is.null(mod$vcov.args)) mod$vcov.args else list(),
              Anova.args   = if (!is.null(mod$Anova.args)) mod$Anova.args else list(),
              unbiased     = if (!is.null(mod$unbiased)) mod$unbiased else TRUE,
              coefficients = FALSE,
              anova        = TRUE
            ),
            error = function(e) {
              stop(paste0("\nFailed to recompute asymptotic CIs at alpha = ", a,
                          ":\n", conditionMessage(e)))
            }
          )
          ci_cols <- paste0(c(a / 2, 1 - a / 2) * 100, '%')
          rn_match <- intersect(rownames(output), rownames(asym_out$anova))
          output[rn_match, ci_cols] <- asym_out$anova[rn_match, ci_cols]
        }
        return(output)
      } else if (is.null(mod$boot.results)) {
        stop('\nresi function was not run with store.boot = TRUE option')
      }
    }
    if(is.null(mod$boot.results)){
      output = mod$anova[c(1:(which(colnames(mod$anova) == 'RESI')),
                              which(colnames(mod$anova)%in%
                                      c(paste(alpha/2*100, '%', sep=''),
                                        paste((1-alpha/2)*100, '%', sep=''))))]
    }
    else{
      output = mod$anova[,1:(which(colnames(mod$anova) == 'RESI'))]
      boot.results = mod$boot.results$t
      CIs = apply(boot.results[,(ncol(boot.results)-
                                          nrow(mod$anova[which(rownames(mod$anova)
                                                                  != "Residuals"),])+1):
                                        ncol(boot.results)], 2,  quantile,
                  probs = sort(c(alpha/2, 1-alpha/2)), na.rm = TRUE)
      CIs = t(CIs)
      output[1:nrow(CIs), c(paste(alpha/2*100, '%', sep=''),
                            paste((1-rev(alpha)/2)*100, '%', sep=''))] = CIs
    }
  }
  return(output)
}

#' Anova method for resi objects
#'
#' After running the \code{\link{resi}} function on a fitted model, this function can be used to print the Anova-style table component. If the resi function was run with the `store.boot = TRUE` option to store the full matrix of bootstrapped estimates, the user can specify a different alpha level for this function's confidence intervals.
#' @param object an object resulting from resi function
#' @param alpha an optional new specification for the confidence level. Can be vector-valued
#' @param ... ignored
#' @return Returns an `anova` object containing the computed Anova-style table
#' @details The resi function uses the car::Anova function to compute the Anova table.
#'@examples
#' # fit a model
#' mod = lm(charges ~ bmi + sex, data = RESI::insurance)
#'
#' # run resi with the store.boot = TRUE option
#' resi.obj = resi(mod, nboot = 100, store.boot = TRUE, alpha = 0.01)
#'
#' # run anova, specifying a different alpha level if desired
#' anova(resi.obj, alpha = 0.05)
#' @export
anova.resi <- function(object, alpha = NULL, ...){
  if(is.null(object$anova)){
    stop("\nresi function was not run with anova = TRUE option")
  }

  if (is.null(alpha)){
    output = object$anova
  }
  else{
    if (!(all(alpha %in% object$alpha))){
      if (!is.null(object$ci.method) && object$ci.method != "boot" &&
          !is.null(object$model.full)) {
        # Asymptotic (qf/normal): recompute anova CIs analytically.
        output <- object$anova[, 1:(which(colnames(object$anova) == "RESI")), drop = FALSE]
        for (a in alpha) {
          asym_out <- tryCatch(
            resi_pe_asymptotic(
              model.full   = object$model.full,
              alpha        = a,
              ci.method    = object$ci.method,
              vcovfunc     = if (!is.null(object$vcovfunc)) object$vcovfunc else sandwich::vcovHC,
              vcov.args    = if (!is.null(object$vcov.args)) object$vcov.args else list(),
              Anova.args   = if (!is.null(object$Anova.args)) object$Anova.args else list(),
              unbiased     = if (!is.null(object$unbiased)) object$unbiased else TRUE,
              coefficients = FALSE,
              anova        = TRUE
            ),
            error = function(e) {
              stop(paste0("\nFailed to recompute asymptotic CIs at alpha = ", a,
                          ":\n", conditionMessage(e)))
            }
          )
          ci_cols <- paste0(c(a / 2, 1 - a / 2) * 100, "%")
          rn_match <- intersect(rownames(output), rownames(asym_out$anova))
          output[rn_match, ci_cols] <- asym_out$anova[rn_match, ci_cols]
        }
        return(output)
      } else if (is.null(object$boot.results)) {
        stop("\nresi function was not run with store.boot = TRUE option")
      }
    }
    if(is.null(object$boot.results)){
      output = object$anova[c(1:(which(colnames(object$anova) == "RESI")),
                              which(colnames(object$anova)%in%
                                      c(paste(alpha/2*100, "%", sep=""),
                                        paste((1-alpha/2)*100, "%", sep=""))))]
    }
    else{
      output = object$anova[,1:(which(colnames(object$anova) == "RESI"))]
      boot.results = object$boot.results$t
      CIs = apply(boot.results[,(ncol(boot.results)-
                                          nrow(object$anova)+1):
                                        ncol(boot.results)], 2,
                  quantile, probs = sort(c(alpha/2, 1-alpha/2)), na.rm = TRUE)
      CIs = t(CIs)
      output[1:nrow(CIs), c(paste(alpha/2*100, "%", sep=""),
                            paste((1-rev(alpha)/2)*100, "%", sep=""))] = CIs
    }
  }
  return(output)
}