Profile-likelihood inference for TVGEV
objects.
# S3 method for TVGEV
profLik(object, fun, level = 0.95, deriv = TRUE, trace = 0, ...)
A TVGEV
object.
A function of the parameter vector for which the
profile-likelihood will be carried over. This function must have
the arguments: psi
for the vector of parameters and
object
for the model object; so the function can use the
some of slots of object
. If needed, a wrapper function can
be used use more arguments, see Details.
Level of confidence. Can be of length > 1
.
Logical. If TRUE
, the function fun
is
assumed to provide a gradient vector as an attribute named
"gradient"
of the result. For now deriv
can only be
TRUE
, which implies that fun
must compute the
gradient.
Level of verbosity; trace = 0
prints nothing.
Not used yet.
An array with the value of the function and the
corresponding Lower and Upper end-points for the given confidence
levels. This array has two attributes with names "diagno"
and "psi"
which both are arrays. The attributes provide
information about the numerical optimisation and the values of the
vector of parameter that maximised or minimised the function
fun
.
Compute the lower and upper end-points of a profile-likelihood based confidence interval. The (apparently new) method used here relies on maximising and minimising the function of interest, say \(\eta(\boldsymbol(\psi)\), under the constraint that the log-likelihood is greater than the maximal log-likelihood minus a positive quantity \(\delta\) depending on the confidence level. This differs from the usual method which relies on an univariate zero-finding for the profile-likelihood function (minus a constant). Remind that each evaluation of the profile requires a \(p-1\) dimensional optimisation. As a major advantage, the new method does not require a re-parameterisation of the model.
For each confidence limit the numerical optimisation may
fail, in which case the limit will be NA
. Using trace
= 1
can be useful to further check the optimisation. The
Optimisation status
value should be 3
or 4
for small changes on the parameter or on the objective. On the
other hand, a value of 5
indicates that the maximal number
of iterations was reached, which is considered here as a
failure. The Constaint check
value should be small because
the constraint must be active at the optimum. The gradDist
is the distance between the two directions of the gradient vectors
(objective and constraint). It should be small as well because the
gradients must be colinear at the optimum (Lagrange conditions).
Deville Y. (2017) "Profile-likelihood using constrained optimisation". Unpublished Tech. Report.
df <- within(TXMax_Dijon, Date <- as.Date(sprintf("%4d-01-01", Year)))
## fit a TVGEV model with constant parameters.
res1 <- TVGEV(data = df, response = "TXMax", date = "Date",
estim = "nloptr")
## define a function of the parameter vector: here the first component
## This is for illustration only since the the result can be obtained
## using the 'confint' method with \code{method = "proflik"}, which
## gives the confidence intervals for each of the parameters.
myfun <- function(psi, object) {
res <- psi[1]
grad <- rep(0.0, object$p)
grad[1] <- 1
attr(res, "gradient") <- grad
res
}
pl <- profLik(object = res1, fun = myfun, deriv = TRUE)
confint(res1, method = "proflik")
#>
#>
#> o Finding CI for "mu_0"
#>
#> o 95%, lower bound: 32.51
#> Constraint check -0.0000000, 184.2918
#> 95%, upper bound: 33.39
#> Constraint & value -0.0000001, 184.2918
#>
#>
#> o Finding CI for "sigma_0"
#>
#> o 95%, lower bound: 1.61
#> Constraint check -0.0000000, 184.2918
#> 95%, upper bound: 2.24
#> Constraint & value -0.0000002, 184.2918
#>
#>
#> o Finding CI for "xi_0"
#>
#> o 95%, lower bound: -0.31
#> Constraint check -0.0000001, 184.2918
#> 95%, upper bound: -0.04
#> Constraint & value -0.0000000, 184.2918
#> , , 95%
#>
#> L U
#> mu_0 32.507 33.387
#> sigma_0 1.611 2.236
#> xi_0 -0.313 -0.044
#>