alogLik: Loglikelihood adjustment of model fits
In paulnorthrop/oola: Objected-oriented Loglikelihood Adjustment

Description Usage Arguments Details Value Methods (by class) References See Also Examples

Description

alogLik(x, ...)

## Default S3 method:
alogLik(x, cluster = NULL, use_sandwich = TRUE,
  use_vcov = TRUE, ...)

## S3 method for class 'glm'
alogLik(x, cluster = NULL, use_sandwich = TRUE,
  use_vcov = TRUE, ...)

`x`	A fitted model object with certain associated S3 methods. See Details.
`...`	Further arguments to be passed to the functions in the sandwich package `meat`, if `cluster = NULL`, or `meatCL`, otherwise.
`cluster`	A vector or factor indicating from which cluster the respective loglikelihood contributions from `loglik` originate. This must have the same length as the vector returned by the `logLikVec` method for object like `x`. If `cluster` is not supplied then it is assumed that each observation forms its own cluster. If the sandwich package (Zeleis, 2006) is used to estimate the quantities required to adjust the loglikelihood (i.e. `use_sanswich = TRUE`) then `cluster` determines whether the variance matrix `V` of the score vector is estimated using `meat` (`cluster` is `NULL`) or `meatCL` (`cluster` is not `NULL`). See `use_sandwich` and Details.
`use_sandwich`	A logical scalar. Should we use the `sandwich` package (Zeleis, 2006) to estimate the variance V of the score function to be passed as the argument `V` to `adjust_loglik`? Otherwise, V is based on numerical derivatives, estimated using the `numDeriv` package. See Details for more information. The main purpose of `use_sandwich` is to enable a check that equivalent results are obtained using the sandwich package and the numerical derivatives.
`use_vcov`	A logical scalar. Should we use the `vcov` S3 method for `x` (if this exists) to estimate the Hessian of the independence loglikelihood to be passed as the argument `H` to `adjust_loglik`? Otherwise, `H` is estimated inside `adjust_loglik` using `optimHess`.

Supported classes, and the functions from which they are returned are:

"lm": lm in the stats package,
"glm": glm in the stats package

Objects of other classes are supported provided that they have certain S3 methods. The class of the object x must have the following S3 methods:

logLikVec: returns a vector of the contributions to the independence loglikelihood from individual observations;
coef: returns a vector of model coefficients, see coef;
nobs: returns the number of (non-missing) observations used in a model fit, see nobs;

and may have the following S3 methods

vcov: returns the estimated variance-covariance matrix of the (main) parameters of a fitted model, see vcov;
estfun: returns an n x k matrix, in which each column gives the derivative of the loglikelihood at each of n observation with respect to the k parameters of the model, see estfun.

If a vcov method is not available then the variance-covariance matrix is estimated inside adjust_loglik. If an estfun method is not available then the matrix is estimated using jacobian.

Loglikelihood adjustment is performed using the adjust_loglik function in the chandwich package. The relevant arguments to adjust_loglik, namely loglik, mle, H and V, are created based on the class of the object x. If use_sandwich = TRUE then H is inferred using bread in the sandwich package and, similarly, V is inferred using either meat, if cluster = NULL and meatCL, otherwise.

If cluster is NULL then arguments of meatCL present in ... will be ignored. Similarly, if cluster is not NULL then arguments of meat present in ... will be ignored.

An object of class inheriting from "oola", which inherits from the class "chandwich". See adjust_loglik. The attribute "name" of the returned object is the elements of class(x) concatenated into a character scalar and separated by _.

default: The default method used for user-supplied objects. These objects must have S3 methods logLikVec, nobs and coef.
glm: Method for objects returned from glm inheriting from "glm".

Chandler, R. E. and Bate, S. (2007). Inference for clustered data using the independence loglikelihood. Biometrika, 94(1), 167-183. http://dx.doi.org/10.1093/biomet/asm015

Zeleis (2006) Object-Oriented Computation and Sandwich Estimators. Journal of Statistical Software, 16, 1-16. http://dx.doi.org/10.18637/jss.v016.i09

summary.chandwich, plot.chandwich, confint.chandwich, anova.chandwich, coef.chandwich, vcov.chandwich and logLik.chandwich for S3 methods for objects of class "chandwich".

adjust_loglik to adjust a user-supplied loglikelihood.

bread, meat, meatCL and sandwich in the sandwich package.

# 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)
}
# 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)