Package {evbsreg}

Type: Package
Title: Local Influence Diagnostics for the Extreme-Value Birnbaum-Saunders Regression Model
Version: 1.0.0
Date: 2026-06-17
Description: Implements local influence diagnostics for the Extreme-Value Birnbaum-Saunders (EVBS) regression model: joint maximum likelihood estimation, conformal normal curvature diagnostics under three perturbation schemes (case-weight, response variable, and explanatory variable), randomized quantile residuals with simulation envelope, Monte Carlo simulation utilities, and publication-quality density and diagnostic plots. The methods are described in Ospina, Lima, Barros, and Macedo (2026, submitted) and are applied to monthly maximum wind gust data from Itajai, Brazil.
License: MIT + file LICENSE
URL: https://github.com/Raydonal/evbsreg, https://raydonal.github.io/evbsreg/
BugReports: https://github.com/Raydonal/evbsreg/issues
Encoding: UTF-8
LazyData: true
Depends: R (≥ 4.0.0)
Imports: stats, graphics, SpatialExtremes, ggplot2
Suggests: grDevices, knitr, rmarkdown, testthat (≥ 3.0.0)
VignetteBuilder: knitr
Config/roxygen2/version: 8.0.0
NeedsCompilation: no
Packaged: 2026-06-25 21:04:36 UTC; raydonal
Author: Raydonal Ospina ORCID iD [aut, cre]
Maintainer: Raydonal Ospina <raydonal@de.ufpe.br>
Repository: CRAN
Date/Publication: 2026-06-30 20:40:16 UTC

evbsreg: Local Influence Diagnostics for the EVBS Regression Model

Description

Implements local influence diagnostics for the Extreme-Value Birnbaum-Saunders (EVBS) regression model: joint maximum likelihood estimation, conformal normal curvature diagnostics under three perturbation schemes (case-weight, response variable, and explanatory variable), randomized quantile residuals with simulation envelope, Monte Carlo simulation utilities, and publication-quality density and diagnostic plots.

Main functions

evbsreg.fit

Joint MLE of the EVBS regression model.

cnc_diagnostics

Conformal normal curvature diagnostics.

plot_cnc

Two-panel diagnostic figure.

rqrandomized, rcoxsnell, envelope_qq

Residuals and envelope.

revbs, evbsreg.fit.mc

Random generation and Monte Carlo study.

Author(s)

Maintainer: Raydonal Ospina raydonal@de.ufpe.br (ORCID)

Authors:

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.

See Also

Useful links:

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])

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)

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)")

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)")

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)")

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)

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)

mirror server hosted at Truenetwork, Russian Federation.