S3 alogLik
method to perform loglikelihood adjustment for fitted
extreme value model objects returned from the functions
gevrFit
and gpdFit
in the eva package.
A fitted model object with certain associated S3 methods. See Details.
A vector or factor indicating from which cluster the
respective log-likelihood contributions from loglik
originate.
The length of cluster
must be consistent with the estfun
method to be used in the estimation of the 'meat' V
of the sandwich
estimator of the covariance matrix of the parameters to be passed to
adjust_loglik
. In most cases, cluster
must have length equal to the number of observations in data. The
exception is the GP (only) model (binom = FALSE
), where the
cluster
may either contain a value for each observation in the raw
data, or for each threshold exceedance in the data.
If cluster
is not supplied (is NULL
) then it is
assumed that each observation forms its own cluster.
See Details for further details.
A logical scalar. Should we use the vcov
S3 method
for x
(if this exists) to estimate the Hessian of the independence
loglikelihood to be passed as the argument H
to
adjust_loglik
?
Otherwise, H
is estimated inside
adjust_loglik
using
optimHess
.
Further arguments to be passed to the functions in the
sandwich package meat
(if cluster = NULL
),
or meatCL
(if cluster
is not
NULL
).
An object inheriting from class "chandwich"
. See
class(x)
is a vector of length 5. The first 3 components are
c("lax", "chandwich", "eva")
.
The 4th component depends on which model was fitted.
"rlarg"
if gevrFit
was used;
"gpd"
if gpdFit
was used.
The 5th component is
"stat"
if there are no covariates in the mode and
"nonstat"
otherwise.
See alogLik
for details.
In the stationary case (no covariates) the function
gevrFit
and gpdFit
in the eva
package offer standard errors based on the expected information or on the
observed information, via the argument information
. In contrast,
alogLik()
always bases calculations on the observed information
matrix. Therefore, unadjusted standard errors resulting from
alogLik()
may be different the corresponding standard errors
from gevrFit
or gpdFit
.
For gevrFit
only GEV fits (gumbel = FALSE
) are
supported.
Chandler, R. E. and Bate, S. (2007). Inference for clustered data using the independence loglikelihood. Biometrika, 94(1), 167-183. doi:10.1093/biomet/asm015
Suveges, M. and Davison, A. C. (2010) Model misspecification in peaks over threshold analysis, The Annals of Applied Statistics, 4(1), 203-221. doi:10.1214/09-AOAS292
Zeileis (2006) Object-Oriented Computation and Sandwich Estimators. Journal of Statistical Software, 16, 1-16. doi:10.18637/jss.v016.i09
alogLik
: loglikelihood adjustment for model fits.
# We need the eva package
got_eva <- requireNamespace("eva", quietly = TRUE)
if (got_eva) {
library(eva)
# An example from the eva::gpdFit documentation
set.seed(7)
x <- eva::rgpd(2000, loc = 0, scale = 2, shape = 0.2)
mle_fit <- eva::gpdFit(x, threshold = 4, method = "mle")
adj_mle_fit <- alogLik(mle_fit)
summary(adj_mle_fit)
# Another example from the eva::gpdFit documentation
# A linear trend in the scale parameter
set.seed(7)
n <- 300
x2 <- eva::rgpd(n, loc = 0, scale = 1 + 1:n / 200, shape = 0)
covs <- as.data.frame(seq(1, n, 1))
names(covs) <- c("Trend1")
result1 <- eva::gpdFit(x2, threshold = 0, scalevars = covs,
scaleform = ~ Trend1)
adj_result1 <- alogLik(result1)
summary(adj_result1)
# An example from the eva::gevrFit documentation
set.seed(7)
x1 <- eva::rgevr(500, 1, loc = 0.5, scale = 1, shape = 0.3)
result1 <- eva::gevrFit(x1, method = "mle")
adj_result1 <- alogLik(result1)
summary(adj_result1)
# Another example from the eva::gevrFit documentation
# A linear trend in the location and scale parameter
n <- 100
r <- 10
x2 <- eva::rgevr(n, r, loc = 100 + 1:n / 50, scale = 1 + 1:n / 300,
shape = 0)
covs <- as.data.frame(seq(1, n, 1))
names(covs) <- c("Trend1")
# Create some unrelated covariates
covs$Trend2 <- rnorm(n)
covs$Trend3 <- 30 * runif(n)
result2 <- eva::gevrFit(data = x2, method = "mle", locvars = covs,
locform = ~ Trend1 + Trend2*Trend3,
scalevars = covs, scaleform = ~ Trend1)
adj_result2 <- alogLik(result2)
summary(adj_result2)
}
#>
#> Attaching package: 'eva'
#> The following objects are masked from 'package:evd':
#>
#> dgpd, pgev, pgpd, qgev, qgpd, rgpd
#> MLE SE adj. SE
#> Location (Intercept) 100.200000 0.199200 0.229000
#> Location Trend1 0.018110 0.003192 0.003304
#> Location Trend2 -0.004186 0.063170 0.066540
#> Location Trend3 -0.003237 0.003934 0.003736
#> Location Trend2:Trend3 -0.001259 0.004247 0.004489
#> Scale (Intercept) 1.123000 0.092090 0.115700
#> Scale Trend1 0.002235 0.001208 0.001218
#> Shape (Intercept) 0.034370 0.030470 0.030090