Package {HZIP}

Envelope Plot for HZIP Models

Description

Produces a normal Q-Q plot with a simulation envelope for randomized quantile residuals from a Hierarchical Zero-Inflated Poisson (HZIP) model fitted via hzip.

Usage

envelope.HZIP(object, nsim = 21, method = "ghq", Q = 15, ...)

Arguments

Details

The envelope is constructed using Monte Carlo simulation based on nsim replications of independent standard normal samples. For each order statistic, the 2.5% and 97.5% empirical quantiles define the envelope limits.

If object is of class HZIP, residuals are computed internally via residuals.HZIP before plotting.

Value

A ggplot2 object containing the Q-Q plot with simulation envelope.

See Also

hzip, residuals.HZIP

Examples

fit.salamander <- hzip(
  y ~ mined | mined + spp,
  data = salamanders
)

# Passing the fitted object directly
envelope.HZIP(fit.salamander)

# Or passing residuals explicitly
res <- residuals(fit.salamander)
envelope.HZIP(res)

Fit a Hierarchical Zero-Inflated Poisson (HZIP) Model

Description

hzip() fits a longitudinal/clustered zero-inflated Poisson model with subject-level random effects following a Generalized Log-Gamma (GLG) distribution, by maximizing a marginal likelihood approximated via adaptive Gauss-Hermite quadrature. The model uses a two-part Formula: y \sim \text{zero part} \mid \text{count part}, where the count intensity (Poisson mean) and the zero-inflation probability are linked to (possibly different) sets of covariates. Initial values are obtained from pscl::zeroinfl(..., dist = "poisson", link = "cloglog").

Usage

hzip(
  formula,
  data,
  hessian = TRUE,
  method = "BFGS",
  Q = 15,
  lower = -Inf,
  upper = Inf,
  control = NULL,
  start = NULL,
  refine = FALSE,
  ep_method = c("optim", "numDeriv"),
  n_starts = 1,
  verbose = FALSE,
  ...
)

Arguments

Details

Let y_{ij} denote the count response for subject i at occasion j. The HZIP model assumes

P(y_{ij}=0 \mid u_i) = \pi_{ij}(u_i) + \{1-\pi_{ij}(u_i)\}\exp\{-\mu_{ij}(u_i)\},

P(y_{ij}=k \mid u_i) = \{1-\pi_{ij}(u_i)\}\frac{\mu_{ij}(u_i)^k e^{-\mu_{ij}(u_i)}}{k!},\quad k\ge 1,

with linear predictors for the count and zero parts (links typically log for the Poisson mean and cloglog for the zero-inflation). Subject-specific random effects u_i follow a GLG distribution and induce within-subject dependence; the marginal likelihood is approximated by adaptive Gauss-Hermite quadrature with Q nodes.

Optimization strategy. The default workflow is:

1. Generate initial values from pscl::zeroinfl (with scales set to a neutral 0.5, ensuring reproducibility);

2. Run optim with the specified method (default BFGS);

3. If refine = TRUE, refine the result with nlminb (off by default; useful for difficult problems);

4. If n_starts > 1, repeat steps 2-3 with perturbed initial values and keep the best log-likelihood;

5. Compute standard errors from the inverse Hessian using the method specified by ep_method.

Value

An object of class "HZIP", a list with elements:

Note

The subject identifier must be named Ind. The internal parameter vector follows the order c(scale_zero, beta_zero, scale_count, beta_count).

For difficult problems (e.g., models with many interaction terms), if the default fit produces large standard errors or the AIC suggests a worse fit than a nested simpler model, try (i) n_starts = 10 for multi-start optimization and (ii) supplying start = ... from a fit of the simpler nested model.

References

Min, Y., & Agresti, A. (2005). Random effect models for repeated measures of zero-inflated count data. Statistical Modelling, 5(1), 1-19.

Jackman, S. (2020). pscl: Classes and Methods for R Developed in the Political Science Computational Laboratory. R package version 1.5.5.

Zeileis, A., & Croissant, Y. (2010). Extended model formulas in R: Journal of Statistical Software, 34(1), 1-13. (Formula)

Examples

# Basic fit
fit.salamander <- hzip(y ~ mined | mined + spp, data = salamanders)
summary(fit.salamander)

# More robust fit for difficult problems (multi-start)
fit.salamander2 <- hzip(y ~ mined | mined + spp + mined:spp,
                        data = salamanders,
                        n_starts = 10, verbose = TRUE)
summary(fit.salamander2)

