---
title: "Local influence diagnostics for the EVBS regression model"
author: "Raydonal Ospina"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Local influence diagnostics for the EVBS regression model}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)
```

## Introduction

The **evbsreg** package implements local influence diagnostics for the
Extreme-Value Birnbaum--Saunders (EVBS) regression model. This vignette
walks through the complete workflow: fitting the model, diagnosing
influential observations, assessing model adequacy, and interpreting the
results, using the bundled Itajaí wind gust dataset. It then covers the
quantities that follow from a fitted model (finite upper endpoint, return
levels and expected shortfall), the prospective monitoring index, and the
same diagnostics under a generalized extreme-value (GEV) regression.

Every result below is reproducible: the computations are deterministic except
the simulated envelope and the block bootstrap, which are run after
`set.seed()`.

```{r setup}
library(evbsreg)
```

## The data

The `itajai` dataset contains 124 monthly maximum wind gust speeds and the
corresponding daily mean atmospheric pressure, recorded at INMET station
A-868 in Itajaí, Brazil, from July 2010 to October 2020.

```{r data}
data(itajai)
str(itajai)
summary(itajai$wind)
```

Observation 82 is the catastrophic event of 26 April 2017 (33.9 m/s).

## Fitting the model

The design matrix must include an intercept column. We model the location
of the log-EVBS distribution as a linear function of atmospheric pressure.

```{r fit}
X <- cbind(1, itajai$pressure)
fit <- evbsreg.fit(X, itajai$wind)

data.frame(
  Parameter = c("beta0", "beta1", "alpha", "gamma"),
  Estimate  = round(fit$coeff, 4),
  SE        = round(c(fit$stderrors, fit$stderroralpha, fit$stderrorgama), 4)
)
```

Both regression coefficients are highly significant:

```{r pvalues}
round(fit$pvalues, 4)
```

## Local influence diagnostics

The `cnc_diagnostics()` function computes the conformal normal curvature
diagnostics from the fitted object.

```{r diag}
diag <- cnc_diagnostics(fit)

## Top four normalized eigenvalues
round(head(diag$eigenvalues_norm, 4), 5)

## Observations flagged at q = 7
which(diag$Bj[7, ] > diag$bq[7])
```

The two-panel diagnostic figure shows the normalized eigenvalues (left) and
the aggregate contributions (right). Observation 82 dominates.

```{r plot-cnc, fig.width = 10, fig.height = 4.5}
plot_cnc(diag, q = 7)
```

## Deletion analysis

We refit the model without the flagged observation and measure the
relative change in each parameter.

```{r deletion}
fit82 <- evbsreg.fit(X[-82, ], itajai$wind[-82])
rc <- 100 * (fit82$coeff - fit$coeff) / abs(fit$coeff)
names(rc) <- c("beta0", "beta1", "alpha", "gamma")
round(rc, 2)
```

The tail-shape parameter $\gamma$ changes by `r round(rc[4], 1)`%, while the
regression coefficients and the shape parameter $\alpha$ change by at most
`r round(max(abs(rc[1:3])), 1)`%.
The influence is therefore concentrated almost entirely on the tail.

## Model adequacy

Randomized quantile residuals should be approximately standard normal under
a correct specification. For a continuous response they are deterministic,
$r_i = \Phi^{-1}\{\hat F(t_i)\}$.

```{r residuals}
r <- rqrandomized(X, itajai$wind)
shapiro.test(r)
```

The simulated envelope is random, so fix the seed to reproduce it.

```{r envelope, fig.width = 6, fig.height = 6}
set.seed(2023)
envelope_qq(X, itajai$wind, nrep = 100)
```

## Endpoint, return levels and expected shortfall

When $\hat\gamma < 0$ the fitted EVBS distribution has a finite upper
endpoint, which depends on the covariate through the linear predictor
$\eta = x^\top\beta$. At the covariate value of observation 82:

```{r endpoint}
eta82 <- sum(X[82, ] * fit$betahat)
c(endpoint = unname(evbs_endpoint(eta82, fit$alphahat, fit$gamahat)),
  observed = itajai$wind[82])
