profile_loglik: Profile loglikelihood

In chandwich: Chandler-Bate Sandwich Loglikelihood Adjustment

Profile loglikelihood

Description

Calculates the profile loglikelihood for a subset of the model parameters. This function is provided primarily so that it can be called by conf_intervals and conf_region.

Usage

profile_loglik(
  object,
  prof_pars = NULL,
  prof_vals = NULL,
  init = NULL,
  type = c("vertical", "cholesky", "spectral", "none"),
  ...
)

Arguments

object	An object of class "chandwich" returned by adjust_loglik.
prof_pars	A vector specifying the subset of the (unfixed) parameters over which to profile. 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. prof_pars must not have any parameters in common with attr(object, "fixed_pars"). prof_pars must not contain all of the unfixed parameters, i.e. there is no point in profiling over all of the unfixed parameters.
prof_vals	A numeric vector. Values of the parameters in prof_pars. If prof_vals = NULL then the MLEs stored in object of the parameters prof_pars are used.
init	A numeric vector of initial estimates of the values of the parameters that are not fixed and are not in prof_pars. Should have length attr(object, "p_current") - length(prof_pars). If init is NULL or is of the wrong length then the relevant components from the MLE stored in object are used.
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

A numeric vector of length 1. The value of the profile loglikelihood. The returned object has the attribute "free_pars" giving the optimal values of the parameters that remain after the parameters prof_pars and attr(object, "fixed_pars") have been removed from the full parameter vector. If there are no such parameters, which happens if an attempt is made to profile over all non-fixed parameters, then this attribute is not present and the value returned is calculated using the function object.

See Also

adjust_loglik to adjust a user-supplied loglikelihood function.

conf_intervals for confidence intervals for individual parameters.

conf_region for a confidence region for a pair of parameters.

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)

# Profile loglikelihood for xi1, evaluated at xi1 = 0
profile_loglik(large, prof_pars = "xi[1]", prof_vals = 0)

# Model with xi1 fixed at 0
medium <- adjust_loglik(larger = large, fixed_pars = "xi[1]")
# Profile loglikelihood for xi0, evaluated at xi0 = -0.1
profile_loglik(medium, prof_pars = "xi[0]", prof_vals = -0.1)

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