conf_intervals: 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"), ...)

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.

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.

...

Further arguments to be passed to optim. These may include gr, method, lower, upper or control.

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].

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.

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.

Examples

Run this code

# NOT RUN {
# ------------------------- 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
conf_intervals(rat_res)
# 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 (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 (any(check <= 0)) return(-Inf)
  check <- 1 + w_pars[3] * (w_data - w_pars[1]) / w_pars[2]
  if (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]")
# }
# NOT RUN {
# 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]")
# }
# NOT RUN {
# --------- 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)
# }

Run the code above in your browser using DataCamp Workspace