```

Return levels come from the exact quantile function `qevbs()`. With monthly
maxima, `period` is in months, so `12 * 50` is the level exceeded on average
once every 50 years (at the given covariate value, here the mean pressure):

```{r return-level}
evbs_return_level(fit, x = c(1, mean(itajai$pressure)),
                  period = 12 * 50, es = TRUE)
## 10-, 50- and 100-year levels at the same covariate value
qevbs(1 - 1 / (12 * c(10, 50, 100)),
      eta = sum(c(1, mean(itajai$pressure)) * fit$betahat),
      alpha = fit$alphahat, gama = fit$gamahat)
```

Note that $\exp(x^\top\hat\beta)$ is a scale, not the conditional median;
conditional quantiles always come from `qevbs()`.

## Standard errors under serial dependence

If the residuals show serial dependence, moving-block bootstrap standard
errors are a safer guide than the information-based ones:

```{r blockboot}
set.seed(2023)
bb <- evbs_block_boot(X, itajai$wind, L = 4, B = 100)
round(bb$se, 4)
```

## Prospective monitoring

`evbs_monitor()` refits the model at each instant after a burn-in, finds the
observation flagged by the curvature diagnostic, and compares it with the
fitted upper endpoint. Values of the index $F_t$ below one mean that the
flagged observation is close to or beyond what the model can accommodate.

```{r monitor, fig.width = 7, fig.height = 4.5}
ct <- evbs_monitor(X, itajai$wind, burn = 40)
ct[ct$F < 1, c("t", "flagged", "gama", "endpoint", "F")][1:4, ]
plot(ct)
```

Apart from a single dip at $t = 43$, the index stays above one until the
April 2017 event enters the sample at $t = 82$, and stays below one from then
on. Setting `use = "max"` gives the naive version that always deletes the
sample maximum instead of the observation flagged by the curvature.

## Comparison with GEV regression

The same diagnostics are available for a GEV regression fitted to the
response scale. Centre the covariates first: the GEV likelihood is poorly
conditioned otherwise.

```{r gev}
Xc <- cbind(1, itajai$pressure - mean(itajai$pressure))
gfit <- gevreg.fit(Xc, itajai$wind)
gfit$convergence
gdiag <- cnc_diagnostics_gev(gfit)
```

To compare influence rankings, use the same statistic for both families:
`colMeans(diag$Bj)` for the EVBS fit and `Bbar` for the GEV fit.

```{r gev-compare}
Bbar_evbs <- colMeans(diag$Bj)
head(order(Bbar_evbs, decreasing = TRUE), 3)
head(order(gdiag$Bbar, decreasing = TRUE), 3)
cor(Bbar_evbs, gdiag$Bbar, method = "spearman")
```

The prospective index under the GEV model is `gev_monitor(Xc, itajai$wind)`.

## Density shapes

The package also provides the density-plotting functions used for the
density figures of the paper.

```{r densities, fig.width = 7, fig.height = 5}
plot_evbs_alpha()
```

## Reproducing the full study

Five standalone scripts reproduce the density figures, the Itajaí
application and the Monte Carlo study. They write only to `tempdir()`. After
installation they are available via `system.file()`:

```{r scripts, eval = FALSE}
# Density figures
source(system.file("scripts/script_01_density_figures.R", package = "evbsreg"))

# Itajai application
source(system.file("scripts/script_02_itajai_application.R", package = "evbsreg"))

# Monte Carlo study; set m <- 500 inside for a quick check
source(system.file("scripts/script_03_simulation_scenario1.R", package = "evbsreg"))
source(system.file("scripts/script_04_simulation_scenario2.R", package = "evbsreg"))
source(system.file("scripts/script_05_simulation_scenario3.R", package = "evbsreg"))
```

## References

Ferreira, M., Gomes, M. I., and Leiva, V. (2012). On an extreme value
version of the Birnbaum--Saunders distribution. *REVSTAT*, 10, 181--210.

Ospina, R., Lima, J. I. C., Barros, M., and Macêdo, A. M. S. (2026).
*Local influence diagnostics for extreme-value Birnbaum--Saunders
regression, with an application to wind-gust extremes.* Submitted.