Compares nested models using the adjusted likelihood ratio test statistic (ALRTS) described in Section 3.5 of Chandler and Bate (2007). The nesting must result from the simple constraint that a subset of the parameters of the larger model is held fixed.
An object of class "chandwich" returned by
adjust_loglik. The larger of the two models.
An object of class "chandwich" returned by
adjust_loglik. The smaller of the two models.
If smaller is supplied then the arguments fixed_pars and
fixed_at described below are ignored.
A logical scalar. If approx = TRUE then the
approximation detailed by equations (18)-(20) of Chandler and Bate (2007)
is used. This option is available only if smaller is supplied.
If smaller is not supplied then approx = TRUE is used,
with no warning.
The approximation doesn't make sense if type = "none". If
type = "none" and approx = TRUE then approx is
set to FALSE with no warning.
If approx = FALSE then the adjusted likelihood is
maximised under the restriction imposed by delta, that is,
equation (17) of Chandler and Bate (2007) is used.
A character scalar. The argument type to the function
returned by adjust_loglik, that is, the type of adjustment
made to the independence loglikelihood function.
A numeric vector. Indices of the components of the
full parameter vector that are restricted to be equal to the
value(s) in fixed_at.
A numeric vector of length 1 or length(fixed_pars).
If length(fixed_at) = 1 then the components fixed_pars
of the parameter vector are all fixed at fixed_at.
If length(fixed_at) = length(fixed_pars) then the component
fixed_pars[i] is fixed at fixed_at[i] for each i.
(Only relevant if approx = FALSE).
A numeric vector of initial values for use in the search for
the MLE under the smaller model. Must have length equal to the number
of parameters in the smaller of the two models being compared. If
init is not supplied, or it is of the wrong length, then
attr(smaller, "MLE") is used if smaller is supplied and
attr(larger, "MLE")[-fixed_pars] otherwise.
Further arguments to be passed to optim.
These may include gr, method, lower, upper
or control.
An object of class "compmod", a list with components
the adjusted likelihood ratio test statistic.
under the null hypothesis that the smaller model is a valid
simplification of the larger model the adjusted likelihood ratio
test statistic has a chi-squared distribution with df
degrees of freedom.
the p-value associated with the test of the null hypothesis.
the MLE of the parameters under the larger model.
the MLE of the parameters under the smaller model.
Numeric vectors of the indices of parameters fixed in the larger and smaller models, respectively.
Numeric vectors of the
values at which the parameters in larger_fixed_pars and
smaller_fixed_pars are fixed.
the argument approx.
The smaller of the two models is specified either by supplying
smaller or fixed_pars. If both are supplied then
smaller takes precedence.
For full details see Section 3.5 of Chandler and Bate (2007).
If approx = FALSE then the a likelihood ratio test of the null
hypothesis that the smaller model is a valid simplification of the larger
model is carried out directly using equation (17) of Chandler and Bate
(2007) based on the adjusted loglikelihood under the larger model,
returned by adjust_loglik. This adjusted loglikelihood is
maximised subject to the constraint that a subset of the parameters
in the larger model are fixed. If smaller is supplied
then this maximisation can be avoided using an approximation
detailed by equations (18)-(20) of Chandler and Bate (2007), which uses
the MLE under the smaller model. The same null distribution (chi-squared
with degrees of freedom equal to the number of parameters that are fixed)
is used in both cases.
Chandler, R. E. and Bate, S. (2007). Inference for clustered data using the independence loglikelihood. Biometrika, 94(1), 167-183. doi:10.1093/biomet/asm015
adjust_loglik to adjust a user-supplied
loglikelihood function.
conf_intervals for confidence intervals for
individual parameters.
conf_region for a confidence region for
pairs of parameters.
# -------------------------- GEV model, owtemps data -----------------------
# ------------ following Section 5.2 of Chandler and Bate (2007) -----------
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)
# Log-likelihood adjustment of 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)
# Log-likelihood adjustment of some smaller models: xi[1] = 0 etc
medium <- adjust_loglik(larger = large, fixed_pars = "xi[1]")
small <- adjust_loglik(larger = medium, fixed_pars = c("sigma[1]", "xi[1]"))
# Tests
# Test xi1 = 0 (2 equivalent ways), vertical adjustment
compare_models(large, fixed_pars = "xi[1]")
#> Model: gev_loglik
#> H0: "xi[1]" = 0
#> HA: unrestricted model
#>
#> test statistic = 6.356, df = 1, p-value = 0.0117
compare_models(large, medium)
#> Model: gev_loglik
#> H0: "xi[1]" = 0
#> HA: unrestricted model
#>
#> test statistic = 6.356, df = 1, p-value = 0.0117
# Test xi1 = 0, using approximation
compare_models(large, medium, approx = TRUE)
#> Model: gev_loglik
#> H0: "xi[1]" = 0
#> HA: unrestricted model
#>
#> Using using approximation (18) of Chandler and Bate (2007):
#> test statistic = 5.245, df = 1, p-value = 0.022
# Horizontal adjustments
compare_models(large, medium, type = "cholesky")$p_value
#> [1] 0.01593336
compare_models(large, medium, type = "spectral")$p_value
#> [1] 0.01229833
# No adjustment (independence loglikelihood)
compare_models(large, medium, type = "none")$p_value
#> [1] 0.074199
# Test sigma1 = 0 for model with xi1 = 0
compare_models(medium, small)
#> Model: gev_loglik
#> H0: ("sigma[1]", "xi[1]") = (0, 0)
#> HA: "xi[1]" = 0
#>
#> test statistic = 4.251, df = 1, p-value = 0.03924
# Test sigma1 = xi1 = 0
compare_models(large, small)
#> Model: gev_loglik
#> H0: ("sigma[1]", "xi[1]") = (0, 0)
#> HA: unrestricted model
#>
#> test statistic = 8.6, df = 2, p-value = 0.01357
# --------- Misspecified Poisson model for negative binomial data ----------
# ... following Section 5.1 of the "Object-Oriented Computation of Sandwich
# Estimators" vignette of the sandwich package
# https://cran.r-project.org/web/packages/sandwich/vignettes/sandwich-OOP.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
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 1.06326821 0.04135723 25.709367 9.184267e-146
#> x 0.99607219 0.05353446 18.606186 2.862861e-77
#> I(x^2) -0.04912373 0.02314608 -2.122335 3.380961e-02
# 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)
#> MLE SE adj. SE
#> alpha 1.06300 0.04136 0.08378
#> beta 0.99610 0.05354 0.10520
#> gamma -0.04913 0.02315 0.03628
pois_lin <- adjust_loglik(larger = pois_quad, fixed_pars = "gamma")
# Test the significance of the quadratic term
compare_models(pois_quad, pois_lin)$p_value
#> [1] 0.1772841
compare_models(pois_quad, pois_lin, approx = TRUE)$p_value
#> [1] 0.1657919