Description Usage Arguments Details Value References See Also Examples
View source: R/adjust_loglik.R
Performs adjustments of a usersupplied independence loglikelihood for the presence of cluster dependence, following Chandler and Bate (2007). The user provides a function that returns a vector of observationspecific loglikelihood contributions and a vector that indicates cluster membership. The loglikelihood of a submodel can be adjusted by fixing a set of parameters at particular values.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 
loglik 
A named function. Returns a vector of the
loglikelihood contributions of individual observations. The first
argument must be the vector of model parameter(s). If any of the model
parameters are outofbounds then 
... 
Further arguments to be passed either to 
cluster 
A vector or factor indicating from which cluster the
respective loglikelihood contributions from 
p 
A numeric scalar. The dimension of the full parameter
vector, i.e. the number of parameters in the full model. Must be
consistent with the lengths of 
init 
A numeric vector of initial values. Must have length equal
to the number of parameters in the full model. If 
par_names 
A character vector. Names of the 
fixed_pars 
A vector specifying which parameters are to be restricted
to be equal to the value(s) in 
fixed_at 
A numeric vector of length 1 or 
name 
A character scalar. A name for the model that gives rise
to 
larger 
Only relevant if An object of class 
alg_deriv 
A function with the vector of model parameter(s) as its
first argument. Returns a 
alg_hess 
A function with the vector of model parameter(s) as its
first argument. Returns a Supplying both 
mle 
A numeric vector. Can only be used if 
H, V 
p by p numeric matrices. Only used if Supplying both 
Three adjustments to the independence loglikelihood described in
Chandler and Bate (2007) are available. The vertical adjustment is
described in Section 6 and two horizontal adjustments are described
in Sections 3.2 to 3.4. See the descriptions of type
and, for the
horizontal adjustments, the descriptions of C_cholesky
and
C_spectral
, in Value.
The adjustments involve first and second derivatives of the loglikelihood
with respect to the model parameters. These are estimated using
jacobian
and optimHess
unless alg_deriv
and/or alg_hess
are supplied.
A function of class "chandwich"
to evaluate an adjusted
loglikelihood, or the independence loglikelihood, at one or more sets
of model parameters, with arguments
x 
A numeric vector or matrix giving values of the 
type 
A character scalar. The type of adjustment to use.
One of 
The latter results in the evaluation of the (unadjusted) independence loglikelihood. The function has (additional) attributes
p_full, p_current 
The number of parameters in the full model and current models, respectively. 
free_pars 
A numeric vector giving the indices of the free
parameters in the current model, with names inferred from

MLE, res_MLE 
Numeric vectors, with names inferred from

SE, adjSE 
The unadjusted and adjusted estimated standard errors, of the free parameters, respectively. 
VC, adjVC 
The unadjusted and adjusted estimated variancecovariance matrix of the free parameters, respectively. 
HI, HA 
The Hessians of the independence and adjusted loglikelihood, respectively. 
C_cholesky, C_spectral 
The matrix C in equation (14) of Chandler and Bate (2007), calculated using Cholesky decomposition and spectral decomposition, respectively. 
full_par_names, par_names 
The names of the parameters in the full and current models, respectively, if these were supplied in this call or a previous call. 
max_loglik 
The common maximised value of the independence and adjusted loglikelihoods. 
loglik, cluster 
The arguments 
loglik_args 
A list containing the further arguments passed to

