Package 'evbsreg'

Conformal Normal Curvature Local Influence Diagnostics

Description

Computes conformal normal curvature (CNC) local influence diagnostics from a fitted EVBS regression model. The CNC approach of Poon and Poon (1999) produces a scale-invariant influence measure bounded in $[0,1]$; the aggregate contribution statistic of Zhu and Lee (2001) attributes curvature to individual observations and is compared against interpretable reference thresholds.

Usage

cnc_diagnostics(fit)

Arguments

fit	A fitted model object returned by evbsreg.fit. The function uses the influence matrix fit$B.

Details

The function performs the symmetric eigendecomposition of the influence matrix $B$, normalizes the eigenvalues to unit norm, and for each threshold $q = 1, \ldots, 7$ identifies the $q$-influential eigenvectors (those whose normalized eigenvalue exceeds $q/\sqrt{n}$). The aggregate contribution of observation $j$ at level $q$, $B_j(q)$, is the sum of the normalized eigenvalues weighted by the squared eigenvector coordinates of observation $j$.

Value

A list with components:

eigenvalues

Numeric vector: the raw eigenvalues of $B$.

eigenvalues_norm

Numeric vector: the absolute normalized eigenvalues $|\lambda_i^*|$.

eigenvectors

Matrix: the eigenvectors of $B$ (columns).

thresholds

Numeric vector of length 7: the reference thresholds $q/\sqrt{n}$ for $q = 1, \ldots, 7$.

Bj

Matrix of dimension $8 \times n$: rows 1–7 hold the aggregate contributions $B_j(q)$ for $q=1,\ldots,7$; row 8 holds the global aggregate over all eigenvectors.

bq

Numeric vector of length 8: the reference values $b(q)$ for flagging influential observations at each level.

n

Integer: the sample size.

References

Poon, W.-Y. and Poon, Y. S. (1999). Conformal normal curvature and assessment of local influence. Journal of the Royal Statistical Society, Series B, 61, 51–61.

Zhu, H. and Lee, S. (2001). Local influence for incomplete-data models. Journal of the Royal Statistical Society, Series B, 63, 111–126.

Ospina, R., Lima, J. I. C., Barros, M., and Macedo, A. M. S. (2026). Local influence diagnostics for the extreme-value Birnbaum-Saunders regression model. Submitted.

See Also

evbsreg.fit, plot_cnc.

Examples

data(itajai)
X <- cbind(1, itajai$pressure)
fit <- evbsreg.fit(X, itajai$wind)
diag <- cnc_diagnostics(fit)

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

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

---

Conformal normal curvature diagnostics for the GEV regression model

Description

Local influence diagnostics under case-weight perturbation for the GEV regression model, obtained from the individual score contributions (see gev_scores) and the observed information.

Usage

cnc_diagnostics_gev(object)

Arguments

object

An object of class gevreg from gevreg.fit.

Value

A list with Bj, the diagonal of the influence matrix normalised to sum to one, the threshold 2/n, the indices of the observations with Bj above it, and Bbar, the mean aggregate contribution over the thresholds $q = 1, \ldots, 7$ and the unthresholded contribution, computed exactly as colMeans(cnc_diagnostics(fit)$Bj) is for the EVBS model. Use Bbar when comparing influence rankings with the EVBS fit.

Examples

data(itajai)
X <- cbind(1, itajai$pressure - mean(itajai$pressure))
fit <- gevreg.fit(X, itajai$wind)
d <- cnc_diagnostics_gev(fit)
head(order(d$Bj, decreasing = TRUE))

---

Density of the Extreme-Value Birnbaum-Saunders distribution

Description

Probability density function of the EVBS distribution on the response (positive) scale. The EVBS variate $T$ satisfies $Z = (2/\alpha)\,\sinh\{(\log T - \eta)/2\} \sim \mathrm{GEV}(0,1,\gamma)$.

Usage

devbs(t, eta, alpha, gama, log = FALSE)

Arguments

t	Vector of positive quantiles.
eta	Location on the log scale, typically $x^\top \beta$.
alpha	Positive shape (scale) parameter.
gama	Extreme-value shape parameter.
log	Logical; if TRUE, the log-density is returned.

Details

The response $T$ lies in the same max-domain of attraction as the GEV variate. Its extreme-value index is $2\gamma$ when $\gamma > 0$, because $T \approx \beta\alpha^2 Z^2$ in the upper tail, and $\gamma$ itself when $\gamma \le 0$: near a finite endpoint the transformation is locally linear (Ferreira, Gomes and Leiva, 2012, Theorems 4.1 and 4.2).

Value

A numeric vector of (log-)density values.

Examples

devbs(c(10, 20, 30), eta = 2.58, alpha = 0.19, gama = -0.16)

