vasicekreg

CRAN status

vasicekreg provides distribution functions and GAMLSS regression families for Vasicek-type distributions on the unit interval. The package supports mean and quantile regression under the standard normal kernel and quantile regression under the logistic and hyperbolic-secant kernels. Normal-kernel families augmented at zero, at one, or at both boundaries are available for responses containing exact boundary values. Dedicated two-part random-intercept functions are also available for zero-augmented longitudinal responses.

Available families

Family Kernel Response range Parameterization Interpretation of mu
NVASIM Standard normal (0, 1) Mean mu = E(Y)
NVASIQ Standard normal (0, 1) Quantile mu = Q_Y(quantile)
LVASIQ Logistic (0, 1) Quantile mu = Q_Y(quantile)
HVASIQ Hyperbolic secant (0, 1) Quantile mu = Q_Y(quantile)
ZANVASIM Standard normal [0, 1) Zero-augmented mean mu = E(Y | Y > 0)
OANVASIM Standard normal (0, 1] One-augmented mean mu = E(Y | Y < 1)
ZOANVASIM Standard normal [0, 1] Zero-and-one-augmented mu = E(Y | 0 < Y < 1)

The documentation uses augmented for boundary mixtures and Vasicek-type for the kernel-based constructions. The established exported names are kept for backward compatibility and retain the familiar GAMLSS abbreviations ZA and OA, historically read as zero- and one-adjusted.

For all seven families, sigma is a shape parameter in (0, 1) that controls dispersion. For ZANVASIM, nu = P(Y = 0) and the marginal mean is E(Y) = (1 - nu) * mu. For OANVASIM, nu = P(Y = 1) and the marginal mean is E(Y) = nu + (1 - nu) * mu. In ZOANVASIM, nu = P(Y = 0) and tau = P(Y = 1 | Y > 0). Its marginal mean is E(Y) = (1 - nu) * (tau + (1 - tau) * mu). The logistic- and hyperbolic-secant-kernel means do not have closed-form expressions and generally differ from mu; consequently, the package does not provide mean-regression families for these kernels.

Each parameterization includes density, cumulative distribution, quantile, and random generation functions:

Longitudinal two-part mixed models

The functions zabr() and zavr() fit two-part models to repeated responses in [0, 1):

Function Presence component Positive component Random effects
zabr() Logistic Beta mean and precision Independent random intercept in each component
zavr() Logistic NVASIM mean and shape Independent random intercept in each component

These are dedicated maximum-likelihood functions rather than gamlss.family objects. Subject likelihoods are evaluated by non-adaptive Gauss–Hermite quadrature. The preferred interface uses one-sided formulas for both components and a random-intercept formula:

fit_longitudinal <- zavr(
  data = longitudinal_data,
  y = "abundance",
  formula_bin = ~ treatment + time,
  formula_cont = ~ treatment + time,
  random = ~ 1 | subject,
  time_ind = "time",
  quad_n = 30
)

coef(fit_longitudinal)
vcov(fit_longitudinal)
logLik(fit_longitudinal)
BIC(fit_longitudinal)

time_ind identifies visits and is used to detect duplicated subject–time pairs. A time trend is fitted only when the time variable is included in formula_bin and/or formula_cont. Joint likelihood-ratio tests require identical design matrices in the two components. The joint AIC and BIC use the complete two-part likelihood. BIC uses the number of subjects in its penalty, following the PROC NLMIXED convention. Component-wise AIC values sum to the joint AIC. Component-wise BIC values do not generally sum to the joint BIC because the positive component uses only subjects with at least one positive response.

For both longitudinal models, gamma denotes the discrete-component fixed effects and beta denotes the continuous-component fixed effects. The names returned by coef() and the corresponding vcov() dimension names use gamma_ and beta_, including the intercepts.

Installation

Install the released version from CRAN:

install.packages("vasicekreg")

Install the development version from GitHub:

install.packages("remotes")
remotes::install_github("jmazucheli/vasicekreg")

Because the package contains compiled C++ code, installation from source requires an appropriate compiler toolchain.

Included datasets

Object Main bounded response Boundary pattern Observations
bodyfat ARMS, LEGS, BODY, ANDROID, and GYNECOID All in (0, 1) 298
aep noinap / los Zero and one 1,383
transport propbiked Zero 60
trees prop One 26

The help pages preserve the original data provenance and the citations associated with each dataset. Use ?bodyfat, ?aep, ?transport, or ?trees for variable definitions, sources, and references.

An application-oriented vignette fits models to every included data set, including all five bodyfat responses and the three boundary patterns:

