R/alogLik_methods.R

In paulnorthrop/oola: Objected-oriented Loglikelihood Adjustment

#' Loglikelihood adjustment of model fits
#'
#' Description
#'
#' @inherit adjust_object params details return references seealso
#' @export
alogLik <- function(x, ...) {
  UseMethod("alogLik")
}

#' @inheritParams adjust_object
#' @describeIn alogLik The default method used for user-supplied objects.
#'   These objects must have S3 methods \code{\link{logLikVec}},
#'   \code{\link[stats]{nobs}} and \code{\link[stats]{coef}}.
#' @examples
#' # An example of a use of alogLik.default based on the evd::fgev() function
#' # In the background are S3 method functions for objects inheriting from
#' # class "gev": loglikVec.gev(), coef.gev(), nobs.gev(), vcov.gev()
#'
#' # We need the evd package
#' got_evd <- requireNamespace("evd", quietly = TRUE)
#'
#' if (got_evd) {
#'   library(evd)
#'   # An example from the evd::fgev documentation
#'   uvdata <- evd::rgev(100, loc = 0.13, scale = 1.1, shape = 0.2)
#'
#'   M1 <- evd::fgev(uvdata, nsloc = (-49:50)/100)
#'   adj_fgev <- alogLik(M1)
#'   summary(adj_fgev)
#'
#'   M2 <- evd::fgev(uvdata, nsloc = (-49:50)/100, shape = 0)
#'   adj_fgev <- alogLik(M2)
#'   summary(adj_fgev)
#'
#'   # An example from Chandler and Bate (2007)
#'   y <- c(chandwich::owtemps[, "Oxford"], chandwich::owtemps[, "Worthing"])
#'   x <- rep(c(-1, 1), each = length(y) / 2)
#'   owfit <- evd::fgev(y, nsloc = x)
#'   year <- rep(1:(length(y) / 2), 2)
#'   adj_owfit <- alogLik(owfit, cluster = year)
#'   summary(adj_owfit)
#' }
#' @export
alogLik.default <- function(x, cluster = NULL, use_sandwich = TRUE,
                            use_vcov = TRUE, ...) {
  adjust_object(x, cluster = cluster, use_sandwich = use_sandwich,
                use_vcov = use_vcov, ...)
}

#' Loglikelihood adjustment of Generalized Linear Model (GLM) fits
#'
#' @inheritParams adjust_object
#' @describeIn alogLik Method for objects returned from
#'   \code{\link[stats]{glm}} inheriting from \code{"glm"}.
#' @examples
#' # Poisson GLM ----------
#'
#' # Example from the help file for stats::glm()
#' counts <- c(18,17,15,20,10,20,25,13,12)
#' outcome <- gl(3,1,9)
#' treatment <- gl(3,3)
#' d.AD <- data.frame(treatment, outcome, counts)
#' glm.D93 <- glm(counts ~ outcome + treatment, family = poisson, data = d.AD)
#' adj_pois <- adjust_object(glm.D93)
#' summary(adj_pois)
#'
#' # Binomial GLMs ----------
#'
#' # An example from Section 5.2 of sandwich vignette at
#' # https://cran.r-project.org/web/packages/sandwich/vignettes/sandwich-OOP.pdf
#'
#' got_AER <- requireNamespace("AER", quietly = TRUE)
#' if (got_AER) {
#'   data("Affairs", package = "AER")
#'   fm_probit <- glm(I(affairs > 0) ~ age + yearsmarried + religiousness +
#'                    occupation + rating, data = Affairs,
#'                    family = binomial(link = "probit"))
#'   adj_binom <- alogLik(fm_probit)
#'   summary(adj_binom)
#' }
#'
#' # An example where the response data are a (two-column) matrix
#' # (number of successes, number of failures)
#' lrfit <- glm(cbind(using, notUsing) ~ age + education + wantsMore,
#'              family = binomial, data = cuse)
#' adj_binom_2 <- alogLik(lrfit)
#' summary(adj_binom_2)
#' @export
alogLik.glm <- function(x, cluster = NULL, use_sandwich = TRUE,
                        use_vcov = TRUE, ...) {
  name_of_class <- switch(x$family$family,
                          poisson = "poisglm",
                          binomial = "binomglm",
                          NULL)
  # Add the model matrix, because it may not be present in x
  x$x <- model.matrix(x)
  if (name_of_class == "binomglm") {
    response <- stats::model.response(stats::model.frame(x$call,
                                                         data = x$data))
    if (is.matrix(response)) {
      x$n_successes <- response[, 1]
      x$n_trials <- rowSums(response)
    } else {
      x$n_successes <- x$y
      x$n_trials <- 1
    }
  }
  class(x) <- c(name_of_class, "glm")
  adjust_object(x, cluster = cluster, use_sandwich = use_sandwich,
                use_vcov = use_vcov, ...)
}

#' Loglikelihood adjustment of Linear Model fits
#'
#' @inheritParams adjust_object
#' @examples
#' data(PetersenCL, package = "sandwich")
#' p_lm <- lm(y ~ x, data = PetersenCL)
#' #adj_p_lm <- alogLik(p_lm) (doesn't work yet)
#' @export
alogLik.lm <- function(x, cluster = NULL, use_sandwich = TRUE,
                      use_vcov = TRUE, ...) {
  adjust_object(x, cluster = cluster, use_sandwich = use_sandwich,
                use_vcov = use_vcov, ...)
}