---

Density of the log-EVBS distribution (GAMLSS parametrization)

Description

Density of the log-EVBS distribution (GAMLSS parametrization)

Usage

dlogEVBS(x, mu = 0, sigma = 1, nu = 0, log = FALSE)

plogEVBS(q, mu = 0, sigma = 1, nu = 0, lower.tail = TRUE, log.p = FALSE)

qlogEVBS(p, mu = 0, sigma = 1, nu = 0, lower.tail = TRUE, log.p = FALSE)

Arguments

x, q, p	Vector of quantiles / probabilities.
mu	Location on the log scale.
sigma	Positive scale parameter.
nu	Extreme-value shape parameter.
log, log.p	Logical; return the log-density / log-probability.
lower.tail	Logical; if TRUE, probabilities are $P(Y \le y)$.

Value

A numeric vector.

---

Normal Probability Plot with Simulation Envelope

Description

Produces a normal probability (QQ) plot of the randomized quantile residuals together with a simulated envelope, reproducing the residual diagnostic figure of the paper. Points falling outside the envelope indicate poor fit.

Usage

envelope_qq(X, t, nrep = 100)

Arguments

X	A numeric design matrix with an intercept column.
t	A numeric vector of strictly positive responses.
nrep	Number of simulated samples used to build the envelope (default 100).

Value

Invisibly, a list with components r (observed residuals), e1 and e2 (lower and upper envelope bounds), and med (envelope median). A base-graphics plot is produced as a side effect.

References

Atkinson, A. C. (1985). Plots, Transformations and Regression. Oxford University Press.

See Also

rqrandomized.

Examples

data(itajai)
X <- cbind(1, itajai$pressure)
envelope_qq(X, itajai$wind, nrep = 100)

---

Block bootstrap standard errors for the EVBS regression model

Description

Nonparametric moving-block bootstrap for the parameters of the EVBS regression model. Intended for series with residual serial dependence, where the standard errors obtained from the observed information under the independence assumption are optimistic.

Usage

evbs_block_boot(X, y, L = 4, B = 500)

Arguments

X	Design matrix (including the intercept column).
y	Response vector.
L	Block length.
B	Number of bootstrap resamples.

Value

A list with the bootstrap standard errors and the matrix of bootstrap estimates.

Examples

data(itajai)
evbs_block_boot(cbind(1, itajai$pressure), itajai$wind, L = 4, B = 100)$se

---

Finite upper endpoint of the EVBS distribution

Description

When $\gamma < 0$ the EVBS lies in the Weibull max-domain of attraction and its support is bounded above. The endpoint is the largest response the fitted model regards as possible, and is therefore the quantity of primary interest in design applications. It is also highly sensitive to individual observations, which is what the local influence diagnostics in this package are designed to detect.

Usage

evbs_endpoint(eta, alpha, gama)

Arguments

eta	Location on the log scale, typically $x^\top \beta$.
alpha	Positive shape (scale) parameter.
gama	Extreme-value shape parameter.

Value

The finite upper endpoint, or Inf when $\gamma \ge 0$.

Examples

evbs_endpoint(eta = 2.58, alpha = 0.19, gama = -0.16)
evbs_endpoint(eta = 2.58, alpha = 0.19, gama = 0.10) # Inf

---

Prospective endpoint-identifiability index

Description

Computes the monitoring index $F_t$ that combines the conformal normal curvature diagnostic with the finite upper endpoint of the fitted model. The index compares the support the model would infer without the observation flagged by the curvature against the largest response actually observed.

Usage

evbs_monitor(X, y, burn = 40, use = c("curvature", "max"))

Arguments

X	Design matrix (including the intercept column).
y	Response vector, in time order.
burn	Number of initial observations before the index is first computed.
use	Either "curvature" (default) or "max".

Details

Let $j^{*}(t)$ be the observation with the largest aggregate contribution when the model is fitted to the first $t$ records, and $y_{\max}(t)$ the largest response up to $t$. Then

$F_t = \hat T^{*}_{(-j^{*}(t))}(t) / y_{\max}(t),$

where $\hat T^{*}_{(-j^{*})}$ is the endpoint estimated after deleting the flagged observation, evaluated at the covariate value of the observation that attains $y_{\max}(t)$, so that numerator and denominator refer to the same covariate value. The flagged observation is the one with the largest mean aggregate contribution, colMeans(cnc_diagnostics(fit)$Bj). See gev_monitor for the same index under the GEV regression.

Values below one indicate that the model fitted without the flagged observation assigns probability zero to a value that was in fact observed: the endpoint is then determined by a single record and should not be quoted as a design value. The control limit at one follows from the definition and is not a tuning constant.

