Skip to contents

Confidence intervals for model parameters

Source: R/confidence_intervals.R
confint.chandwich.Rd

confint method for objects of class "chandwich". Computes confidence intervals for one or more model parameters based on an object returned from adjust_loglik.

Usage

# S3 method for class 'chandwich'
confint(
  object,
  parm,
  level = 0.95,
  type = c("vertical", "cholesky", "spectral", "none"),
  profile = TRUE,
  ...
)

Arguments

object

An object of class "chandwich", returned by adjust_loglik.

parm

A vector specifying the (unfixed) parameters for which confidence intervals are required. If missing, all parameters are included.

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.

parm must not have any parameters in common with attr(object, "fixed_pars"). parm must not contain all of the unfixed parameters, i.e. there is no point in profiling over all the unfixed parameters.

level

The confidence level required. A numeric scalar in (0, 1).

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 TRUE then confidence intervals based on an (adjusted) profile loglikelihood are returned. If FALSE then intervals based on approximate large sample normal theory, which are symmetric about the MLE, are returned.

...

Further arguments to be passed to conf_intervals.

Value

A matrix with columns giving lower and upper confidence limits for each parameter. These will be labelled as (1 - level)/2 and 1 - (1 - level)/2 in % (by default 2.5% and 97.5%). The row names are the names of the model parameters, if these are available.

Details

For details see the documentation for the function conf_intervals, on which confint.chandwich is based.

See also

The underlying function conf_intervals. If you would like to plot the profile loglikelihood function for a parameter then call conf_intervals directly and then use the associated plot method. Note that in conf_intervals() a parameter choice is specified using an argument called which_pars, not parm.

conf_region for a confidence region for pairs of parameters.

compare_models for an adjusted likelihood ratio test of two models.

adjust_loglik to adjust a user-supplied loglikelihood function.

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)
confint(large)
#> Waiting for profiling to be done...
#>               2.5 %      97.5 %
#> mu[0]    80.3729113 81.96020651
#> mu[1]     2.2435715  3.08287072
#> sigma[0]  3.2992230  4.25997602
#> sigma[1]  0.1611821  0.94337677
#> xi[0]    -0.2740991 -0.11571585
#> xi[1]    -0.1651864 -0.02001604
confint(large, profile = FALSE)
#>               2.5 %      97.5 %
#> mu[0]    80.3790976 81.96120586
#> mu[1]     2.2512784  3.08552167
#> sigma[0]  3.2536759  4.20458854
#> sigma[1]  0.1565970  0.90573598
#> xi[0]    -0.2762342 -0.12166268
#> xi[1]    -0.1593869 -0.01730812