In RuoqingYin/oolax: Object-orientated Loglikelihood Adjustment for Extreme Value Models

options(tinytex.verbose = TRUE)
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

The oolax package performs loglikelihood adjustments of extreme value models with clustered data. It uses sandwich estimator of the parameter covariate matrix, based on the methodology in @CB2007. This adjustment is built from (or, say, limited to) Generalised Extreme Value model (GEV) and generalised Pareto distribution (GPD) from packages ismev and evd.

Loglikelihood adjustment using the alogLik function

The alogLik function returns an object of class "oolax", "chandwich". It can be used to adjust the standard error of the estimated parameters in extreme value models given clustered data. This function is particular designed for fitted object from evd and ismev package.

It is an S3 method which uses coef, vcov, and nobs method for the fitted extreme value objects which have class gev or pot if they are from evd package, or class ismev_gev or ismev_gpd if they are from ismev package. If some of the above methods are not available, it will be created in oolax. Some work in other functions in the package is to link the methods to classes, and make sure in alogLik function, it would check if the object is from the correct class, and suitable methods can be applied to that class, then use the correct information to do the adjustment based on chandwich::adjust_loglik. Since we use chandwich::adjust_loglik to perform the loglikelihood adjustment, we created a function called logLikVec which provides a vector of the contributions to the independence loglikelihood from individual observations for GEV or GP distribution. The default is cluster=NULL then each observation is assumed to be in its own cluster.

We illustrate the loglikelihood adjustments and the use of the function in oolax using one example, starting from a GEV model for annual maximum temperatures from two geographical locations. The clustering arises because a given pair of observations from these two locations occur in the same year. Therefore, year of observation is the cluster variable.

Extreme value modelling of maximum temperatures

We consider the example shown in Section 5.2 of @CB2007. The owtemps data (which can be found in chandwich package) contains annual maximum temperatures in Oxford and Worthing in the U.K. from 1901 to 1980. Year is considered as cluster, therefore there are 80 clusters of independent observations from a bivariate distribution.

Model fitting developed from ismev

The loglikelihood adjustment is based on generalized extreme value (GEV) models which fitted by oogev.fit. This function is a slightly modified version of gev.fit from ismev.For a non-stationary fitted, we need some additional information which has not created in the original gev.fit nor gpd.fit to perform loglikelihood adjustment, therefore we modified those and create oogev.fit and oogpd.fit.

The model is parameterized so that the marginal distribution at Oxford is GEV($\mu_0 + \mu_1, \sigma_0 + \sigma_1, \xi_0 + \xi_1$) and at Worthing is GEV($\mu_0 - \mu_1, \sigma_0 - \sigma_1, \xi_0 - \xi_1$), where GEV($\mu, \sigma, \xi$) denotes a GEV distribution with location $\mu$, scale $\sigma$ and shape $\xi$.

We perform loglikelihood adjustment for the full six-parameter model and reproduce the relevant rows of Table 2 in @CB2007.

library(oolax)
y <- c(chandwich::owtemps[, "Oxford"], chandwich::owtemps[, "Worthing"])
x <- as.matrix(rep(c(1, -1), each = length(y) / 2))
# Fit the model by oogev.fit which allows non-stationary parameters
owfit <- oogev.fit(y, x, mul = 1, sigl = 1, shl = 1, method = "BFGS" )
year <- rep(1:(length(y) / 2), 2)
# Perform the loglikelihood adjustment of the full model
adj_owfit <- alogLik(owfit, cluster = year, cadjust = FALSE)
# Provides MLE, standard error of MLE, and adjusted standard error of MLE.
summary(adj_owfit)

Confidence intervals

We use confint to calculate confidence intervals for one or more model parameters based on an object returned from alogLik.

library(chandwich)
# 95% confidence intervals, vertically adjusted
conf_intervals(adj_owfit)

Confidence regions

We compute the confidence region for a set of parameter ($\sigma, \sigma_1$), and reproduce partial of Figure 4(b) of @CB2007, adding a 95\% confidence region for the vertical adjustment.

which_pars <- c("scale", "scale1")
gev_none <- conf_region(adj_owfit, which_pars = which_pars, type = "none")
gev_vertical <- conf_region(adj_owfit, which_pars = which_pars)
plot(gev_none, gev_vertical, lwd = 2, xlim = c(3.0, 4.5), ylim = c(-0.1, 1.25))

The 95\% contours of the profile adjusted loglikelihoods are similar.

Comparing nested models

Suppose that we want to test the null hypothesis that both Oxford and Worthing share a common GEV shape parameter, that is, $shape1 = 0$. The nested model can be fitted by oogev.fit, and we perform loglikelihood adjustment under it. Then anova is used to conduct an adjusted likelihood ratio test.

owfit_small <- oogev.fit(y, x, mul = 1, sigl = 1, method = "BFGS" )
adj_owfit_small <- alogLik(owfit_small, cluster = year, cadjust = FALSE)
anova(adj_owfit_small, adj_owfit)

Model fitting developed from evd

This time, the loglikelihood adjustment is based on generalized extreme value (GEV) models which fitted by fgev, which is originally from evd. The model is parameterized so that the marginal distribution at Oxford is GEV($\mu_0 + \mu_1, \sigma_0, \xi_0$) and at Worthing is GEV($\mu_0 - \mu_1, \sigma_0, \xi_0$), where GEV($\mu, \sigma, \xi$) denotes a GEV distribution with location $\mu$, scale $\sigma$ and shape $\xi$.

y <- c(chandwich::owtemps[, "Oxford"], chandwich::owtemps[, "Worthing"])
x <- rep(c(-1, 1), each = length(y) / 2)
owfit <- evd::fgev(y, nsloc = x)
year <- rep(1:(length(y) / 2), 2)
oola_small <- alogLik(owfit, cluster = year)
summary(oola_small)

Confidence intervals

We use confint to calculate confidence intervals for one or more model parameters based on an object returned from alogLik.

library(chandwich)
# 95% confidence intervals, vertically adjusted
conf_intervals(oola_small)

Confidence regions

We compute the confidence region for a set of parameter ($\mu_0, \mu_1$), and adding 95\% confidence region for the vertical adjustment.

which_pars <- c("loc", "loctrend")
gev_none <- conf_region(oola_small, which_pars = which_pars, type = "none")
gev_vertical <- conf_region(oola_small, which_pars = which_pars)
plot(gev_none, gev_vertical, xlim = c(80, 82.3), ylim = c(-3.5, -1.5))

The 95\% contours of the profile adjusted loglikelihoods are similar.

Comparing nested models

Suppose that we want to test the null hypothesis that both Oxford and Worthing share a common GEV location parameter, that is, $loctrend = 0$. The nested model can be fitted by fgev, and we perform loglikelihood adjustment under it. Then anova is used to conduct an adjusted likelihood ratio test.

owfit <- evd::fgev(y)
year <- rep(1:(length(y) / 2), 2)
oola_tiny <- alogLik(owfit, cluster = year)
anova(oola_small, oola_tiny)

RuoqingYin/oolax documentation built on May 28, 2019, 12:20 p.m.