loglikVecMLE 
a vector containing the contributions of individual observations to the independence loglikelihood evaluated at the MLE. 
name 
The argument 
nobs 
The number of observations. 
call 
The call to 
If fixed_pars
is not NULL
then there are further attributes
fixed_pars 
The argument 
fixed_at 
The argument 
If alg_deriv
and/or alg_hess
were supplied then these are
returned as further attributes.
To view an individual attribute called att_name
use
attr(x, "att_name")
or attributes(x)$att_name
.
Chandler, R. E. and Bate, S. (2007). Inference for clustered data using the independence loglikelihood. Biometrika, 94(1), 167183. doi: 10.1093/biomet/asm015
summary.chandwich
for maximum likelihood estimates
and unadjusted and adjusted standard errors.
plot.chandwich
for plots of onedimensional adjusted
loglikelihoods.
confint.chandwich
, anova.chandwich
,
coef.chandwich
, vcov.chandwich
and logLik.chandwich
for other chandwich
methods.
conf_intervals
for confidence intervals for
individual parameters.
conf_region
for a confidence region for
a pair of parameters.
compare_models
to compare nested models using an
(adjusted) likelihood ratio test.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110  #  Binomial model, rats data 
# Contributions to the independence loglikelihood
binom_loglik < function(prob, data) {
if (prob < 0  prob > 1) {
return(Inf)
}
return(dbinom(data[, "y"], data[, "n"], prob, log = TRUE))
}
rat_res < adjust_loglik(loglik = binom_loglik, data = rats, par_names = "p")
# Plot the loglikelihoods
plot(rat_res, type = 1:4, legend_pos = "bottom", lwd = 2, col = 1:4)
# MLE, SEs and adjusted SEs
summary(rat_res)
#  GEV model, owtemps data 
#  following Section 5.2 of Chandler and Bate (2007) 
# Contributions to the independence loglikelihood
gev_loglik < function(pars, data) {
o_pars < pars[c(1, 3, 5)] + pars[c(2, 4, 6)]
w_pars < pars[c(1, 3, 5)]  pars[c(2, 4, 6)]
if (isTRUE(o_pars[2] <= 0  w_pars[2] <= 0)) return(Inf)
o_data < data[, "Oxford"]
w_data < data[, "Worthing"]
check < 1 + o_pars[3] * (o_data  o_pars[1]) / o_pars[2]
if (isTRUE(any(check <= 0))) return(Inf)
check < 1 + w_pars[3] * (w_data  w_pars[1]) / w_pars[2]
if (isTRUE(any(check <= 0))) return(Inf)
o_loglik < log_gev(o_data, o_pars[1], o_pars[2], o_pars[3])
w_loglik < log_gev(w_data, w_pars[1], w_pars[2], w_pars[3])
return(o_loglik + w_loglik)
}
# Initial estimates (method of moments for the Gumbel case)
sigma < as.numeric(sqrt(6 * diag(var(owtemps))) / pi)
mu < as.numeric(colMeans(owtemps)  0.57722 * sigma)
init < c(mean(mu), diff(mu) / 2, mean(sigma), diff(sigma) / 2, 0, 0)
# Loglikelihood adjustment for the full model
par_names < c("mu[0]", "mu[1]", "sigma[0]", "sigma[1]", "xi[0]", "xi[1]")
large < adjust_loglik(gev_loglik, data = owtemps, init = init,
par_names = par_names)
# Rows 1, 3 and 4 of Table 2 of Chandler and Bate (2007)
t(summary(large))
# Loglikelihood adjustment of some smaller models: xi[1] = 0 etc
# Starting from a larger model
medium < adjust_loglik(larger = large, fixed_pars = "xi[1]")
small < adjust_loglik(larger = large, fixed_pars = c("sigma[1]", "xi[1]"))
small < adjust_loglik(larger = medium, fixed_pars = c("sigma[1]", "xi[1]"))
# Starting from scratch
medium < adjust_loglik(gev_loglik, data = owtemps, init = init,
par_names = par_names, fixed_pars = "xi[1]")
small < adjust_loglik(gev_loglik, data = owtemps, init = init,
par_names = par_names, fixed_pars = c("sigma[1]", "xi[1]"))
#  Misspecified Poisson model for negative binomial data 
# ... following Section 5.1 of the "ObjectOriented Computation of Sandwich
# Estimators" vignette of the sandwich package
# https://cran.rproject.org/web/packages/sandwich/vignettes/sandwichOOP.pdf
# Simulate data
set.seed(123)
x < rnorm(250)
y < rnbinom(250, mu = exp(1 + x), size = 1)
# Fit misspecified Poisson model
fm_pois < glm(y ~ x + I(x^2), family = poisson)
summary(fm_pois)$coefficients
# Contributions to the independence loglikelihood
pois_glm_loglik < function(pars, y, x) {
log_mu < pars[1] + pars[2] * x + pars[3] * x ^ 2
return(dpois(y, lambda = exp(log_mu), log = TRUE))
}
pars < c("alpha", "beta", "gamma")
pois_quad < adjust_loglik(pois_glm_loglik, y = y, x = x, par_names = pars)
summary(pois_quad)
# Providing algebraic derivatives and Hessian
pois_alg_deriv < function(pars, y, x) {
mu < exp(pars[1] + pars[2] * x + pars[3] * x ^ 2)
return(cbind(y  mu, x * (y  mu), x ^2 * (y  mu)))
}
pois_alg_hess < function(pars, y, x) {
mu < exp(pars[1] + pars[2] * x + pars[3] * x ^ 2)
alg_hess < matrix(0, 3, 3)
alg_hess[1, ] < c(sum(mu), sum(x * mu), sum(x ^ 2 * mu))
alg_hess[2, ] < c(sum(x * mu), sum(x ^ 2 * mu), sum(x ^ 3 * mu))
alg_hess[3, ] < c(sum(x ^ 2 * mu), sum(x ^ 3 * mu), sum(x ^ 4 * mu))
return(alg_hess)
}
pois_quad < adjust_loglik(pois_glm_loglik, y = y, x = x, p = 3,
alg_deriv = pois_alg_deriv, alg_hess = pois_alg_hess)
summary(pois_quad)
got_sandwich < requireNamespace("sandwich", quietly = TRUE)
if (got_sandwich) {
# Providing MLE, H and V
# H and V calculated using bread() and meat() from sandwich package
n_obs < stats::nobs(fm_pois)
pois_quad < adjust_loglik(pois_glm_loglik, y = y, x = x, p = 3,
mle = fm_pois$coefficients,
H = solve(sandwich::bread(fm_pois) / n_obs),
V = sandwich::meat(fm_pois) * n_obs)
}

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.