conf_intervals: Confidence intervals
In chandwich: Chandler-Bate Sandwich Loglikelihood Adjustment

conf_intervals

R Documentation

Confidence intervals

Description

Calculates confidence intervals for individual parameters.

Usage

conf_intervals(
  object,
  which_pars = NULL,
  init = NULL,
  conf = 95,
  mult = 1.5,
  num = 10,
  type = c("vertical", "cholesky", "spectral", "none"),
  profile = TRUE,
  ...
)

Arguments

`object`	An object of class `"chandwich"` returned by `adjust_loglik`.
`which_pars`	A vector specifying the (unfixed) parameters for which confidence intervals are required. Can be either a numeric vector, specifying indices of the components of the full parameter vector, or a character vector of parameter names, which must be a subset of those supplied in `par_names` in the call to `adjust_loglik` that produced `object`. `which_pars` must not have any parameters in common with `attr(object, "fixed_pars")`. `which_pars` must not contain all of the unfixed parameters, i.e. there is no point in profiling over all the unfixed parameters. If missing, all parameters are included.
`init`	A numeric vector of initial estimates of the values of the parameters that are not fixed and are not in `which_pars`. Should have length `attr(object, "p_current") - length(which_pars)`. If `init` is `NULL` or is of the wrong length then the relevant components from the MLE stored in `object` are used.
`conf`	A numeric scalar in (0, 100). Confidence level for the intervals.
`mult`	A numeric vector of length 1 or the same length as `which_pars`. The search for the profile loglikelihood-based confidence limits is conducted over the corresponding symmetric confidence intervals (based on approximate normal theory), extended by a factor of the corresponding component of `mult`.
`num`	A numeric scalar. The number of values at which to evaluate the profile loglikelihood either side of the MLE. Increasing `num` increases the accuracy of the confidence limits, but the code will take longer to run.
`type`	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.
`profile`	A logical scalar. If `FALSE` then only intervals based on approximate large sample normal theory, which are symmetric about the MLE, are returned (in `sym_CI`) and `prof_CI` in the returned object will contain `NA`s.
`...`	Further arguments to be passed to `optim`. These may include `gr`, `method`, `lower`, `upper` or `control`.

Details

Calculates (profile, if necessary) likelihood-based confidence intervals for individual parameters, and also provides symmetric intervals based on a normal approximation to the sampling distribution of the estimator. See also the S3 confint method confint.chandwich.

Value

An object of class "confint", a list with components

`conf`	The argument `conf`.
`cutoff`	A numeric scalar. For values inside the confidence interval the profile loglikelihood lies above `cutoff`.
`parameter_vals, prof_loglik_vals`	`2 * num + 1` by `length(which_pars)` numeric matrices. Column i of `parameter_vals` contains the profiled values of parameter `which_par[i]`. Column i of `prof_loglik_vals` contains the corresponding values of the profile loglikelihood.
`sym_CI, prof_CI`	`length(which_pars)` by 2 numeric matrices. Row i of `sym_CI` (`prof_CI`) contains the symmetric (profile loglikelihood-based) confidence intervals for parameter `which_pars[i]`.

If a value in prof_CI is NA then this means that the search for the confidence limit did no extend far enough. A remedy is to increase the value of mult, or the relevant component of mult, and perhaps also increase num.

`max_loglik`	The value of the adjusted loglikelihood at its maximum, stored in `object$max_loglik`.
`type`	The argument `type` supplied in the call to `conf_intervals`, i.e. the type of loglikelihood adjustment.
`which_pars`	The argument `which_pars`.
`name`	A character scalar. The name of the model, stored in `attr(object, "name")`.
`p_current`	The number of free parameters in the current model.
`fixed_pars, fixed_at`	`attr(object, "fixed_pars")` and `attr(object, "fixed_at")`, the arguments `fixed_pars` and `fixed_at` to `adjust_loglik`, if these were supplied.

Examples

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

# 95% likelihood-based confidence intervals, vertically adjusted
ci <- conf_intervals(rat_res)
plot(ci)
# Unadjusted
conf_intervals(rat_res, type = "none")

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

# 95% likelihood-based confidence intervals, vertically adjusted
large_v <- conf_intervals(large, which_pars = c("xi[0]", "xi[1]"))
large_v
plot(large_v)
plot(large_v, which_par = "xi[1]")

# Unadjusted
large_none <- conf_intervals(large, which_pars = c("xi[0]", "xi[1]"),
                             type = "none")
large_none
plot(large_v, large_none)
plot(large_v, large_none, which_par = "xi[1]")


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

# 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)
conf_intervals(pois_quad)

chandwich documentation built on Aug. 26, 2023, 1:07 a.m.