Covariate-driven extremes with margEVT

This vignette walks through a complete margEVT analysis: penalized estimation of a non-homogeneous Poisson point process (NHPP) for threshold exceedances, model diagnostics, marginal (covariate-integrated) return levels, bootstrap uncertainty, and walk-forward validation.

Model

Exceedances of a threshold $u$ by a series $Y_t$ are modelled as a point process with intensity \[ \Lambda_t\{(y, \infty)\} = \frac{1}{n_y} \left[1 + \xi_t \frac{y - \mu_t}{\sigma_t}\right]_+^{-1/\xi_t}, \qquad y > u, \] where $n_y$ is the number of observations per year, so that $\mu_t$, $\sigma_t$ and $\xi_t$ are the parameters of the annual-maximum GEV distribution. Each parameter is linear in covariates ($\log\sigma_t$ for the scale), and the covariate coefficients are estimated by penalized maximum likelihood with an elastic-net penalty weighted coordinate-wise; with penalty_weights = "adaptive" the weights are $1/|\tilde\theta_k|^\gamma$, the adaptive LASSO.

Since $G(z \mid \mathbf{v}) = \exp\{-\sum_t \Lambda_t(z)\}$ is the annual non-exceedance probability given one year of covariates $\mathbf{v}$, the marginal return level for period $T$ solves \[ G_\Pi(z_T) = \mathrm{E}_\Pi\{G(z_T \mid \mathbf{V})\} = 1 - 1/T , \] where $\Pi$ is a law for annual covariate trajectories: fixed scenarios (approach A), a fitted VAR generator (approach B), or resampled observed years (approach C).

Simulated data

We simulate 40 years of weekly data with an annual cycle, an autocorrelated climate index x that drives the location, and an irrelevant covariate z.

set.seed(2024)
n_years <- 40; n_y <- 52
n <- n_years * n_y
week <- rep(seq_len(n_y), n_years)
df <- cbind(
  data.frame(year = rep(1981:2020, each = n_y),
             x = as.numeric(arima.sim(list(ar = 0.8), n)) * 0.6,
             z = rnorm(n)),
  seasonal_harmonics(week, period = n_y)
)
df$y <- 20 + 3 * df$x + 4 * df$cos1 + 3 * rexp(n)
u <- quantile(df$y, 0.9)

seasonal_harmonics() builds the cos1, sen1, … columns with the same convention used later when annual covariate blocks are rebuilt.

Estimation

Harmonics are declared in free_vars, so they are never penalized. The penalty level is selected by BIC along a regularization path, and the reported model is the path solution at the selected penalty.

fit <- fit_nhpp(df, threshold = u,
                loc_vars   = c("cos1", "sen1", "x", "z"),
                scale_vars = c("x", "z"),
                free_vars  = c("cos1", "sen1"),
                penalty = "lasso", lambda = "bic",
                penalty_weights = "adaptive",
                obs_per_year = n_y, verbose = FALSE)
fit
#> -- nhpp_fit --------------------------------------
#>   Threshold    : 29.11
#>   Penalty      : lasso
#>   Lambda       : mu=2.9444, sigma=2.9444, xi=2.9444
#>   Weights      : adaptive: mle
#>   Alpha        : 1
#>   Converged    : TRUE
#>   nllh (raw)   : 162.3987
#>   nllh (pen)   : 165.7571
#>   obs/year     : 52.00
#>   Active params: 6 of 9
#>   Coefficients (non-zero):
#>    mu.(Intercept)           mu.cos1           mu.sen1              mu.x 
#>          31.60851           4.04716           0.33637           2.48863 
#> sigma.(Intercept)    xi.(Intercept) 
#>           1.01718          -0.03196

summary(fit)
#> -- nhpp_fit summary --------------------------------------
#>   Threshold      : 29.11
#>   Penalty        : lasso
#>   Lambda         : mu=2.9444, sigma=2.9444, xi=2.9444
#>   Weights        : adaptive: mle
#>   Alpha          : 1
#>   Penalize shape : TRUE
#>   obs/year       : 52.00
#>   Converged      : TRUE
#>   KKT violation  : 8.46e-08 (per exceedance)
#>   nllh (raw)     : 162.3987
#>   nllh (pen)     : 165.7571
#>   Exceedances    : 208 | BIC: 356.888
#> 
#>   Active coefficients (6 of 9):
#>                   Estimate
#> mu.(Intercept)    31.60851
#> mu.cos1            4.04716
#> mu.sen1            0.33637
#> mu.x               2.48863
#> sigma.(Intercept)  1.01718
#> xi.(Intercept)    -0.03196
#> 
#>   Shrunk to zero : mu.z, sigma.x, sigma.z
active_covariates(fit)
#> [1] "x"

