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.
profile_loglik(
object,
prof_pars = NULL,
prof_vals = NULL,
init = NULL,
type = c("vertical", "cholesky", "spectral", "none"),
...
)An object of class "chandwich" returned by
adjust_loglik.
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.
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.
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.
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.
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.
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.
# -------------------------- 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)
#> [1] -446.0887
#> attr(,"free_pars")
#> mu[0] mu[1] sigma[0] sigma[1] xi[0]
#> 80.9917215 2.5067225 3.7828562 0.2433040 -0.2010567
# 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)
#> [1] -446.3349
#> attr(,"free_pars")
#> mu[0] mu[1] sigma[0] sigma[1]
#> 80.9256806 2.4551350 3.6790743 0.4088366