# 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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.