Deleting the sample maximum instead of the flagged observation (use = "max") lowers the endpoint essentially by construction and produces an index that signals almost everywhere; it is provided for comparison only.

Value

A data frame with the time index, the running maximum, the flagged observation, the fitted shape parameter and the index F.

Examples

data(itajai)
ct <- evbs_monitor(cbind(1, itajai$pressure), itajai$wind, burn = 40)
plot(ct$t, ct$F, type = "l"); abline(h = 1, lty = 2)

---

Return level and expected shortfall for the EVBS regression model

Description

Computes the period-observation return level and, optionally, the conditional tail expectation (expected shortfall) at a given linear predictor value.

Usage

evbs_return_level(object, x, period, es = FALSE)

Arguments

object	An object of class evbsreg from evbsreg.fit.
x	Covariate vector (including the intercept) at which to evaluate.
period	Return period, expressed in the units of one observation. For monthly maxima and a $T$-year return level, use period = 12 * T.
es	Logical; if TRUE, also return the expected shortfall.

Details

The return level is the quantile of order $1 - 1/\code{period}$ of the fitted conditional distribution of one observation at covariate value x: it is exceeded on average once in period observations with that covariate value, not a return level of the marginal distribution. With monthly maxima, period = 12 * T gives the level exceeded once in $12T$ months; reading it as the $T$-year level of the annual maximum assumes the monthly maxima of a year are independent. Periods much longer than the record are extrapolations.

Value

A named numeric vector with the return level and, if requested, the expected shortfall.

Examples

data(itajai)
fit <- evbsreg.fit(cbind(1, itajai$pressure), itajai$wind)
# 50-year return level at mean pressure, from monthly maxima
evbs_return_level(fit, x = c(1, mean(itajai$pressure)),
                  period = 12 * 50, es = TRUE)

---

Fit the Extreme-Value Birnbaum-Saunders Regression Model

Description

Fits the log-Extreme-Value Birnbaum-Saunders (log-EVBS) regression model by joint maximum likelihood estimation. All parameters $(\beta^\top, \alpha, \gamma)^\top$ are estimated simultaneously using the BFGS algorithm with an analytic score function. Standard errors are computed from the analytic observed Fisher information matrix, which is numerically stable even for ill-conditioned design matrices.

Usage

evbsreg.fit(x, t)

Arguments

x	A numeric design matrix of dimension $n \times p$. It must include an intercept column (a column of ones). Each row corresponds to one observation and each column to one covariate.
t	A numeric vector of length $n$ containing the strictly positive responses (e.g. wind gust speeds). Internally the model is fit to log(t).

Details

The EVBS regression model links the location parameter of the log-EVBS distribution to a linear predictor $\mu_i = x_i^\top \beta$. The shape parameter $\alpha > 0$ controls dispersion and the tail-shape parameter $\gamma$ governs the Generalized Extreme Value tail behaviour (Frechet for $\gamma>0$, Gumbel for $\gamma=0$, Weibull for $\gamma<0$).

Initial values are obtained from lm.fit applied to the log response, together with the moment-based starting point $\alpha_0 = \sqrt{(4/n)\sum \sinh^2((y_i - x_i^\top \beta_0)/2)}$ and $\gamma_0 = 0.01$. Optimization uses optim with method = "BFGS", the analytic score, and tolerance reltol = 1e-12.

