Health care expenditures and other nonnegative outcomes often contain a substantial proportion of zeros together with a right-skewed positive distribution. A single Gamma regression excludes the zero observations and therefore addresses only the positive-outcome mean. A two-part model describes both the probability of a positive outcome and its conditional mean, allowing inference about the overall mean.
Let \(\pi(z)=\Pr(Y>0\mid z)\) and \(m(x)=E(Y\mid Y>0,x)\). We use a binomial regression for \(\pi(z)\) and a Gamma regression with a log link for \(m(x)\). The resulting marginal mean is \[ E(Y\mid x,z)=\pi(z)m(x). \] The two regressions describe different quantities: whether an outcome is positive and how large it is when positive. Under the stated likelihood model, their coefficient uncertainty separates into two covariance blocks, with zero covariance between the blocks. Uncertainty in the overall mean still depends on both regressions. The delta method approximates this uncertainty by calculating how small changes in either set of coefficients change the product \(\pi(z)m(x)\).
library(fastglm)
set.seed(619)
d <- data.frame(age=rnorm(800,50,10), severity=rnorm(800))
probability <- plogis(-.4+.02*(d$age-50)+.5*d$severity)
positive_mean <- exp(5+.015*(d$age-50)+.4*d$severity)
d$cost <- rbinom(800,1,probability)*rgamma(800,shape=2,scale=positive_mean/2)
fit <- fastglm_twopart(cost~age+severity|age+severity,data=d)
summary(fit)
#> Two-part binomial / Gamma regression
#>
#> Call:
#> fastglm_twopart(formula = cost ~ age + severity | age + severity,
#> data = d)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> count_(Intercept) 4.220789 0.209400 20.157 < 1e-04 ***
#> count_age 0.015317 0.003991 3.838 0.000124 ***
#> count_severity 0.412469 0.038743 10.646 < 1e-04 ***
#> zero_(Intercept) -1.713377 0.375465 -4.563 < 1e-04 ***
#> zero_age 0.024144 0.007305 3.305 0.000950 ***
#> zero_severity 0.515819 0.077785 6.631 < 1e-04 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Log-likelihood: -2313
#> Observations: 800
#> Convergence: convergedThe first right-hand side specifies the positive-outcome regression, and the second specifies the probability of a positive outcome. A formula with a single right-hand side uses the same predictors in both regressions. The two regressions can therefore reflect different mechanisms when scientific knowledge supports that specification.
An association can be tested within either component or jointly
across both. The positive component describes the amount
among positive outcomes, whereas the binary component
describes the probability of a positive outcome. The default tests the
selected term wherever it occurs in the two regressions.
fastglm_wald_test(fit,terms="severity",component="positive")
#> Joint Wald test
#> Terms: severity
#> Estimation: fitted model coefficients
#> Covariance: model-based
#> Component: Positive amount
#>
#> Statistic df p-value
#> 113.3 1 < 1e-04
#>
#> Tested contrasts:
#> Estimate Std. Error Null
#> count_severity 0.41247 0.03874 0.00000
#>
#> Note: First-order chi-square approximation; validity depends on coefficient
#> estimability and the covariance assumptions.
fastglm_wald_test(fit,terms="severity",component="binary")
#> Joint Wald test
#> Terms: severity
#> Estimation: fitted model coefficients
#> Covariance: model-based
#> Component: Positive-outcome probability
#>
#> Statistic df p-value
#> 43.97 1 < 1e-04
#>
#> Tested contrasts:
#> Estimate Std. Error Null
#> zero_severity 0.51582 0.07779 0.00000
#>
#> Note: First-order chi-square approximation; validity depends on coefficient
#> estimability and the covariance assumptions.
fastglm_wald_test(fit,terms="severity")
#> Joint Wald test
#> Terms: severity
#> Estimation: fitted model coefficients
#> Covariance: model-based
#> Component: Positive amount, Positive-outcome probability
#>
#> Statistic df p-value
#> 157.3 2 < 1e-04
#>
#> Tested contrasts:
#> Estimate Std. Error Null
#> count_severity 0.41247 0.03874 0.00000
#> zero_severity 0.51582 0.07779 0.00000
#>
#> Note: First-order chi-square approximation; validity depends on coefficient
#> estimability and the covariance assumptions.These Wald tests use a first-order chi-square approximation. They assess the specified regression coefficients and do not test whether the entire two-part distribution is appropriate. A factor with several levels is tested as a complete formula term, avoiding separate tests of each reference-level comparison.
A coefficient in the positive-outcome regression describes a conditional association among observations with positive outcomes. Its effect on the overall mean also depends on the binary regression. Predictions make this distinction explicit.
new <- data.frame(age=c(40,60),severity=c(0,0))
cbind(probability_positive=1-predict(fit,new,type="zero"),
positive_mean=predict(fit,new,type="count"),
marginal_mean=predict(fit,new,type="response"))
#> probability_positive positive_mean marginal_mean
#> 1 0.3213365 125.6466 40.37485
#> 2 0.4341924 170.6843 74.10982
predict(fit,new,type="response",se.fit=TRUE)
#> Predictions with standard errors
#> Scale: response mean
#>
#> Row Estimate Std. Error
#> 1 40.37 3.844
#> 2 74.11 5.885The standard errors condition on the supplied covariates and use the model-based covariance of both regressions. They do not account for model selection or misspecification of the positive-outcome distribution. Simulation provides a complementary assessment of the fitted zero frequency and positive tail.
replicates <- simulate(fit,nsim=2,seed=620)
vapply(replicates,mean,numeric(1))
#> sim_1 sim_2
#> 64.91743 77.56146
mean(d$cost)
#> [1] 69.3556The firth=TRUE option modifies the estimating equations
to reduce the leading small-sample bias in both regressions. For the
Gamma regression, it also adjusts the equation used to estimate
dispersion, which determines the spread of positive outcomes around
their conditional mean. It can be useful when the positive subset is
small, although numerical convergence and the adequacy of both mean
models still require examination.
A change in a predictor can affect both the probability of a positive
outcome and its amount. An average prediction contrast combines these
changes on the outcome scale. The function supplied through
contrast below increases severity by one unit in both
regression designs and leaves the other predictors unchanged. The
reported estimate is the average change in the fitted overall mean.
change_severity <- function(design) {
design$count[,"severity"] <- design$count[,"severity"]+1
design$zero[,"severity"] <- design$zero[,"severity"]+1
design
}
fastglm_average(fit,contrast=change_severity)
#> estimate std.error conf.low conf.high
#> 1 63.9238 7.876662 48.48582 79.36177This comparison conditions on the observed predictor distribution and
describes the fitted association. A causal interpretation requires
additional assumptions about treatment assignment, confounding, and the
meaning of the intervention. If only the positive-outcome process is to
be changed, the transformation can instead modify
design$count alone.
Repeated outcomes within the same patient or site can make the model-based standard errors inappropriate. Cluster-robust covariance sums the estimating contributions from both regressions within each independent group. Their estimation errors may then be related across components even though the likelihood-based covariance has separate blocks.
cluster <- rep(seq_len(160),each=5)
V <- fastglm_vcov(fit,type="HC1",cluster=cluster,full=TRUE)
fastglm_average(fit,contrast=change_severity,vcov.=V)
#> estimate std.error conf.low conf.high
#> 1 63.9238 7.714562 48.80353 79.04406
fastglm_wald_test(fit,terms="severity",vcov.=V)
#> Joint Wald test
#> Terms: severity
#> Estimation: fitted model coefficients
#> Covariance: supplied coefficient covariance
#> Component: Positive amount, Positive-outcome probability
#>
#> Statistic df p-value
#> 195 2 < 1e-04
#>
#> Tested contrasts:
#> Estimate Std. Error Null
#> count_severity 0.41247 0.03199 0.00000
#> zero_severity 0.51582 0.07754 0.00000
#>
#> Note: First-order chi-square approximation; validity depends on coefficient
#> estimability and the covariance assumptions.The joint test uses the same cross-component covariance as the
average comparison. The group identifiers in this example are
illustrative. Actual identifiers should follow the sampling structure,
and inference requires enough independent groups for the approximation
to be useful. Robust covariance for mean-bias-reduced two-part fits is
not currently supported; their Wald tests use the model-based
first-order convention. A cluster bootstrap through
fastglm_bootstrap() provides a refitting-based alternative
under the same independence assumption; the count-model vignette gives a
worked example.
The same two-part model can be fitted when the complete predictor
matrix cannot be kept in memory. A reader supplies blocks of rows, and
fastglm_source() records the common factor coding and
transformations. Both component regressions then revisit the same data
as their coefficients are updated.
read_cost <- function(k) d[((k-1)*200+1):(k*200),]
cost_source <- fastglm_source(read_cost,4,
formula=cost~age+severity|age+severity,prototype=d)
stream_fit <- fastglm_twopart(cost_source)
max(abs(unname(coef(stream_fit))-unname(coef(fit))))
#> [1] 7.638334e-14
fastglm_average(stream_fit)
#> estimate std.error conf.low conf.high
#> 1 69.62791 4.260844 61.27681 77.97901The fit retains the reader and coefficient summaries, allowing
predictions and average outcomes to be calculated by another pass
through the blocks. This reduces design memory while requiring the
observations to remain accessible and unchanged. Sparse predictor
matrices are also supported through
fastglm_twopart(X, response=y, z=Z), where X
specifies the positive-outcome regression and Z the
probability of a positive outcome.