Poisson Regression Inference for Count Responses
Source:R/inference_count_poisson.R
InferenceCountPoisson.RdFits a Poisson log-link regression for count responses:
\(Y_i \mid w_i, x_i \sim \mathrm{Poisson}(\mu_i)\), \(\log \mu_i =
\beta_0 + \beta_T w_i + x_i^\top \gamma\)
(fast_poisson_regression_cpp/
fast_poisson_regression_with_var_cpp). \(\hat\beta_T\) is a
log-rate-ratio: \(\exp(\hat\beta_T)\) is the estimated rate ratio.
likelihood_tier = "full": Wald, score, gradient, and
likelihood-ratio tests are all available, plus parametric-likelihood
bootstrap calibration of the likelihood-ratio test (simulating new Poisson
responses under the null).
Design-conservative testing. Every asymptotic/likelihood test
method on this class does not report the raw model-based result directly.
Instead, it also computes a design-based jackknife-Wald test
(compute_jackknife_wald_two_sided_pval()/
compute_jackknife_wald_confidence_interval(), which do not assume
the Poisson mean-variance relationship) and combines the two
conservatively: p-values report \(\max\) of the model-based and
design-based p-values, and confidence intervals report the union of
the model-based and design-based intervals. This guards against the
model-based test being anti-conservative when the Poisson equidispersion
assumption (\(\mathrm{Var}(Y_i) = E[Y_i]\)) fails — a real risk for count
data, which is frequently overdispersed (see
InferenceCountNegBin for a model
that estimates dispersion directly instead). If either component is
unavailable, the available one is used alone; if neither is available, the
result is NA. Validity requires the usual correctly-specified
linear predictor on the log scale; unlike the raw Poisson likelihood
alone, this class's actual reported inference degrades gracefully (rather
than becoming anti-conservative) under mean-variance misspecification.
Estimand. Composes
MarginalEstimand
(set_estimand()/get_estimand()/get_supported_estimands()).
Under the default estimand = "conditional", \(\hat\beta_T\) is the
log-rate-ratio above. Under estimand = "marginal_mean_diff", the
reported quantity is instead the g-computation marginal rate difference
\(\frac{1}{n}\sum_i \{\exp(\hat\beta_0 + \hat\beta_T + X_i^\top
\hat\gamma) - \exp(\hat\beta_0 + X_i^\top \hat\gamma)\}\). Under
estimand = "marginal_ratio", the log of the corresponding marginal
rate ratio — which is numerically identical to the conditional
\(\hat\beta_T\) for this family: because the log link is linear in
\(w_i\) with no treatment-by-covariate interaction term, every subject's
treated-vs-control mean ratio is \(\exp(\hat\beta_0 + \hat\beta_T +
X_i^\top\hat\gamma) / \exp(\hat\beta_0 + X_i^\top\hat\gamma) =
\exp(\hat\beta_T)\) exactly, so averaging over subjects before or after
taking the ratio makes no difference. "marginal_ratio" is provided
for estimand-API consistency with the other model families, not because it
differs numerically from "conditional" here; "marginal_mean_diff"
is the estimand where g-computation actually changes the reported number
for a Poisson GLM, since a difference (unlike a ratio) does not collapse
under a nonlinear (log) mean function. Because there is no latent submodel
for this family (unlike e.g.
InferenceCountZeroInflatedPoisson's
excess-zero mixture), the marginal mean function is exactly the model's
own fitted mean; no separate standardization step beyond the g-computation
average is needed. Standard errors under a marginal estimand use the delta
method against the model's coefficient covariance (degrees of freedom
Inf), including in the design-conservative union/max combination
above (the design-based jackknife-Wald component also refits under the
active estimand); testing_type is restricted to "wald"
whenever the estimand is non-conditional. The underlying model fit is
identical regardless of estimand — switching estimand is a pure
post-fit transform, never a refit.
References
Cameron, A. C., and Trivedi, P. K. (2013). Regression Analysis of Count Data (2nd ed.). Cambridge University Press, for the Poisson regression model and its maximum-likelihood theory.
See also
Comparable Python API:
statsmodels
discrete models (Poisson). See also:
Poisson
regression (Wikipedia).
Super class
Inference -> InferenceCountPoisson
Methods
Public methods
InferenceCountPoisson$compute_lik_ratio_confidence_interval()InferenceCountPoisson$compute_gradient_confidence_interval()InferenceCountPoisson$compute_lik_ratio_bootstrap_two_sided_pval()InferenceCountPoisson$compute_lik_ratio_bootstrap_confidence_interval()InferenceCountPoisson$compute_estimate_with_bootstrap_weights()
+ inherited public methods from Inference
Inference$capabilities()Inference$compute_exact_confidence_interval()Inference$compute_exact_two_sided_pval_for_treatment_effect()Inference$duplicate()Inference$get_analysis_data()Inference$get_covariates()Inference$get_design_object()Inference$get_model_formula()Inference$get_nonestimable_reason()Inference$get_nonestimable_stage()Inference$get_optimization_alg()Inference$get_response()Inference$get_response_type()Inference$get_treatment()Inference$is_nonestimable()Inference$set_optimization_alg()Inference$set_seed()Inference$supports()
InferenceCountPoisson$new()
Initialize inference for the Poisson regression model
\(Y_i \mid w_i, x_i \sim \mathrm{Poisson}(\mu_i)\), \(\log \mu_i =
\beta_0 + \beta_T w_i + x_i^\top \gamma\); see
InferenceCountPoisson for the
model form and the design-conservative testing mechanism. Does not fit
the model; the fit is deferred to the first call to
compute_estimate() or a method that requires it.
Usage
InferenceCountPoisson$new(
des_obj,
model_formula = NULL,
verbose = FALSE,
smart_cold_start_default = NULL,
harden = TRUE
)Arguments
des_objA completed
Designobject with a count response.model_formulaOptional formula for covariate adjustment. If
NULL(default), the formula from the design object is used and its pre-computed design matrix is reused. If a formula is provided, a new design matrix is constructed from the design's imputed covariates.verboseWhether to print progress messages.
smart_cold_start_defaultWhether to use smart cold start values by default.
hardenWhether to apply robustness measures.
InferenceCountPoisson$compute_estimate()
Fits the Poisson regression model by maximum likelihood.
Under the default estimand = "conditional", returns
\(\hat\beta_T\), the treatment log-rate-ratio. Under
estimand = "marginal_mean_diff"/"marginal_ratio" (set
via set_estimand()), returns the g-computation marginal
rate difference/log-rate-ratio instead — see the class-level
@details for the formula. The underlying model fit is
identical either way (a pure post-fit transform of the same
cached fit, no refit).
InferenceCountPoisson$compute_asymp_confidence_interval()
Design-conservative confidence interval for \(\beta_T\)
using whichever test type is configured
(private$testing_type: "wald", "score",
"gradient", or "lik_ratio"); see
InferenceCountPoisson for the
union-with-jackknife-Wald combination rule.
InferenceCountPoisson$compute_asymp_two_sided_pval()
Design-conservative two-sided p-value for \(H_0: \beta_T =
\code{delta}\) using whichever test type is configured
(private$testing_type); see
InferenceCountPoisson for the
max-with-jackknife-Wald combination rule.
InferenceCountPoisson$compute_wald_confidence_interval()
Wald confidence interval for \(\beta_T\) using the fitted
Poisson model's Fisher-information-based standard error, unioned with
the design-based jackknife-Wald interval; see
InferenceCountPoisson for the
combination rule and InferenceAsymp
for the underlying Wald contract.
InferenceCountPoisson$compute_wald_two_sided_pval()
Wald test of \(H_0: \beta_T = \code{delta}\) using the
fitted Poisson model's Fisher-information-based standard error, taking
the max with the design-based jackknife-Wald p-value; see
InferenceCountPoisson for the
combination rule.
InferenceCountPoisson$compute_score_confidence_interval()
Score-test confidence interval for \(\beta_T\) (inverting
the Poisson score test at each candidate null, no full re-fit needed at
the observed information), unioned with the design-based jackknife-Wald
interval; see InferenceCountPoisson
for the combination rule.
InferenceCountPoisson$compute_score_two_sided_pval()
Score test of \(H_0: \beta_T = \code{delta}\), taking the
max with the design-based jackknife-Wald p-value; see
InferenceCountPoisson for the
combination rule.
InferenceCountPoisson$compute_lik_ratio_confidence_interval()
Likelihood-ratio-test confidence interval for \(\beta_T\)
(test inversion, requiring a null refit at each candidate value),
unioned with the design-based jackknife-Wald interval; see
InferenceCountPoisson for the
combination rule.
InferenceCountPoisson$compute_lik_ratio_two_sided_pval()
Likelihood-ratio test of \(H_0: \beta_T = \code{delta}\),
taking the max with the design-based jackknife-Wald p-value; see
InferenceCountPoisson for the
combination rule.
InferenceCountPoisson$compute_gradient_confidence_interval()
Gradient-test confidence interval for \(\beta_T\) (a
score-test variant using the observed rather than expected
information), unioned with the design-based jackknife-Wald interval;
see InferenceCountPoisson for
the combination rule.
InferenceCountPoisson$compute_gradient_two_sided_pval()
Gradient test of \(H_0: \beta_T = \code{delta}\), taking
the max with the design-based jackknife-Wald p-value; see
InferenceCountPoisson for the
combination rule.
InferenceCountPoisson$compute_lik_ratio_bootstrap_two_sided_pval()
Parametric-likelihood-bootstrap-calibrated likelihood-ratio
test of \(H_0: \beta_T = \code{delta}\) (simulating new Poisson
responses from the null-constrained fit to calibrate the LR statistic's
null distribution), taking the max with the design-based jackknife-Wald
p-value; see InferenceCountPoisson
for the combination rule.
InferenceCountPoisson$compute_lik_ratio_bootstrap_confidence_interval()
Parametric-likelihood-bootstrap-calibrated likelihood-ratio
confidence interval for \(\beta_T\) (test inversion using the
bootstrap-calibrated null distribution), unioned with the design-based
jackknife-Wald interval; see
InferenceCountPoisson for the
combination rule.
Usage
InferenceCountPoisson$compute_lik_ratio_bootstrap_confidence_interval(
alpha = 0.05,
B = 199,
show_progress = FALSE,
min_number_usable_samples = 5L,
max_attempts_per_replicate = 2L,
root_tolerance = NULL,
max_root_iterations = 8L
)Arguments
alphaSignificance level. Default 0.05.
BNumber of bootstrap replicates.
show_progressWhether to show progress.
min_number_usable_samplesMinimum usable bootstrap samples.
max_attempts_per_replicateMaximum attempts per replicate.
root_toleranceRoot tolerance.
max_root_iterationsMaximum root iterations.
InferenceCountPoisson$compute_estimate_with_bootstrap_weights()
Refits the Poisson model with subject/block-level weights
applied to the fitting log-likelihood (Bayesian-bootstrap or
nonparametric-bootstrap draw weights, expanded to row level via
private$expand_subject_or_block_weights_to_row_weights()) via
fast_poisson_regression_weighted_cpp, and returns the
reweighted log-rate-ratio estimate \(\hat\beta_T^{(w)}\). Uses the
same QR column-dropping hardening as the unweighted fit; a hardened fit
with a non-finite treatment coefficient is cached as nonestimable and
returns NA.
Examples
# \donttest{
seq_des = DesignSeqOneByOneBernoulli$new(n = 10, response_type = 'count')
for (i in 1:10) {
seq_des$add_one_subject_to_experiment_and_assign(data.frame(x1 = rnorm(1)))
}
seq_des$add_all_subject_responses(rpois(10, 2))
inf = InferenceCountPoisson$new(seq_des)
inf$compute_estimate()
#> [1] -0.1362724
# }
# \donttest{
inf$set_seed(1)
inf$compute_lik_ratio_bootstrap_two_sided_pval(delta = 0, B = 9, show_progress = FALSE)
#> [1] 0.8911806
# }