R/glm_logliks.R

In paulnorthrop/oola: Objected-oriented Loglikelihood Adjustment

# GLM (included so that a direct call to adjust_object() will work)

#' @export
logLikVec.glm <- function(object, pars = NULL, ...) {
  temp <- object
  # Add the model matrix, because it may not be present in x
  # This enables a direct call to adjust_object(), rather than using alogLik()
  temp$x <- model.matrix(object)
  class(temp) <- switch(object$family$family,
                        poisson = "poisglm",
                        binomial = "binomglm",
                        NULL)
  return(logLikVec(temp, pars = pars, ...))
}

# Poisson

#' @export
logLikVec.poisglm <- function(object, pars = NULL, ...) {
  if (object$family$family != "poisson") {
    stop("logLik.poisglm can only be used for objects in the Poisson family")
  }
  if (!missing(...)) {
    warning("extra arguments discarded")
  }
  # If the parameter estimates have been provided in pars then recalculate
  # the linear predictor
  if (!is.null(pars)) {
    object$linear.predictors <- drop(object$x %*% pars)
  }
  # offset
  off <- object$offset
  if (is.null(off)) {
    off <- 0
  }
  # linear predictor
  eta <- object$linear.predictors
  # Poisson mean.
  mu <- object$family$linkinv(off + eta)
  # responses
  y <- object$y
  # weights supplied to call to glm()
  w <- object$prior.weights
  # Calculate the weighted loglikelihood contributions
  val <- w * stats::dpois(x = y, lambda = mu, log = TRUE)
  # Remove attributes
  val <- as.vector(val)
  # Return the usual attributes for a "logLik" object
  attr(val, "nobs") <- sum(!is.na(object$residuals))
  attr(val, "df") <- object$rank
  class(val) <- "logLikVec"
  return(val)
}

# Binomial

#' @export
logLikVec.binomglm <- function(object, pars = NULL, ...) {
  if (object$family$family != "binomial") {
    stop("logLik.binomglm can only be used for objects in the Binomial family")
  }
  if (!missing(...)) {
    warning("extra arguments discarded")
  }
  # If the parameter estimates have been provided in pars then recalculate
  # the linear predictor
  if (!is.null(pars)) {
    object$linear.predictors <- drop(object$x %*% pars)
  }
  # offset
  off <- object$offset
  if (is.null(off)) {
    off <- 0
  }
  # linear predictor
  eta <- object$linear.predictors
  # Binomial mean
  mu <- object$family$linkinv(off + eta)
  # If the number of successes (and trials) are not available then create them
  # This enables a direct call to adjust_object(), rather than using alogLik()
  if (is.null(object$n_successes)) {
    response <- stats::model.response(stats::model.frame(object$call,
                                                         data = object$data))
    if (is.matrix(response)) {
      n_successes <- response[, 1]
      n_trials <- rowSums(response)
    } else {
      n_successes <- object$y
      n_trials <- 1
    }
  } else {
    n_successes <- object$n_successes
    n_trials <- object$n_trials
  }
  # weights supplied to call to glm()
  w <- object$prior.weights
  # Adjustment because prior.weights are used to store the number of trials
  w <- w / n_trials
  # Calculate the weighted loglikelihood contributions
  val <- w * stats::dbinom(x = n_successes, size = n_trials, prob = mu,
                           log = TRUE)
  # Remove attributes
  val <- as.vector(val)
  # Return the usual attributes for a "logLik" object
  attr(val, "nobs") <- sum(!is.na(object$residuals))
  attr(val, "df") <- object$rank
  class(val) <- "logLikVec"
  return(val)
}