vignette("vasicekreg", package = "vasicekreg")

Distribution functions

The following example uses the normal-kernel mean parameterization:

library(vasicekreg)

set.seed(123)
y <- rNVASIM(n = 1000, mu = 0.50, sigma = 0.25)

dNVASIM(x = 0.50, mu = 0.50, sigma = 0.25)
pNVASIM(q = 0.50, mu = 0.50, sigma = 0.25)
qNVASIM(p = 0.50, mu = 0.50, sigma = 0.25)

For the quantile parameterizations, the fixed level is supplied through the quantile argument:

qNVASIQ(p = 0.25, mu = 0.60, sigma = 0.25, quantile = 0.25)
qLVASIQ(p = 0.25, mu = 0.60, sigma = 0.25, quantile = 0.25)
qHVASIQ(p = 0.25, mu = 0.60, sigma = 0.25, quantile = 0.25)

All three calls return mu = 0.60 because mu represents the quantile at the fixed level quantile = 0.25.

GAMLSS mean regression

library(gamlss)
library(vasicekreg)

set.seed(123)
dat_mean <- data.frame(
  y = rNVASIM(n = 300, mu = 0.60, sigma = 0.25)
)

fit_mean <- gamlss(
  y ~ 1,
  data = dat_mean,
  family = NVASIM(
    mu.link = "logit",
    sigma.link = "logit"
  ),
  control = gamlss.control(trace = FALSE)
)

fitted(fit_mean, what = "mu")[1]

GAMLSS zero-augmented mean regression

set.seed(123)
dat_zero <- data.frame(
  y = r0NVASIM(n = 500, mu = 0.60, sigma = 0.25, nu = 0.20)
)

fit_zero <- gamlss(
  y ~ 1,
  sigma.formula = ~ 1,
  nu.formula = ~ 1,
  data = dat_zero,
  family = ZANVASIM(
    mu.link = "logit",
    sigma.link = "logit",
    nu.link = "logit"
  ),
  control = gamlss.control(trace = FALSE)
)

fitted(fit_zero, what = "mu")[1]
fitted(fit_zero, what = "sigma")[1]
fitted(fit_zero, what = "nu")[1]

Here mu is the conditional mean of the positive component. The fitted marginal mean is (1 - fitted(fit_zero, what = "nu")) * fitted(fit_zero, what = "mu").

GAMLSS one-augmented and zero-and-one-augmented regression

dat_one <- data.frame(
  y = r1NVASIM(500, mu = 0.60, sigma = 0.25, nu = 0.20)
)

fit_one <- gamlss(
  y ~ 1,
  sigma.formula = ~ 1,
  nu.formula = ~ 1,
  data = dat_one,
  family = OANVASIM(),
  control = gamlss.control(trace = FALSE)
)

dat_boundary <- data.frame(
  y = r01NVASIM(
    500, mu = 0.60, sigma = 0.25, nu = 0.20, tau = 0.25
  )
)

fit_boundary <- gamlss(
  y ~ 1,
  sigma.formula = ~ 1,
  nu.formula = ~ 1,
  tau.formula = ~ 1,
  data = dat_boundary,
  family = ZOANVASIM(),
  control = gamlss.control(trace = FALSE)
)

For ZOANVASIM, p0 = nu, p1 = (1 - nu) * tau, and pc = (1 - nu) * (1 - tau). This parameterization guarantees valid probabilities while retaining logit links for both boundary parameters. Here tau is an estimated model parameter and is unrelated to the fixed quantile argument of NVASIQ(), LVASIQ(), and HVASIQ().

The included aep data provide a real example with both boundary values:

data("aep", package = "vasicekreg")
aep$inappropriate <- with(aep, noinap / los)

fit_aep <- gamlss(
  inappropriate ~ sex + ward + year + age + loglos,
  sigma.formula = ~ loglos,
  nu.formula = ~ loglos,
  tau.formula = ~ loglos,
  family = ZOANVASIM(),
  data = aep,
  control = gamlss.control(n.cyc = 200, trace = FALSE)
)

GAMLSS quantile regression

For NVASIQ(), LVASIQ(), and HVASIQ(), the fixed quantile level is passed directly to the family constructor through quantile. Each fitted object is self-contained, so models at different levels can coexist without changing global variables.

Normal kernel

library(gamlss)
library(vasicekreg)

set.seed(123)
quantile_level <- 0.50

dat_normal <- data.frame(
  y = rNVASIQ(
    n = 300,
    mu = 0.60,
    sigma = 0.25,
    quantile = quantile_level
  )
)