# Using the simpler model as starting point for the harder one
init <- with(fit.salamander, c(scale_zero, coefficients_zero,
                               scale_count, coefficients_count,
                               rep(0, 6)))  # 6 zeros for the interaction
fit.salamander3 <- hzip(y ~ mined | mined + spp + mined:spp,
                        data = salamanders, start = init)
summary(fit.salamander3)

Print Method for hzip_test Objects

Description

Prints a formatted summary of a simulation-based test returned by testDisp.HZIP or testZI.HZIP, including the null hypothesis, observed statistic, p-value, and decision.

Usage

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

Arguments

Simulate Data from a Hierarchical Zero-Inflated Poisson (HZIP) Model

Description

rHZIP() generates panel/longitudinal data from a two–part hierarchical zero-inflated Poisson model with subject-specific random effects for both the zero-inflation and the count components. Random effects are drawn from a generalized log-gamma (GLG) distribution via rgengamma() (user must provide/attach a function with this name; see Dependencies).

Usage

rHZIP(n, m, para, x1, x2)

Arguments

Details

For subject i=1,\dots,n with m_i observations, the model draws two subject-level random effects b^{(0)}_i and b^{(1)}_i independently from GLG distributions parameterized by lambda1 and lambda2. Conditional on these effects, outcomes are generated as

Y_{ij} = Z_{ij}\times C_{ij},

where Z_{ij}\sim \mathrm{Bernoulli}(p_{ij}) controls structural zeros and C_{ij}\sim \mathrm{Poisson}(\mu_{ij}) controls the count size.

The returned data frame contains subject IDs, the response y, and the supplied covariates. An attribute "propZeros" stores a small summary with the percentage of structural zeros, additional Poisson zeros, and total zeros.

Value

A tibble with columns:

The object has an attribute "propZeros": a 3×1 data.frame with rows "Zeros" (structural zeros, %), "Count" (extra zeros from the Poisson part, %), and "Total" (overall zero percentage).

Column naming

Columns of x1 are renamed to x1, x2, ..., xp1. Columns of x2 are copied except the first column is dropped and the remaining are renamed w1, w2, ..., w_{p2-1}. (This mirrors the current implementation that excludes the first x2 column from the output.)

Dependencies

This function calls rgengamma() to draw GLG random effects. Ensure that such a function is available on the search path (e.g., from a package that provides a generalized log-gamma RNG) or provide your own implementation with the signature rgengamma(n, mu, sigma, lambda). It also uses dplyr and tibble.

See Also

hzip for model fitting.

Examples

set.seed(123)

n <- 50
m <- rep(4, n)
N <- sum(m)

# design matrices (with intercepts)
x1 <- cbind(1, rnorm(N))
x2 <- cbind(1, rnorm(N), rbinom(N, 1, 0.5))

p1 <- ncol(x1); p2 <- ncol(x2)
lambda1 <- 0.7
beta1   <- c(-0.2, 0.6)
lambda2 <- 0.9
beta2   <- c( 0.3, 0.5, -0.4)

para <- c(lambda1, beta1, lambda2, beta2)

sim <- rHZIP(n, m, para, x1, x2)
head(sim)
attr(sim, "propZeros")

Residuals for HZIP Models

Description

Computes randomized quantile residuals for fitted HZIP models.

Usage

## S3 method for class 'HZIP'
residuals(object, method = c("ghq", "laplace"), Q = 21, ...)

Arguments

Details

The residuals are obtained by integrating out the latent random effects using either Gauss-Hermite quadrature or a Laplace approximation.

Let Y_{ij} denote the response variable for the j-th observation of the i-th subject. The residuals are computed using the randomized quantile residual approach proposed by Dunn and Smyth (1996).

For each observation, the conditional cumulative distribution function is evaluated as

F^*(y_{ij}) = F(y_{ij} - 1) + U_{ij} P(Y_{ij} = y_{ij}),

where U_{ij} \sim \mathrm{Uniform}(0,1).

The residuals are then obtained as

r_{ij} = \Phi^{-1}(F^*(y_{ij})),

where \Phi^{-1}(\cdot) is the standard normal quantile function.

Random effects are approximated through either:

• Gauss-Hermite quadrature (method = "ghq");

• Laplace approximation (method = "laplace").

The Laplace approximation computes posterior modes using optim, while the quadrature approach relies on gauss.quad.

Value

A numeric vector containing the randomized quantile residuals.

References

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

Examples

fit.salamander <- hzip(
  y ~ mined | mined + spp,
  data = salamanders
)

res1 <- residuals(fit.salamander, method = "ghq")
res2 <- residuals(fit.salamander, method = "laplace")

head(res1)
head(res2)

Salamanders data