The observed Fisher information is assembled analytically (see the package vignette and the paper's appendix) and inverted via a Cholesky factorization, falling back to solve if the matrix is not positive definite.

Value

A list of class "evbsreg" with components:

coeff

Numeric vector of length $p+2$: the full parameter vector $(\beta_0,\ldots,\beta_{p-1}, \alpha, \gamma)$.

betahat

Numeric vector of length $p$: the regression coefficients.

alphahat

Scalar: the estimated shape parameter $\hat\alpha$.

gamahat

Scalar: the estimated tail-shape parameter $\hat\gamma$.

stderrors

Numeric vector of length $p$: standard errors of the regression coefficients.

stderroralpha

Scalar: standard error of $\hat\alpha$.

stderrorgama

Scalar: standard error of $\hat\gamma$.

zstats

Numeric vector: Wald z-statistics for the regression coefficients.

pvalues

Numeric vector: two-sided p-values for the regression coefficients.

muhat

Numeric vector of length $n$: fitted linear predictor on the log scale.

xi1, xi2

Numeric vectors of length $n$: helper quantities $\xi_{i1}, \xi_{i2}$ evaluated at the MLE.

observmatrix

Numeric matrix of dimension $(p+2)\times(p+2)$: the analytic Hessian $\ddot\ell(\hat\theta)$ of the log-likelihood (negative definite at the maximum).

hessian

Numeric matrix: identical to observmatrix; provided under the conventional name. The observed Fisher information is -hessian.

inv

Numeric matrix of dimension $(p+2)\times(p+2)$: the inverse of the observed Fisher information -hessian, i.e. the asymptotic variance-covariance matrix of $\hat\theta$.

B

Numeric matrix of dimension $n\times n$: the influence matrix $B = \Delta^\top (-\ddot\ell)^{-1} \Delta$ for the case-weight perturbation scheme, consumed by cnc_diagnostics.

nobs

Integer: the number of observations $n$.

npar

Integer: the number of parameters $p+2$.

The returned object has class "evbsreg" and has a print.evbsreg method that displays a coefficient table.

References

Ospina, R., Lima, J. I. C., Barros, M., and Macedo, A. M. S. (2026). Local influence diagnostics for the extreme-value Birnbaum-Saunders regression model: methodology, validation, and application to anomalous wind gusts. Submitted.

Leiva, V., Ferreira, M., Gomes, M. I., and Lillo, C. (2016). Extreme value Birnbaum-Saunders regression models applied to environmental data. Stochastic Environmental Research and Risk Assessment, 30, 1045–1058.

Cook, R. D. (1986). Assessment of local influence. Journal of the Royal Statistical Society, Series B, 48, 133–169.

See Also

cnc_diagnostics for influence diagnostics, rqrandomized for residuals, itajai for the example dataset.

Examples

data(itajai)
X <- cbind(1, itajai$pressure)
fit <- evbsreg.fit(X, itajai$wind)

## Parameter estimates and standard errors
round(fit$coeff, 4)
round(c(fit$stderrors, fit$stderroralpha, fit$stderrorgama), 4)

---

Monte Carlo Simulation Study for the EVBS Regression Model

Description

Runs a Monte Carlo simulation that evaluates the finite-sample properties of the joint maximum likelihood estimator of the EVBS regression model under one of three scenarios.

Usage

evbsreg.fit.mc(
  m,
  n,
  beta0,
  beta1,
  alpha,
  gama,
  scenario = c("canonical", "leverage", "robustness"),
  semente = 2023
)

Arguments

m	Number of Monte Carlo replicates (the paper uses 5000; set to a smaller value such as 500 for a quick check).
n	Sample size for each replicate (the paper uses 60, 120, 180).
beta0	True intercept $\beta_0$.
beta1	True slope $\beta_1$.
alpha	True shape parameter $\alpha$.
gama	True tail-shape parameter $\gamma$.
scenario	Character string selecting the design: "canonical" covariate $x \sim U(0,1)$, no contamination (baseline). "leverage" 10% of covariate values drawn from $U(5,10)$ to introduce high-leverage points. "robustness" 10% of observations generated with the shape parameter shifted by $-0.5$ (alpha contamination).
semente	Integer RNG seed for reproducibility (default 2023).

Value

A numeric matrix of dimension $10 \times 4$. Columns are the four parameters Beta0, Beta1, Alpha, Gama; rows are: Parametro (true value), EMV (mean estimate), VIES-ABS (absolute bias), VIES-REL (relative bias), VAR (empirical variance), EQM (mean squared error), E-PADRAO (empirical standard error), EP-FISHER (mean Fisher standard error), RAIZ-EQM (root mean squared error), and TAXA-COB (empirical 95% coverage rate).

References

Ospina, R., Lima, J. I. C., Barros, M., and Macedo, A. M. S. (2026). Local influence diagnostics for the extreme-value Birnbaum-Saunders regression model. Submitted.

See Also

revbs, evbsreg.fit.

Examples

## Quick check with m = 50 replicates
res <- evbsreg.fit.mc(m = 50, n = 60,
                       beta0 = 0.5, beta1 = 0.5,
                       alpha = 0.5, gama = 0.20,
                       scenario = "canonical")
print(res)

---

Generate EVBS Density Values

Description

Computes $(t, f(t))$ pairs of the Extreme-Value Birnbaum-Saunders (EVBS) density over its support, for given parameters.

Usage

generate_evbs_data(alpha, beta, gama)

Arguments

alpha	Shape parameter $\alpha > 0$.
beta	Scale parameter $\beta > 0$.
gama	Tail-shape parameter $\gamma$ (real). The support is bounded below when $\gamma > 0$ and above when $\gamma < 0$.

Value

A data.frame with columns t (support points) and y (density values).

References

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

See Also

plot_evbs_alpha, generate_logevbs_data.

Examples

d <- generate_evbs_data(alpha = 0.5, beta = 1, gama = 0.5)
plot(d$t, d$y, type = "l", xlab = "t", ylab = "f(t)")

---

Generate log-EVBS Density Values

Description

Computes $(y, f(y))$ pairs of the log-EVBS density over its support, for given parameters.

Usage

generate_logevbs_data(alpha, eta, gama)

Arguments

alpha	Shape parameter $\alpha > 0$.
eta	Location parameter $\eta$ (real).
gama	Tail-shape parameter $\gamma$ (real).

Value

A data.frame with columns y_vals (support points) and fdp (density values).

References

Leiva, V., Ferreira, M., Gomes, M. I., and Lillo, C. (2016). Extreme value Birnbaum-Saunders regression models applied to environmental data. Stochastic Environmental Research and Risk Assessment, 30, 1045–1058.

See Also

plot_logevbs_alpha, generate_evbs_data.

Examples

d <- generate_logevbs_data(alpha = 1, eta = 0, gama = 0.5)
plot(d$y_vals, d$fdp, type = "l", xlab = "y", ylab = "f(y)")

---

Prospective endpoint-identifiability index under the GEV regression

Description

The index $F_t$ of evbs_monitor, computed with the GEV regression of gevreg.fit in place of the EVBS model. At each $t \ge$ burn the model is fitted to the first $t$ records, the observation with the largest mean aggregate contribution (cnc_diagnostics_gev(fit)$Bbar) is deleted, the model is refitted, and the endpoint $\mu(x) - \sigma/\xi$ of the refit, evaluated at the covariate value of the running maximum, is divided by the running maximum. The endpoint is infinite, and so is the index, when the refitted $\xi$ is not negative.

Usage

gev_monitor(X, y, burn = 40, use = c("curvature", "max"))

Arguments

X	Design matrix (including the intercept column). Centre the covariates: the GEV likelihood is poorly conditioned otherwise.
y	Response vector, in time order.
burn	Number of initial observations before the index is first computed.
use	Either "curvature" (default) or "max".

Value

A data frame of class evbs_monitor with the time index, the running maximum, the flagged observation, the fitted shape parameter xi, the endpoint and the index F. Instants at which a fit fails or does not converge are omitted.

Examples

data(itajai)
X <- cbind(1, itajai$pressure - mean(itajai$pressure))
gt <- gev_monitor(X, itajai$wind, burn = 40)
plot(gt)

---

Individual score contributions for the GEV regression model

Description

Returns the matrix whose $i$th row is the score vector of the $i$th observation, $s_i(\theta) = \partial \ell_i / \partial \theta$, for the generalized extreme-value regression model with linear location, $Y_i \sim \mathrm{GEV}(\mu_i, \sigma, \xi)$ and $\mu_i = x_i^\top \beta$.

Usage

gev_scores(y, X, beta, sigma, xi)

Arguments

y	Response vector.
X	Design matrix (including the intercept column).
beta	Regression coefficients for the location parameter.
sigma	Positive scale parameter.
xi	Extreme-value shape parameter.

Details

This matrix is the key to generalizing local influence diagnostics beyond a single model class. Under the case-weight perturbation scheme the perturbed log-likelihood is $\ell(\theta \mid \omega) = \sum_i \omega_i \ell_i(\theta)$, so the perturbation matrix is

$\Delta = \frac{\partial^2 \ell(\theta \mid \omega)}{\partial \theta \, \partial \omega^\top} = [\, s_1(\theta) \; \cdots \; s_n(\theta) \,],$

evaluated at the maximum likelihood estimate. That is, $\Delta$ is simply the matrix of individual score contributions. This holds for any likelihood model, which is why the conformal normal curvature framework requires nothing beyond the score and the observed information.

With $z_i = (y_i - \mu_i)/\sigma$, $u_i = 1 + \xi z_i > 0$ and $w_i = u_i^{-1/\xi}$, the components are

$\partial \ell_i / \partial \mu_i = (\xi + 1 - w_i)/(\sigma u_i),$

$\partial \ell_i / \partial \sigma = -1/\sigma + z_i(\xi + 1 - w_i)/(\sigma u_i),$

$\partial \ell_i / \partial \xi = \xi^{-2}(1 - w_i)\log u_i - (z_i/u_i)\{1 + (1 - w_i)/\xi\}.$

Value

An $n \times (p+2)$ matrix of individual score contributions.

Examples

data(itajai)
X <- cbind(1, itajai$pressure - mean(itajai$pressure))
fit <- gevreg.fit(X, itajai$wind)
S <- gev_scores(itajai$wind, X, fit$beta, fit$sigma, fit$xi)
colSums(S) # approximately zero at the MLE

---

Maximum likelihood fit of the GEV regression model

Description

Fits $Y_i \sim \mathrm{GEV}(\mu_i, \sigma, \xi)$ with $\mu_i = x_i^\top \beta$ by maximum likelihood, using the analytic gradient and several starting values for the shape parameter.

Usage

gevreg.fit(X, y, xi_start = c(-0.2, -0.05, 0.05, 0.2))

Arguments

X	Design matrix (including the intercept column).
y	Response vector.
xi_start	Vector of starting values for the shape parameter.

Details

Centring the covariates is strongly recommended: the GEV likelihood is poorly conditioned when covariates are far from the origin, and an uncentred fit may fail to converge without any warning.

Value

A list with the estimates, the log-likelihood, the observed information (Hessian of the negative log-likelihood) and the convergence code.

Examples

data(itajai)
X <- cbind(1, itajai$pressure - mean(itajai$pressure))
gevreg.fit(X, itajai$wind)

---

Monthly Maximum Wind Gusts at Itajai, Brazil

Description

Monthly maximum wind gust speed and the daily mean atmospheric pressure on the day of each maximum, recorded at INMET station A-868 (latitude $-26.95083^\circ$, longitude $-48.76194^\circ$) in Itajai, Santa Catarina, Brazil, from July 2010 to October 2020.

Usage

itajai

Format

A data frame with 124 rows and 3 variables:

month

Integer index of the monthly observation (1 = July 2010, ..., 124 = October 2020).

wind

Numeric. Monthly maximum wind gust speed (m/s).

pressure

Numeric. Daily mean atmospheric pressure (mb) on the day of the monthly maximum.

Details

Observation 82 corresponds to the catastrophic event of April 26, 2017 (wind = 33.9 m/s), a severe cold front that struck the coast of Santa Catarina, with fatalities reported by local civil-defense authorities. This observation is identified as highly influential on the tail-shape parameter by the diagnostics implemented in this package.

Source

Instituto Nacional de Meteorologia (INMET), Brazil. Data access: https://www.inmet.gov.br/

References

Ospina, R., Lima, J. I. C., Barros, M., and Macedo, A. M. S. (2026). Local influence diagnostics for the extreme-value Birnbaum-Saunders regression model. Submitted.

Examples

data(itajai)
str(itajai)
plot(itajai$pressure, itajai$wind, pch = 16, col = "darkgreen",
     xlab = "Pressure (mb)", ylab = "Wind gust (m/s)")

---

The log-EVBS distribution for GAMLSS models

Description

Defines the log-Extreme-Value Birnbaum-Saunders distribution as a three-parameter gamlss family, so that each parameter may depend on covariates through a linear predictor or a smooth term. This is the natural way to let the tail-shape parameter vary with covariates, which the fixed three-parameter fit of evbsreg.fit does not allow.

Usage

logEVBS(mu.link = "identity", sigma.link = "log", nu.link = "identity")

Arguments

mu.link	Link for the location parameter $\mu$ (default "identity"; the predictor is on the log-response scale).
sigma.link	Link for the scale parameter $\alpha > 0$ (default "log").
nu.link	Link for the extreme-value shape parameter $\gamma$ (default "identity").

Details

The response is $Y = \log T$ where $T$ follows the EVBS distribution, so that $Z = (2/\sigma)\sinh\{(Y-\mu)/2\} \sim \mathrm{GEV}(0,1,\nu)$. The parameter map is $\mu = \eta = x^\top\beta$ (identity link), $\sigma = \alpha$ (log link, positive) and $\nu = \gamma$ (identity link, real-valued).

With a constant predictor the fitted $(\mu, \sigma, \nu)$ reproduce the $(\eta, \alpha, \gamma)$ of evbsreg.fit; this equivalence is the recommended check after loading the family. The score and second-order derivatives required by gamlss are obtained by automatic numerical differentiation of the analytic log-density, which is exact to working precision and avoids transcription error in the lengthy closed forms.

Value

A gamlss.family object.

Examples

## Not run: 
  library(gamlss)
  data(itajai)
  y <- log(itajai$wind)
  # constant model: should match evbsreg.fit on the log scale
  m0 <- gamlss(y ~ 1, family = logEVBS())
  # tail shape depending on pressure:
  m1 <- gamlss(y ~ pressure, nu.formula = ~ pressure,
               data = itajai, family = logEVBS())

## End(Not run)

---

Distribution function of the Extreme-Value Birnbaum-Saunders distribution

Description

Distribution function of the Extreme-Value Birnbaum-Saunders distribution

Usage

pevbs(t, eta, alpha, gama, lower.tail = TRUE)

Arguments

t	Vector of positive quantiles.
eta	Location on the log scale, typically $x^\top \beta$.
alpha	Positive shape (scale) parameter.
gama	Extreme-value shape parameter.
lower.tail	Logical; if TRUE (default), probabilities are $P(T \le t)$, otherwise $P(T > t)$.

Value

A numeric vector of probabilities.

Examples

pevbs(30, eta = 2.58, alpha = 0.19, gama = -0.16, lower.tail = FALSE)

---

Plot Aggregate Contributions B_j(q)

Description

Produces panel (b) of the diagnostic figure: the aggregate contributions $B_j(q)$ of each observation, with a horizontal reference line at $b(q)$ and automatic labelling of the most influential points.

Usage

plot_aggregate_contributions(
  diag,
  q = 7,
  label.flagged = 5,
  pch = 16,
  cex = 0.7,
  main = ""
)

Arguments

diag	A list returned by cnc_diagnostics.
q	Integer influence threshold in $1, \ldots, 7$. The paper uses q = 7.
label.flagged	Number of largest observations to label.
pch	Plotting character (default 16).
cex	Point size expansion (default 0.7).
main	Plot title (default empty).

Value

The indices of the labelled observations, returned invisibly. Also produces a base-graphics plot as a side effect.

See Also

plot_cnc, cnc_diagnostics.

Examples

data(itajai)
fit <- evbsreg.fit(cbind(1, itajai$pressure), itajai$wind)
diag <- cnc_diagnostics(fit)
plot_aggregate_contributions(diag, q = 7, main = "(b)")

---

Combined Conformal Normal Curvature Diagnostic Plot

Description

Produces the two-panel diagnostic figure of the paper: normalized eigenvalues (left) and aggregate contributions $B_j(q)$ (right), side by side.

Usage

plot_cnc(diag, q = 7, label.flagged = 5)

Arguments

diag	A list returned by cnc_diagnostics.
q	Integer influence threshold in $1, \ldots, 7$ (default 7).
label.flagged	Number of largest observations to label in the right panel (default 5).

Value

Called for its side effect (a two-panel base-graphics figure). Returns NULL invisibly.

See Also

cnc_diagnostics, plot_normalized_eigenvalues, plot_aggregate_contributions.

Examples

data(itajai)
fit <- evbsreg.fit(cbind(1, itajai$pressure), itajai$wind)
diag <- cnc_diagnostics(fit)
plot_cnc(diag, q = 7)

---

Plot EVBS Densities for Varying alpha

Description

Reproduces Figure 1(a) of the paper: EVBS densities for $\alpha \in \{0.25, 0.5, 0.75, 1\}$, $\beta = 1$, $\gamma = 0.5$. Uses ggplot2 with the Dark2 colour palette.

Usage

plot_evbs_alpha()

Value

A ggplot object.

See Also

plot_evbs_gama, generate_evbs_data.

Examples

print(plot_evbs_alpha())

---

Plot EVBS Densities for Varying gamma

Description

Reproduces Figure 1(b) of the paper: EVBS densities for $\gamma \in \{-1.25, -0.75, 0.75, 1.25\}$, $\alpha = 0.5$, $\beta = 1.5$. Uses ggplot2 with the Dark2 colour palette.

Usage

plot_evbs_gama()

Value

A ggplot object.

See Also

plot_evbs_alpha.

Examples

print(plot_evbs_gama())

---

Plot log-EVBS Densities for Varying alpha

Description

Reproduces Figure 2(a) of the paper: log-EVBS densities for $\alpha \in \{0.25, 0.5, 2, 4\}$, $\eta = 0$, $\gamma = 0$. Uses ggplot2 with the Dark2 colour palette.

Usage

plot_logevbs_alpha()

Value

A ggplot object.

See Also

plot_logevbs_gama.

Examples

print(plot_logevbs_alpha())

---

Plot log-EVBS Densities for Varying gamma

Description

Reproduces Figure 2(b) of the paper: log-EVBS densities for $\gamma \in \{-1.05, -0.5, 0.5, 1.05\}$, $\alpha = 1$, $\eta = 0$. Uses ggplot2 with the Dark2 colour palette.

Usage

plot_logevbs_gama()

Value

A ggplot object.

See Also

plot_logevbs_alpha.

Examples

print(plot_logevbs_gama())

---

Plot Normalized Eigenvalues of the Influence Matrix

Description

Produces panel (a) of the diagnostic figure: the normalized eigenvalues of the influence matrix, with horizontal reference thresholds for $q = 1, \ldots, 7$.

Usage

plot_normalized_eigenvalues(diag, pch = 16, cex = 1, main = "")

Arguments

diag	A list returned by cnc_diagnostics.
pch	Plotting character (default 16).
cex	Point size expansion (default 1).
main	Plot title (default empty).

Value

Called for its side effect (a base-graphics plot). Returns NULL invisibly.

See Also

plot_cnc, cnc_diagnostics.

Examples

data(itajai)
fit <- evbsreg.fit(cbind(1, itajai$pressure), itajai$wind)
diag <- cnc_diagnostics(fit)
plot_normalized_eigenvalues(diag, main = "(a)")

---

Plot an endpoint-identifiability control chart

Description

Plot an endpoint-identifiability control chart

Usage

## S3 method for class 'evbs_monitor'
plot(x, ...)

Arguments

x	An object returned by evbs_monitor.
...	Further arguments passed to plot.

Value

Called for its side effect; invisibly returns x.

---

Print Method for EVBS Regression Fits

Description

Compactly prints the parameter estimates, standard errors, and Wald tests of an object returned by evbsreg.fit.

Usage

## S3 method for class 'evbsreg'
print(x, digits = 4, ...)

Arguments

x	An object of class "evbsreg".
digits	Number of significant digits (default 4).
...	Further arguments passed to print.default.

Value

The object x, invisibly. Called for the side effect of printing a coefficient table to the console.

See Also

evbsreg.fit.

Examples

data(itajai)
fit <- evbsreg.fit(cbind(1, itajai$pressure), itajai$wind)
print(fit)

---

Quantile function of the Extreme-Value Birnbaum-Saunders distribution

Description

Exact quantile function on the response scale. This is the function needed to compute return levels; the generalized extreme-value return-level formula does not apply to the EVBS response, since the response is a monotone transformation of a GEV variate rather than a GEV variate itself.

Usage

qevbs(p, eta, alpha, gama)

Arguments

p	Vector of probabilities in $(0,1)$.
eta	Location on the log scale, typically $x^\top \beta$.
alpha	Positive shape (scale) parameter.
gama	Extreme-value shape parameter.

Value

A numeric vector of quantiles.

Examples

# 50-year return level from monthly maxima: p = 1 - 1/(12*50)
qevbs(1 - 1 / 600, eta = 2.58, alpha = 0.19, gama = -0.16)

---

Cox-Snell Residuals for the Log-EVBS Regression Model

Description

Computes Cox-Snell residuals for a log-EVBS regression fit. Under a correctly specified model these residuals form an approximately standard exponential sample.

Usage

rcoxsnell(X, t)

Arguments

X	A numeric design matrix with an intercept column.
t	A numeric vector of strictly positive responses.

Value

A numeric vector of length length(t) containing the Cox-Snell residuals.

References

Cox, D. R. and Snell, E. J. (1968). A general definition of residuals. Journal of the Royal Statistical Society, Series B, 30, 248–275.

See Also

rqrandomized, envelope_qq.

Examples

data(itajai)
X <- cbind(1, itajai$pressure)
cs <- rcoxsnell(X, itajai$wind)
summary(cs)

---

Random Number Generation from the EVBS Distribution

Description

Generates random variates from the Extreme-Value Birnbaum-Saunders (EVBS) distribution by transforming Generalized Extreme Value (GEV) variates.

Usage

revbs(n, alpha, beta, gama)

Arguments

n	Sample size (number of variates to generate).
alpha	Shape parameter $\alpha > 0$.
beta	Scale parameter $\beta > 0$.
gama	Tail-shape parameter $\gamma$ (real).

Details

If $Z \sim \mathrm{GEV}(0, 1, \gamma)$, then $T = \beta\{1 + \alpha^2 Z^2/2 + \alpha Z \sqrt{1 + \alpha^2 Z^2/4}\}$ follows the EVBS distribution. GEV variates are drawn via rgev from the SpatialExtremes package.

Value

A numeric vector of n EVBS variates. Returns NA if n is NA.

References

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

See Also

evbsreg.fit.mc, evbsreg.fit.

Examples

set.seed(2023)
x <- revbs(100, alpha = 0.5, beta = 1, gama = 0.2)
summary(x)

---

Randomized Quantile Residuals for the Log-EVBS Regression Model

Description

Computes randomized quantile residuals (Dunn and Smyth, 1996) for a log-EVBS regression fit. Under a correctly specified model these residuals are approximately standard normal, so departures from normality indicate lack of fit.

Usage

rqrandomized(X, t)

Arguments

X	A numeric design matrix with an intercept column.
t	A numeric vector of strictly positive responses.

Value

A numeric vector of length length(t) containing the randomized quantile residuals.

References

Dunn, P. K. and Smyth, G. K. (1996). Randomized quantile residuals. Journal of Computational and Graphical Statistics, 5, 236–244.

See Also

rcoxsnell, envelope_qq, evbsreg.fit.

Examples

data(itajai)
X <- cbind(1, itajai$pressure)
r <- rqrandomized(X, itajai$wind)
shapiro.test(r)

---