chandwich (version 1.1.6)

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

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.

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.

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)