Counts can exhibit variation beyond that allowed by a Poisson model, and the process generating a zero may differ from the process determining a positive count. These features motivate negative-binomial, hurdle, and zero-inflated models. We describe their mean structures and use predictions to distinguish the quantities they estimate.
library(fastglm)
set.seed(544)
d <- data.frame(x=rnorm(500), exposure=runif(500,.8,1.2))
d$y <- rnbinom(500,mu=d$exposure*exp(.3+.5*d$x),size=1.5)The negative-binomial model used here, often called NB2, has variance
\(\mu+\mu^2/\theta\), where larger
\(\theta\) corresponds to less
overdispersion. fastglm_nb() alternates between estimating
the regression coefficients with \(\theta\) held fixed and estimating \(\theta\) with the fitted means held fixed.
A fixed value of \(\theta\) can instead
be supplied through negbin(theta).
X <- model.matrix(~x,d)
nb <- fastglm_nb(X,d$y,offset=log(d$exposure))
c(theta=nb$theta,SE.theta=nb$SE.theta)
#> theta SE.theta
#> 1.5243381 0.2398187
summary(nb)
#> Negative-binomial regression
#>
#> Call:
#> fastglm_nb(x = X, y = d$y, offset = log(d$exposure))
#> Family: Negative Binomial(1.5243) (log link)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 0.25871 0.05441 4.755 <1e-04 ***
#> x 0.37117 0.05316 6.983 <1e-04 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Dispersion Null deviance Null df Deviance Residual df AIC
#> 1 580.8 499 530.6 498 1570
#> Theta Std. Error
#> 1.524 0.2398
#> Convergence: converged (10 iterations)For a fixed-theta family, the negative-binomial likelihood is fully
specified and its scale parameter is one. When theta is estimated,
SE.theta accounts for estimating the coefficients at the
same time. The default coefficient covariance retains the conventional
GLM calculation with theta fixed; vcov(nb, full=TRUE)
instead returns the joint covariance obtained from the observed
likelihood curvature for coefficients and theta. The estimated theta in
fastglm_nb() contributes an additional parameter to AIC.
fastglm_confint_theta() obtains an interval by holding
theta at a sequence of candidate values and refitting the coefficients.
The resulting profile likelihood measures how well each candidate fits
after allowing the coefficients to adjust. The interval retains theta
values that are not rejected by a likelihood-ratio comparison. Its upper
endpoint can be infinite when the data remain compatible with the
Poisson limit, in which theta grows without bound.
fastglm_confint_theta(nb)
#> lower upper
#> 1.136223 2.114785
#> attr(,"level")
#> [1] 0.95
#> attr(,"profile")
#> theta logLik
#> 1 1.524338 -782.0222
#> 2 0.762169 -793.9905
#> 3 0.762169 -793.9905
#> 4 1.524338 -782.0222
#> 5 1.363864 -782.2817
#> 6 1.019556 -785.7401
#> 7 1.185980 -783.4048
#> 8 1.126418 -784.0633
#> 9 1.137083 -783.9326
#> 10 1.136237 -783.9427
#> 11 1.136223 -783.9429
#> 12 1.136229 -783.9428
#> 13 1.136223 -783.9429
#> 14 1.524338 -782.0222
#> 15 3.048676 -789.4464
#> 16 1.524338 -782.0222
#> 17 3.048676 -789.4464
#> 18 1.823730 -782.6314
#> 19 2.357957 -785.2901
#> 20 2.070129 -783.7144
#> 21 2.117989 -783.9596
#> 22 2.114698 -783.9424
#> 23 2.114785 -783.9429
#> 24 2.114795 -783.9429
#> 25 2.114785 -783.9429A zero-truncated model describes counts conditional on being positive. If the underlying count distribution has mean \(\mu\) and zero probability \(p_0\), its conditional mean is \(\mu/(1-p_0)\). The regression coefficients describe \(\log(\mu)\), whereas response predictions report the conditional mean.
positive <- subset(d,y>0)
truncated <- fastglm_truncated(y~x+offset(log(exposure)),data=positive,
dist="negbin")
head(predict(truncated,se.fit=TRUE)$fit)
#> 1 3 4 6 7 13
#> 1.916825 2.666691 3.908014 1.983223 3.253662 2.363961A hurdle model additionally estimates \(\pi=\Pr(Y>0)\), giving marginal mean \(\pi\mu/(1-p_0)\). Its likelihood separates into a binary component and a positive-count component. Within the positive-count component, the regression coefficients and theta are estimated from the same outcomes and their estimation errors can be related. The covariance calculation therefore uses the curvature of their joint log likelihood, called the observed information. This accounts for uncertainty in theta when reporting coefficient standard errors.
hurdle <- fastglm_hurdle(y~x+offset(log(exposure))|x,data=d,dist="negbin")
summary(hurdle)
#> Hurdle regression
#>
#> Call:
#> fastglm_hurdle(formula = y ~ x + offset(log(exposure)) | x, data = d,
#> dist = "negbin")
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> count_(Intercept) 0.30076 0.11648 2.582 0.00982 **
#> count_x 0.36578 0.07098 5.153 < 1e-04 ***
#> zero_(Intercept) 0.42498 0.09362 4.540 < 1e-04 ***
#> zero_x 0.41755 0.09699 4.305 < 1e-04 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Theta Std. Error
#> 1.746 0.6344
#> Log-likelihood: -782.1
#> Observations: 500
#> Convergence: converged
predict(hurdle,newdata=d[1:4,],type="response",se.fit=TRUE)
#> Predictions with standard errors
#> Scale: response mean
#>
#> Row Estimate Std. Error
#> 1 1.1480 0.06080
#> 2 0.4857 0.08546
#> 3 1.7494 0.10545
#> 4 3.0366 0.35039
predict(hurdle,newdata=d[1:4,],type="zero")
#> 1 2 3 4
#> 0.4010736 0.6524941 0.3439664 0.2229778
predict(hurdle,newdata=d[1:4,],type="count")
#> 1 2 3 4
#> 1.916825 1.397675 2.666691 3.908014The first right-hand side specifies the count model, and the second specifies the binary model. An offset inside a right-hand side applies only to that component. An offset supplied as a vector through the function argument applies to the count component; a two-column matrix supplies count and binary offsets separately.
vcov(fit) has one row per regression coefficient.
vcov(fit, full=TRUE) also includes theta when it is
estimated. Prediction standard errors use this full covariance to
account for how a small change in each fitted parameter changes the
prediction. This local approximation is called the delta method; it
includes dependence between the regression coefficients and theta.
For these count models, the joint information is computed from analytic likelihood derivatives and accumulated over blocks of observations. This avoids repeatedly perturbing each parameter and recalculating all score contributions to approximate the curvature. The calculation is useful when many coefficients are fitted or observations are read in blocks, while retaining the same joint uncertainty target. Numerical derivatives remain useful as an independent check of the implementation.
A zero-inflated model assigns probability \(\pi\) to a structural zero and probability \(1-\pi\) to an ordinary count distribution. Its marginal zero probability is \(\pi+(1-\pi)p_0\), and its marginal mean is \((1-\pi)\mu\). For an observed zero, the data do not reveal which process generated it. Estimation therefore involves uncertainty about both components, and their coefficient errors can be related. The covariance calculation includes this dependence when reporting uncertainty in the overall mean or zero probability.
set.seed(545)
d$inflated <- ifelse(rbinom(nrow(d),1,plogis(-.8+.4*d$x)),0,d$y)
zi <- fastglm_zi(inflated~x+offset(log(exposure))|x,data=d,
dist="poisson",em.maxit=300)
predict(zi,newdata=d[1:4,],type="prob",at=0:4)
#> 0 1 2 3 4
#> 1 0.5984625 0.18071284 0.12915008 0.06153314 0.021987951
#> 2 0.7316940 0.20074527 0.05566819 0.01029147 0.001426952
#> 3 0.5339135 0.09641854 0.12592323 0.10963770 0.071593771
#> 4 0.5619508 0.01746113 0.04189834 0.06702397 0.080412716
predict(zi,newdata=d[1:4,],type="response",se.fit=TRUE)
#> Predictions with standard errors
#> Scale: response mean
#>
#> Row Estimate Std. Error
#> 1 0.7546 0.04888
#> 2 0.3496 0.05674
#> 3 1.3138 0.09052
#> 4 2.1197 0.30496The same type="zero" prediction reports the marginal
probability of observing zero in both hurdle and zero-inflated models.
The meaning of type="count" differs: it is the conditional
positive mean for a hurdle model and the ordinary count-component mean
for a zero-inflated model.
A predictor can contribute to the count regression, the binary
regression, or both. fastglm_wald_test() tests the
corresponding coefficients jointly and prints the hypothesis, statistic,
degrees of freedom, and p-value. By default, a term is tested in each
component where it appears. Selecting component="count" or
component="zero" restricts the hypothesis to that
component. For a categorical predictor, all coefficients belonging to
the formula term are included.
fastglm_wald_test(hurdle,terms="x")
#> Joint Wald test
#> Terms: x
#> Estimation: fitted model coefficients
#> Covariance: model-based
#> Component: Count distribution, Positive-outcome probability
#>
#> Statistic df p-value
#> 45.09 2 < 1e-04
#>
#> Tested contrasts:
#> Estimate Std. Error Null
#> count_x 0.36578 0.07098 0.00000
#> zero_x 0.41755 0.09699 0.00000
#>
#> Note: First-order chi-square approximation; validity depends on coefficient
#> estimability and the covariance assumptions.
fastglm_wald_test(hurdle,terms="x",component="count")
#> Joint Wald test
#> Terms: x
#> Estimation: fitted model coefficients
#> Covariance: model-based
#> Component: Count distribution
#>
#> Statistic df p-value
#> 26.56 1 < 1e-04
#>
#> Tested contrasts:
#> Estimate Std. Error Null
#> count_x 0.36578 0.07098 0.00000
#>
#> Note: First-order chi-square approximation; validity depends on coefficient
#> estimability and the covariance assumptions.
fastglm_wald_test(zi,terms="x",component="zero")
#> Joint Wald test
#> Terms: x
#> Estimation: fitted model coefficients
#> Covariance: model-based
#> Component: Structural-zero probability
#>
#> Statistic df p-value
#> 1.947 1 0.1629
#>
#> Tested contrasts:
#> Estimate Std. Error Null
#> zero_x 0.1706 0.1223 0.0000
#>
#> Note: First-order chi-square approximation; validity depends on coefficient
#> estimability and the covariance assumptions.
fastglm_wald_test(truncated,terms="x")
#> Joint Wald test
#> Terms: x
#> Estimation: fitted model coefficients
#> Covariance: model-based
#>
#> Statistic df p-value
#> 26.56 1 < 1e-04
#>
#> Tested contrasts:
#> Estimate Std. Error Null
#> x 0.36578 0.07098 0.00000
#>
#> Note: First-order chi-square approximation; validity depends on coefficient
#> estimability and the covariance assumptions.The binary coefficient describes positive-outcome probability in a hurdle model and structural-zero probability in a zero-inflated model. The printed component label makes this distinction explicit. These are first-order chi-square tests, so sparse outcomes, extreme coefficients, or weakly identified mixture components can make their approximation inaccurate. An estimated negative-binomial size at its numerical boundary returns an unavailable test with an explanation.
Named contrasts can also combine coefficients from different
components using the count_ and zero_ names
returned by coef(). Such a hypothesis needs a substantive
interpretation: equality between a log-mean coefficient and a log-odds
coefficient does not mean equal changes in the outcome. For comparisons
on the outcome scale, the average predictions below are often easier to
interpret. The test retains cross-component covariance, including
covariance supplied through fastglm_vcov().
For truncated, hurdle, and zero-inflated negative-binomial fits, the
default regression covariance accounts for estimating size jointly.
fastglm_nb() retains its conventional GLM covariance with
size fixed. Its joint observed covariance can instead be used by
selecting the regression rows and columns:
regression <- names(coef(nb))
joint_covariance <- vcov(nb,full=TRUE)[regression,regression,drop=FALSE]
fastglm_wald_test(nb,terms="x",vcov.=joint_covariance)
#> Joint Wald test
#> Terms: x
#> Estimation: fitted model coefficients
#> Covariance: supplied coefficient covariance
#>
#> Statistic df p-value
#> 47.81 1 < 1e-04
#>
#> Tested contrasts:
#> Estimate Std. Error Null
#> x 0.37117 0.05368 0.00000
#>
#> Note: First-order chi-square approximation; validity depends on coefficient
#> estimability and the covariance assumptions.Simulation from a fitted model can help assess whether its zero frequency and upper tail resemble those in the observed data. This comparison describes model adequacy and does not replace uncertainty assessment or validation on new observations.
replicates <- simulate(hurdle,nsim=3,seed=546)
c(observed=mean(d$y==0),vapply(replicates,function(z) mean(z==0),numeric(1)))
#> observed sim_1 sim_2 sim_3
#> 0.402 0.408 0.414 0.380For nonnegative continuous outcomes, the two-part model in the semicontinuous-outcomes vignette uses a binary component together with a positive Gamma distribution.
A large count analysis may require predictors to be stored sparsely
or read in blocks. The model specification and uncertainty target should
remain the same when the data representation changes.
fastglm_nb() accepts a sparse matrix directly, while
fastglm_hurdle() accepts separate count and binary
predictor matrices through its first argument and z. The
response for the matrix hurdle interface is supplied as
response.
Xs <- as(model.matrix(~x,d),"dgCMatrix")
nb_sparse <- fastglm_nb(Xs,d$y,offset=log(d$exposure))
hurdle_sparse <- fastglm_hurdle(Xs,response=d$y,z=Xs,
offset=log(d$exposure),dist="negbin")
max(abs(unname(coef(nb_sparse))-unname(coef(nb))))
#> [1] 7.545836e-11When data are read in blocks, fastglm_source() records a
function that returns the requested rows and the formula information
needed to construct the same predictor columns in each block. A formula
with two right-hand sides specifies separate count and binary
regressions. The next fit reads five blocks of 100 observations, and its
coefficient summary requires no further reads.
read_block <- function(k) d[((k-1)*100+1):(k*100),]
count_source <- fastglm_source(read_block,5,
formula=y~x+offset(log(exposure))|x,prototype=d)
hurdle_stream <- fastglm_hurdle(count_source,dist="negbin")
max(abs(unname(coef(hurdle_stream))-unname(coef(hurdle))))
#> [1] 4.748764e-08
summary(hurdle_stream)
#> Hurdle regression
#>
#> Call:
#> fastglm_hurdle(formula = count_source, dist = "negbin")
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> count_(Intercept) 0.30076 0.11648 2.582 0.00982 **
#> count_x 0.36578 0.07098 5.153 < 1e-04 ***
#> zero_(Intercept) 0.42498 0.09362 4.540 < 1e-04 ***
#> zero_x 0.41755 0.09699 4.305 < 1e-04 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Theta Std. Error
#> 1.746 0.6344
#> Log-likelihood: -782.1
#> Observations: 500
#> Convergence: convergedWhile fitting, these models hold one block of predictor rows and coefficient-sized matrices in memory, and fitted objects keep small records of the block sizes. The source must return the same observations throughout fitting and later inference. The large-data vignette explains how to use file or database readers and how to save predictions block by block when the full prediction vector would be too large for memory.
An overall average count can be more directly interpretable than the
coefficients of the two regressions. fastglm_average()
averages the fitted response means over the supplied predictor
distribution. Its standard error combines uncertainty in both components
and, for negative-binomial counts, theta. The calculation conditions on
the predictor values used for averaging.
If observations are correlated within independent groups, a cluster-robust covariance can be used. The estimating contributions from both components are summed within each group before uncertainty is calculated; estimating the two components separately does not justify omitting their cross-component covariance under this sampling structure.
cluster <- rep(seq_len(100),each=5)
V <- fastglm_vcov(hurdle,type="HC1",cluster=cluster,full=TRUE)
fastglm_average(hurdle,vcov.=V)
#> estimate std.error conf.low conf.high
#> 1 1.368745 0.07139844 1.228807 1.508684Here the group labels illustrate the interface. In an analysis, they
should identify the independent sampling groups, such as sites or
patients, and a sufficiently large number of groups is needed for the
covariance approximation. full=TRUE includes theta so that
uncertainty in the average response can use the complete parameter
covariance.
A bootstrap provides another way to assess uncertainty by repeatedly resampling the independent units and refitting the model. The cluster bootstrap below samples whole groups, preserving observations within each selected group. Its statistic is the average prediction on the original predictor distribution, so the interpretation matches the preceding conditional average.
boot <- fastglm_bootstrap(hurdle,R=49,seed=547,type="cluster",cluster=cluster,
statistic=function(model) mean(predict(model,newdata=d,type="response")))
#> Warning: 1 bootstrap refits failed; uncertainty summaries are unavailable
boot
#> Refit bootstrap (cluster)
#> Refits: 48 of 49 successful
#>
#> Estimates and confidence intervals:
#> Estimate Std. Error lower upper
#> statistic1 1.369 NA NA NA
#> Inference: Unavailable; inspect the refit diagnosticsThe small number of repetitions keeps this example short; a
substantive analysis needs enough repetitions to estimate the interval
endpoints adequately. The reported percentile interval uses quantiles of
the refitted statistics. Failed or nonconverged refits are recorded in
diagnostics, and uncertainty summaries are withheld if any
refit fails. The procedure keeps the original predictor coding and
transformations fixed, so it does not repeat model selection or
data-dependent preprocessing.