S3 `alogLik`

method to perform loglikelihood adjustment for fitted
extreme value model objects returned from the functions
`gev.fit`

, `gpd.fit`

,
`pp.fit`

and `rlarg.fit`

in the
`ismev`

package. If regression modelling is used then
the model will need to be re-fitted, see `ismev_refits`

.

```
# S3 method for gev.fit
alogLik(x, cluster = NULL, use_vcov = TRUE, ...)
# S3 method for pp.fit
alogLik(x, cluster = NULL, use_vcov = TRUE, ...)
# S3 method for gpd.fit
alogLik(
x,
cluster = NULL,
use_vcov = TRUE,
binom = FALSE,
k,
inc_cens = TRUE,
...
)
# S3 method for rlarg.fit
alogLik(x, cluster = NULL, use_vcov = TRUE, ...)
```

- x
A fitted model object with certain associated S3 methods. See

**Details**.- cluster
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.- use_vcov
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`

).- binom
A logical scalar. This option is only relevant to

**GP models**and is only available in the**stationary**(no covariates) case. If`binom = FALSE`

then loglikelihood adjustment is only performed using the GP model. If`binom = TRUE`

then loglikelihood adjustment is also performed for inferences about the probability of threshold exceedance, using a Bernoulli model for the instances of threshold exceedance.- k
A non-negative integer scalar. This option is only relevant to

**GP models**and is only available in the**stationary**(no covariates) case. If`k`

is supplied then it is passed as the run parameter \(K\) to`kgaps`

for making inferences about the extremal index \(\theta\) using the \(K\)-gaps model of Suveges and Davison (2010).- inc_cens
A logical scalar. This argument is only relevant if

`k`

is supplied. Passed to`kgaps`

to indicate whether or not to include censored inter-exceedance times, relating to the first and last observations.

An object inheriting from class `"chandwich"`

. See

`class(x)`

is a vector of length 5. The first 3 components are

`c("lax", "chandwich", "ismev")`

.
The remaining 2 components depend on the model that was fitted.
The 4th component is:

`"gev"`

if `gev.fit`

(or `gev_refit`

) was used;

`"gpd"`

if `gpd.fit`

(or `gpd_refit`

) was used;

`"pp"`

(or `pp_refit`

) was used;

`"rlarg"`

(or `rlarg_refit`

) was used.
The 5th component is

`"stat"`

if `x$trans = FALSE`

and

`"nonstat"`

if `x$trans = TRUE`

.

See `alogLik`

for details.

If regression modelling is used then the ismev functions
`gev.fit`

, `gpd.fit`

,
`pp.fit`

and `rlarg.fit`

return residuals but `alogLik`

needs the raw data.
The model will need to be re-fitted, using one of the functions in
`ismev_refits`

, and the user will be prompted to do this
by an error message produced by `alogLik`

.

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 ismev package
got_ismev <- requireNamespace("ismev", quietly = TRUE)
if (got_ismev) {
library(ismev)
# GEV model -----
# An example from the ismev::gev.fit documentation
gev_fit <- gev.fit(revdbayes::portpirie, show = FALSE)
adj_gev_fit <- alogLik(gev_fit)
summary(adj_gev_fit)
# An example from chapter 6 of Coles (2001)
data(fremantle)
xdat <- fremantle[, "SeaLevel"]
# Set year 1897 to 1 for consistency with page 113 of Coles (2001)
ydat <- cbind(fremantle[, "Year"] - 1896, fremantle[, "SOI"])
gev_fit <- gev_refit(xdat, ydat, mul = 1:2, show = FALSE)
adj_gev_fit <- alogLik(gev_fit)
summary(adj_gev_fit)
# An example from Chandler and Bate (2007)
gev_fit <- gev_refit(ow$temp, ow, mul = 4, sigl = 4, shl = 4,
show = FALSE)
adj_gev_fit <- alogLik(gev_fit, cluster = ow$year)
summary(adj_gev_fit)
# Get closer to the values reported in Table 2 of Chandler and Bate (2007)
gev_fit <- gev_refit(ow$temp, ow, mul = 4, sigl = 4, shl = 4,
show = FALSE, method = "BFGS")
# Call sandwich::meatCL() with cadjust = FALSE
adj_gev_fit <- alogLik(gev_fit, cluster = ow$year, cadjust = FALSE)
summary(adj_gev_fit)
# GP model -----
# An example from the ismev::gpd.fit documentation
# \donttest{
data(rain)
rain_fit <- gpd.fit(rain, 10, show = FALSE)
adj_rain_fit <- alogLik(rain_fit)
summary(adj_rain_fit)
# Continuing to the regression example on page 119 of Coles (2001)
ydat <- as.matrix((1:length(rain)) / length(rain))
reg_rain_fit <- gpd_refit(rain, 30, ydat = ydat, sigl = 1, siglink = exp,
show = FALSE)
adj_reg_rain_fit <- alogLik(reg_rain_fit)
summary(adj_reg_rain_fit)
# }
# Binomial-GP model -----
# Use Newlyn seas surges data from the exdex package
surges <- exdex::newlyn
u <- quantile(surges, probs = 0.9)
newlyn_fit <- gpd.fit(surges, u, show = FALSE)
# Create 5 clusters each corresponding approximately to 1 year of data
cluster <- rep(1:5, each = 579)[-1]
adj_newlyn_fit <- alogLik(newlyn_fit, cluster = cluster, binom = TRUE,
cadjust = FALSE)
summary(adj_newlyn_fit)
summary(attr(adj_newlyn_fit, "pu_aloglik"))
# Add inference about the extremal index theta, using K = 1
adj_newlyn_theta <- alogLik(newlyn_fit, cluster = cluster, binom = TRUE,
k = 1, cadjust = FALSE)
summary(attr(adj_newlyn_theta, "theta"))
# PP model -----
# An example from the ismev::pp.fit documentation
data(rain)
# Start from the mle to save time
init <- c(40.55755732, 8.99195409, 0.05088103)
muinit <- init[1]
siginit <- init[2]
shinit <- init[3]
rain_fit <- pp_refit(rain, 10, muinit = muinit, siginit = siginit,
shinit = shinit, show = FALSE)
adj_rain_fit <- alogLik(rain_fit)
summary(adj_rain_fit)
# An example from chapter 7 of Coles (2001).
# Code from demo ismev::wooster.temps
data(wooster)
x <- seq(along = wooster)
usin <- function(x, a, b, d) {
return(a + b * sin(((x - d) * 2 * pi) / 365.25))
}
wu <- usin(x, -30, 25, -75)
ydat <- cbind(sin(2 * pi * x / 365.25), cos(2 * pi *x / 365.25))
# Start from the mle to save time
init <- c(-15.3454188, 9.6001844, 28.5493828, 0.5067104, 0.1023488,
0.5129783, -0.3504231)
muinit <- init[1:3]
siginit <- init[4:6]
shinit <- init[7]
wooster.pp <- pp_refit(-wooster, threshold = wu, ydat = ydat, mul = 1:2,
sigl = 1:2, siglink = exp, method = "BFGS",
muinit = muinit, siginit = siginit, shinit = shinit,
show = FALSE)
adj_pp_fit <- alogLik(wooster.pp)
summary(adj_pp_fit)
# r-largest order statistics model -----
# An example based on the ismev::rlarg.fit() documentation
vdata <- revdbayes::venice
rfit <- rlarg.fit(vdata, muinit = 120.54, siginit = 12.78,
shinit = -0.1129, show = FALSE)
adj_rfit <- alogLik(rfit)
summary(adj_rfit)
# \donttest{
# Adapt this example to add a covariate
set.seed(30102019)
ydat <- matrix(runif(nrow(vdata)), nrow(vdata), 1)
rfit2 <- rlarg_refit(vdata, ydat = ydat, mul = 1,
muinit = c(120.54, 0), siginit = 12.78,
shinit = -0.1129, show = FALSE)
adj_rfit2 <- alogLik(rfit2)
summary(adj_rfit2)
# }
}
#> Loading required package: mgcv
#> Loading required package: nlme
#> This is mgcv 1.8-42. For overview type 'help("mgcv-package")'.
#> MLE SE adj. SE
#> loc 122.8000 1.77100 3.17100
#> loc1 -5.5400 2.65700 4.46100
#> scale 12.6500 0.53290 0.81270
#> shape -0.1172 0.01948 0.02741
```