Convergence is certified by the Karush-Kuhn-Tucker conditions of the penalized problem (fit$kkt, the largest violation per exceedance). The regularization path and BIC curve are stored in fit$path:

head(fit$path$table)
#>      lambda  phase      bic k_active converged selected
#> 1 2000.0000 coarse 485.3724        5      TRUE    FALSE
#> 2 1524.1526 coarse 485.3724        5      TRUE    FALSE
#> 3 1161.5206 coarse 485.3724        5      TRUE    FALSE
#> 4  885.1673 coarse 485.3724        5      TRUE    FALSE
#> 5  674.5651 coarse 485.3724        5      TRUE    FALSE
#> 6  514.0701 coarse 485.3724        5      TRUE    FALSE
plot(fit, type = "path")

Diagnostics

The time-rescaling theorem gives interarrival residuals that are standard exponential under a correct model; magnitude residuals check the fitted generalized Pareto tail. Q-Q plots show exact pointwise 95% bands.

op <- par(mfrow = c(1, 2))
plot(fit, type = "qq", residuals = "interarrival")
plot(fit, type = "qq", residuals = "magnitude")

par(op)

Standard generics are available: coef(), logLik(), AIC(), BIC(), nobs(), predict(), and residuals().

logLik(fit)
#> 'log Lik.' -162.4315 (df=6)
predict(fit, newdata = data.frame(cos1 = 1, sen1 = 0, x = c(-1, 0, 1),
                                  z = 0))
#> $mu
#> [1] 33.16704 35.65567 38.14430
#> 
#> $sigma
#> [1] 2.765376 2.765376 2.765376
#> 
#> $xi
#> [1] -0.03196123 -0.03196123 -0.03196123

Marginal return levels

Approach B needs a model for whole years of covariates. The generator removes the seasonal cycle, fits a VAR(p) with order selected by the Schwarz criterion, and simulates stationary trajectories.

gen <- fit_var_generator(fit, df, vars = "x", lag_max = 4)
gen$p_opt
#> [1] 1
sims <- simulate_covariates(gen, n_mc = 300, seed = 1)

Approach A evaluates fixed scenarios, B integrates over the simulated trajectories, and C resamples observed years. Monte Carlo standard errors of the roots are reported for B and C.

rl <- marginalize(fit, df, TRs = c(2, 10, 50, 100),
                  approaches = c("A", "B", "C"),
                  scenarios = list(neutral = list(x = 0),
                                   positive = list(x = 1)),
                  mc_sample = sims, n_boot = 300, seed = 1)
rl_table(rl)
#>   approach            scenario      T2     T10     T50    T100
#> 1        A             neutral 33.9255 38.9870 43.1849 44.8952
#> 2        A            positive 36.4142 41.4756 45.6735 47.3839
#> 3        B       parametric_mc 34.7541 40.0346 44.3141 46.0518
#> 4        C empirical_bootstrap 34.8285 40.0566 44.3210 46.0551
plot(rl)

Scenario-based levels can be far from the marginal ones: the marginal law averages over within-year covariate variability, which a single fixed value does not represent.

Uncertainty

bootstrap_rl() refits the selected specification (same penalty level and weights) on parametric bootstrap samples and recomputes the return levels. Each replicate uses its own L’Ecuyer-CMRG random-number stream, so results are identical for any number of cores.

ci <- bootstrap_rl(fit, df, TRs = c(10, 50), R = 20, approach = "C",
                   marginalize_args = list(n_boot = 100, seed = 1),
                   verbose = FALSE)
ci
#>   TR  RL_est  CI_low CI_high n_ok
#> 1 10 40.1088 38.4451 41.6862   20
#> 2 50 44.3752 42.2021 47.5903   20

The intervals are conditional on the selected model and do not include model-selection uncertainty. In applications, use several hundred replicates and cores > 1.

Walk-forward validation

backtest() refits the model on expanding training windows, computes return levels from the training years, and compares the observed annual maxima of the following years with the nominal exceedance rates. Adaptive weights are recomputed inside each window, so no validation data enters the training fit; reselect_lambda = TRUE also re-selects the penalty, and trend_col handles a linear time trend.

bt <- backtest(fit, df, varname = "y", TRs = c(2, 5),
               window_years = 10, min_train_years = 20, n_boot = 100,
               verbose = FALSE)
bt$binom_tests
#>   TR expected observed n_exceedances n_total p_value not_rejected calibrated
#> 1  2      0.5     0.45             9      20  0.8238         TRUE       TRUE
#> 2  5      0.2     0.10             2      20  0.4019         TRUE       TRUE