fit_normal <- gamlss(
  y ~ 1,
  data = dat_normal,
  family = NVASIQ(
    quantile = quantile_level,
    mu.link = "logit",
    sigma.link = "logit"
  ),
  control = gamlss.control(trace = FALSE)
)

fitted(fit_normal, what = "mu")[1]

Logistic kernel

set.seed(123)
quantile_level <- 0.25

dat_logistic <- data.frame(
  y = rLVASIQ(
    n = 300,
    mu = 0.60,
    sigma = 0.25,
    quantile = quantile_level
  )
)

fit_logistic <- gamlss(
  y ~ 1,
  data = dat_logistic,
  family = LVASIQ(
    quantile = quantile_level,
    mu.link = "logit",
    sigma.link = "logit"
  ),
  control = gamlss.control(trace = FALSE)
)

fitted(fit_logistic, what = "mu")[1]

Hyperbolic-secant kernel

set.seed(123)
quantile_level <- 0.25

dat_hs <- data.frame(
  y = rHVASIQ(
    n = 300,
    mu = 0.60,
    sigma = 0.25,
    quantile = quantile_level
  )
)

fit_hs <- gamlss(
  y ~ 1,
  sigma.formula = ~ 1,
  data = dat_hs,
  family = HVASIQ(
    quantile = quantile_level,
    mu.link = "logit",
    sigma.link = "logit"
  ),
  control = gamlss.control(trace = FALSE)
)

fitted(fit_hs, what = "mu")[1]

The fixed level is embedded in the family and residual definitions. vasicek_envelope() also inserts that stored value into every bootstrap refit, so it does not depend on a global variable or on a symbol used in the original family call. The generic update() method retains the usual R call semantics; when direct later use of update() is planned, use a literal level in the original call (for example, NVASIQ(quantile = 0.25)) or keep the referenced symbol available and unchanged.

Simulated residual envelopes

The function vasicek_envelope() constructs pointwise simulated envelopes for normalized randomized quantile residuals and generalized Cox–Snell residuals. In every bootstrap replication, a new response is generated under the fitted model, the same GAMLSS model is re-estimated, and the residuals are recalculated and ordered. This accounts for parameter-estimation variability.

For augmented families, the probability integral transform is randomized over the fitted CDF jump at zero and/or one. Quantile residuals are compared with the standard normal distribution, whereas Cox–Snell residuals are compared with the unit-mean exponential distribution. Cox–Snell residuals are evaluated as -pnorm(rq, lower.tail = FALSE, log.p = TRUE), which is equivalent to -log(1 - pnorm(rq)) but remains stable in the upper tail.

envelope_fit <- vasicek_envelope(
  object = fit_mean,
  residual = c("quantile", "cox-snell"),
  nsim = 200,
  level = 0.95,
  envelope = "quantile",
  seed = 123
)

plot(envelope_fit, which = "quantile")
plot(envelope_fit, which = "cox-snell")

The plots are full quantile–quantile plots rather than half-normal plots. The default envelope uses pointwise empirical percentiles and is not a simultaneous confidence band. To use the minimum and maximum of the ordered simulated residuals at each position, set envelope = "minmax". The blue mean curve is calibrated to the fitted model and sample size; for Cox–Snell residuals it may be more informative than the theoretical identity line in the finite-sample upper tail. Failed or nonconverged bootstrap fits are discarded and replaced up to the limit specified by max.attempts. Automatic refitting requires the supplied data to contain exactly the observations used in the original fit.

Implementation

The density, cumulative distribution, and quantile functions call compiled C++ routines through Rcpp. Random generation for NVASIM and NVASIQ uses inverse transformation with the corresponding compiled quantile functions, whereas rLVASIQ() and rHVASIQ() call compiled random-generation routines directly. The boundary-augmented functions reuse the compiled NVASIM functions and add the required point masses in R.

The GAMLSS family definitions and analytical log-likelihood derivatives are implemented in R. Numerical differentiation is not used. The mean and variance components of LVASIQ() and HVASIQ() are evaluated by numerical quadrature because these moments have no closed-form expressions.

Testing

From the package root directory:

Rscript -e 'testthat::test_local(path = ".", reporter = "summary")'

cd ..
R CMD build vasicekreg
R CMD check vasicekreg_1.3.0.tar.gz

Citation

To cite the package, use:

citation("vasicekreg")

References

Vasicek distribution and regression

Longitudinal two-part mixed models

GAMLSS framework

Boundary-augmented distributions

Residuals and simulated envelopes

Included data sources

Computational implementation

License

vasicekreg is distributed under the MIT License.