enrichwith (version 0.3.1)

enrich.glm: Enrich objects of class glm

Description

Enrich objects of class glm with any or all of a set of auxiliary functions, the maximum likelihood estimate of the dispersion parameter, the expected or observed information at the maximum likelihood estimator, and the first term in the expansion of the bias of the maximum likelihood estimator.

Usage

# S3 method for glm
enrich(object, with = "all", ...)

Arguments

object

an object of class glm

with

a character vector of options for the enrichment of object

...

extra arguments to be passed to the compute_* functions

Value

The object object of class glm with extra components. See get_enrichment_options.glm() for the components and their descriptions.

Details

The auxiliary_functions component consists of any or all of the following functions:

• score: the log-likelihood derivatives as a function of the model parameters; see get_score_function.glm

• information: the expected or observed information as a function of the model parameters; see get_information_function.glm

• bias: the first-order term in the expansion of the bias of the maximum likelihood estimator as a function of the model parameters; see get_bias_function.glm

• simulate: a simulate function for glm objects that can simulate variates from the model at user-supplied parameter values for the regression parameters and the dispersion (default is the maximum likelihood estimates); see get_simulate_function.glm

• dmodel: computes densities or probability mass functions under the model at user-supplied data.frames and at user-supplied values for the regression parameters and the dispersion, if any (default is at the maximum likelihood estimates); see get_dmodel_function.glm

• pmodel: computes distribution functions under the model at user-supplied data.frames and at user-supplied values for the regression parameters and the dispersion, if any (default is at the maximum likelihood estimates); see get_pmodel_function.glm

• qmodel: computes quantile functions under the model at user-supplied data.frames and at user-supplied values for the regression parameters and the dispersion, if any (default is at the maximum likelihood estimates); see get_qmodel_function.glm

Examples

# NOT RUN {
# }
# NOT RUN {
# A Gamma example, from McCullagh & Nelder (1989, pp. 300-2)
clotting <- data.frame(
   u = c(5,10,15,20,30,40,60,80,100, 5,10,15,20,30,40,60,80,100),
   time = c(118,58,42,35,27,25,21,19,18,69,35,26,21,18,16,13,12,12),
   lot = factor(c(rep(1, 9), rep(2, 9))))
cML <- glm(time ~ lot*log(u), data = clotting, family = Gamma)

# The simulate method for the above fit would simulate at coef(cML)
# for the regression parameters and MASS::gamma.dispersion(cML) for
# the dispersion. It is not possible to simulate at different
# parameter values than those, at least not, without "hacking" the
# cML object.

# A general simulator for cML results via its enrichment with
# auxiliary functions:
cML_functions <- get_auxiliary_functions(cML)
# which is a shorthand for
# enriched_cML <- enrich(cML, with = "auxiliary functions")
# cML_functions <- enriched_cML$auxiliary_functions

# To simulate 2 samples at the maximum likelihood estimator do
dispersion_mle <- MASS::gamma.dispersion(cML)
cML_functions$simulate(coef = coef(cML),
                       dispersion = dispersion_mle,
                       nsim = 2, seed = 123)
# To simulate 5 samples at c(0.1, 0.1, 0, 0) and dispersion 0.2 do
cML_functions$simulate(coef = c(0.1, 0.1, 0, 0),
                       dispersion = 0.2,
                       nsim = 5, seed = 123)

# }
# NOT RUN {
# }
# NOT RUN {
## Reproduce left plot in Figure 4.1 in Kosimdis (2007)
## (see http://www.ucl.ac.uk/~ucakiko/files/ikosmidis_thesis.pdf)
mod <- glm(1 ~ 1, weights = 10, family = binomial())
enriched_mod <- enrich(mod, with = "auxiliary functions")
biasfun <- enriched_mod$auxiliary_functions$bias
probabilities <- seq(1e-02, 1 - 1e-02, length = 100)
biases <- Vectorize(biasfun)(qlogis(probabilities))
plot(probabilities, biases, type = "l", ylim = c(-0.5, 0.5),
     xlab = expression(pi), ylab = "first-order bias")
abline(h = 0, lty = 2); abline(v = 0.5, lty = 2)
title("First-order bias of the MLE of the log-odds", sub = "m = 10")
# }