conf_region: Two-dimensional confidence regions
In chandwich: Chandler-Bate Sandwich Loglikelihood Adjustment

conf_region

R Documentation

Two-dimensional confidence regions

Description

Calculates the (profile, if necessary) loglikelihood for a pair of parameters from which confidence regions can be plotted using plot.confreg.

Usage

conf_region(
  object,
  which_pars = NULL,
  range1 = c(NA, NA),
  range2 = c(NA, NA),
  conf = 95,
  mult = 2,
  num = c(10, 10),
  type = c("vertical", "cholesky", "spectral", "none"),
  ...
)

Arguments

`object`	An object of class `"chandwich"` returned by `adjust_loglik`.
`which_pars`	A vector of length 2 specifying the 2 (unfixed) parameters for which confidence region is 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 `which_pars` is not supplied but the current model has exactly two free parameters, i.e. `attr(object, "p_current") = 2` then `which_pars` is set to `attr(object, "free_pars") = 2`.
`range1, range2`	Numeric vectors of length 2. Respective ranges (of the form `c(lower, upper)`) of values of `which_pars[1]` and `which_pars[2]` over which to profile. Missing values in `range1` and/or `range2` are filled in using `conf` and `mult`. See below for details.
`conf`	A numeric scalar in (0, 100). The highest confidence level of interest. This is only relevant if `range1` and/or `range2` are not completely specified. In that event `conf` is used, in combination with `mult`, to try to set up the grid of parameter values to include the largest confidence region of interest.
`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 vector of length 1 or 2. The numbers of values at which to evaluate the profile loglikelihood either side of the MLE. `num[i]` relates to `which_pars[i]`. If `num` has length 1 then `num` is replicated to have length 2.
`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`. Any arguments that are not appropriate for `optim`, i.e. not in `methods::formalArgs(stats::optim)`, will be removed without warning.

Value

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

`grid1, grid2`	Numeric vectors. Respective values of `which_pars[1]` and `which_pars[2]` in the grid over which the (profile) loglikelihood is evaluated.
`max_loglik`	A numeric scalar. The value value of the loglikelihood at its maximum.
`prof_loglik`	An 2 `num` + 1 by 2 `num` + 1 numeric matrix containing the values of the (profile) loglikelihood.
`type`	A character scalar. The input `type`.
`which_pars`	A numeric or character vector. The input `which_pars`. If the `which_pars` was numeric then it is supplemented by the parameter names, if these are available in `object`.
`name`	A character scalar. The name of the model, stored in `attr(object, "name")`.

Examples

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


# Plots like those in Figure 4 of Chandler and Bate (2007)
# (a)
which_pars <- c("mu[0]", "mu[1]")
reg_1 <- conf_region(large, which_pars = which_pars)
reg_none_1 <- conf_region(large, which_pars = which_pars, type = "none")
plot(reg_1, reg_none_1)
# (b)
which_pars <- c("sigma[0]", "sigma[1]")
reg_2 <- conf_region(large, which_pars = which_pars)
reg_none_2 <- conf_region(large, which_pars = which_pars, type = "none")
plot(reg_2, reg_none_2)
# (c)
# Note: the naive and bivariate model contours are the reversed in the paper
which_pars <- c("sigma[0]", "xi[0]")
reg_3 <- conf_region(large, which_pars = which_pars)
reg_none_3 <- conf_region(large, which_pars = which_pars, type = "none")
plot(reg_3, reg_none_3)


# --------- 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")
# Linear model (gamma fixed at 0)
pois_lin <- adjust_loglik(pois_glm_loglik, y = y, x = x, par_names = pars,
                          fixed_pars = "gamma")
pois_vertical <- conf_region(pois_lin)
pois_none <- conf_region(pois_lin, type = "none")
plot(pois_none, pois_vertical, conf = c(50, 75, 95, 99), col = 2:1, lwd = 2,
     lty = 1)

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