Description

This dataset is adapted from the glmmTMB package and contains salamander counts with information on mining status and species. It is intended for illustrating zero-inflated Poisson models with random effects using the hzip() function.

Usage

salamanders

Format

A data frame with 644 rows and 4 variables:

Ind

Individual identifier (factor).

y

Count response variable (integer).

mined

Mining status: "yes" or "no".

spp

Species factor with multiple levels (e.g., GP, PR, DM, etc.).

Details

The dataset was originally included in the glmmTMB package (Brooks et al., 2017), and has been slightly modified for testing the HZIP package.

Source

Adapted from the glmmTMB package.

Examples

data(salamanders, package = "HZIP")

## Fit zero-inflated Poisson with random effects
fit.salamander <- hzip(y ~ mined | mined + spp + mined * spp,
                       data = salamanders)
summary(fit.salamander)

Simulate Responses from a Fitted HZIP Model

Description

Generates nsim replicated response vectors from a fitted HZIP model by sampling latent random effects from the Generalized Log-Gamma distribution and then drawing from the hurdle-Poisson component.

Usage

## S3 method for class 'HZIP'
simulate(object, nsim = 1, seed = NULL, ...)

Arguments

Details

For each replicate and each observation i, the algorithm:

1. Samples b_{1i} \sim \mathrm{LGG}(0, \sigma_1, \lambda_1) and b_{2i} \sim \mathrm{LGG}(0, \sigma_2, \lambda_2) independently via rgengamma().

2. Computes \pi_i = 1 - \exp\{-\exp(\eta_i + b_{1i})\} and \mu_i = \exp(\rho_i + b_{2i}).

3. Draws Y_i: with probability \pi_i sets Y_i = 0; otherwise draws from a zero-truncated \mathrm{Poisson}(\mu_i) via rejection sampling.

Value

A data.frame with n rows and nsim columns, where n is the number of observations in object. Each column is one simulated response vector, named sim_1, sim_2, ..., sim_nsim. This follows the convention of the S3 generic simulate.

See Also

hzip, testDisp.HZIP, testZI.HZIP

Examples

fit.salamander <- hzip(
  y ~ mined | mined + spp,
  data = salamanders
)

sims <- simulate(fit.salamander, nsim = 10, seed = 42)
dim(sims)   # n x 10
head(sims)

Test for Overdispersion

Description

Generic function for simulation-based dispersion tests.

Usage

testDisp(object, ...)

Arguments

Simulation-Based Dispersion Test for HZIP Models

Description

Compares the observed variance of the response with the distribution of variances obtained from parametric bootstrap replicates simulated under the fitted HZIP model.

Usage

## S3 method for class 'HZIP'
testDisp(
  object,
  nsim = 500,
  alternative = c("two.sided", "greater", "less"),
  alpha = 0.05,
  seed = NULL,
  plot = TRUE,
  ...
)

Arguments

Value

An object of class "hzip_test" (invisibly), a list with components:

statistic

Observed variance of the response.

p.value

Simulation-based p-value.

alternative

The alternative hypothesis used.

alpha

Significance level used for the decision.

simulated

Numeric vector of simulated variances.

method

Description of the test.

h0

Statement of the null hypothesis.

See Also

testZI.HZIP, simulate.HZIP

Examples

fit.salamander <- hzip(
  y ~ mined | mined + spp,
  data = salamanders
)

testDisp(fit.salamander, nsim = 200, seed = 42)

Test for Zero-Inflation

Description

Generic function for simulation-based zero-inflation tests.

Usage

testZI(object, ...)

Arguments

Simulation-Based Zero-Inflation Test for HZIP Models

Description

Compares the observed proportion of zeros with the distribution of zero proportions obtained from parametric bootstrap replicates simulated under the fitted HZIP model.

Usage

## S3 method for class 'HZIP'
testZI(
  object,
  nsim = 500,
  alternative = c("two.sided", "greater", "less"),
  alpha = 0.05,
  seed = NULL,
  plot = TRUE,
  ...
)

Arguments

Value

An object of class "hzip_test" (invisibly), a list with components:

statistic

Observed proportion of zeros.

p.value

Simulation-based p-value.

alternative

The alternative hypothesis used.

alpha

Significance level used for the decision.

simulated

Numeric vector of simulated zero proportions.

method

Description of the test.

h0

Statement of the null hypothesis.

See Also

testDisp.HZIP, simulate.HZIP

Examples

fit.salamander <- hzip(
  y ~ mined | mined + spp,
  data = salamanders
)

testZI(fit.salamander, nsim = 200, seed = 42)