Cookbook: Count Outcome, End to End

The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

One complete, runnable script for a count outcome — number of events per subject (visits, relapses, defects). Same four steps as every cookbook (design → assign → record → infer; see vignette("cookbook-continuous")). The natural estimand is a log rate ratio; the inference classes are the InferenceCount* family, covering Poisson, negative binomial, hurdle and zero-inflated variants for over-dispersed or zero-heavy counts.

Setup

EDI is not on CRAN yet, so install.packages("EDI") fails — install from R-universe (fallback: GitHub, subdir = "R/EDI"). Not evaluated here.

install.packages("EDI", repos = c("https://kapelner.r-universe.dev", "https://cloud.r-project.org"))
# or: remotes::install_github("kapelner/EDI", subdir = "R/EDI")

library(EDI)
set.seed(20260916)

n = 80
X = data.frame(
  baseline_rate = round(rgamma(n, 4, 1), 1),
  urban         = rbinom(n, 1, 0.5)
)
true_log_rr = -0.4   # treatment reduces the event rate by ~33%

Fixed design, Poisson regression

des = DesignFixedBernoulli$new(n = n, response_type = "count", verbose = FALSE)
des$add_all_subjects_to_experiment(X)
des$assign_w_to_all_subjects()
w = des$get_w()

mu = exp(0.5 + true_log_rr * w + 0.15 * X$baseline_rate + 0.3 * X$urban)
y = rpois(n, mu)
des$add_all_subject_responses(y)

inf = InferenceCountPoisson$new(des, verbose = FALSE)
inf$num_cores = 1L
inf$compute_estimate()                         # log rate ratio for treatment
#> [1] -0.2878601
inf$compute_asymp_confidence_interval(alpha = 0.05)
#>         2.5%        97.5% 
#> -0.567025852 -0.008694318
inf$compute_asymp_two_sided_pval()
#> [1] 0.04327925

inf$set_seed(1)
inf$compute_rand_two_sided_pval(r = 200, show_progress = FALSE)
#> [1] 0.03
inf$set_seed(1)
inf$compute_bootstrap_confidence_interval(alpha = 0.05, B = 200, show_progress = FALSE)
#>        2.5%       97.5% 
#> -0.50039002 -0.04860621

Over-dispersion: negative binomial on the same design

Real counts are usually over-dispersed relative to Poisson. Swap the class; the design object is unchanged.

inf_nb = InferenceCountNegBin$new(des, verbose = FALSE)
inf_nb$num_cores = 1L
inf_nb$compute_estimate()
#> [1] -0.2876102
inf_nb$compute_asymp_confidence_interval(alpha = 0.05)
#>       2.5%      97.5% 
#> -0.5471982 -0.0280222

Everything at once

suite = InferenceSuite$new(des)
res = suite$run_all_inference(screen = TRUE, plots = FALSE, num_cores = 1L,
                              methods = c("wald", "score", "lik_ratio"), max_secs_per_class = 15)
#> inference       cov      estimand    est       se        pval       pval method     status 
#> class           mod                                                                        
#> ===========================================================================================
#> Classes 0/23  [                0%                 ] Status: Estimating...[KAvg Δ                    mean Δ      -0.742    0.449     1.02e-01   wald            ok     
#> Classes 1/23  [=              4%                ] Estimated Time Left: 0s[KAvg Δ Pooled …           mean Δ      -0.742    0.445     9.95e-02   wald            ok     
#> Classes 2/23  [==             8%                ] Estimated Time Left: 0s[KWilcox                   HL shift    -1.00     0.510     9.72e-02   wald            ok     
#> Classes 3/23  [====           13%               ] Estimated Time Left: 0s[KHurd Neg Bin    ~.       log rate …  -0.360    0.147     1.45e-02   wald            ok     
#> Classes 4/23  [=====          17%               ] Estimated Time Left: 0s[KHurd Neg Bin    ~.       log rate …  -0.360    0.147     1.38e-02   score           ok     
#> Classes 5/23  [=======        21%               ] Estimated Time Left: 0s[KHurd Neg Bin    ~.       log rate …  -0.360    0.147     1.38e-02   lik_ratio       ok     
#> Classes 6/23  [========       26%               ] Estimated Time Left: 0s[KHurd Poisson    ~.       log rate …  -0.360    0.133     6.96e-03   wald            ok     
#> Classes 7/23  [==========     30%               ] Estimated Time Left: 0s[KHurd Poisson    ~.       log rate …  -0.360    0.133     1.38e-02   score           ok     
#> Classes 8/23  [===========    34%               ] Estimated Time Left: 0s[KHurd Poisson    ~.       log rate …  -0.360    0.133     1.38e-02   lik_ratio       ok     
#> Classes 9/23  [============   39%               ] Estimated Time Left: 0s[KNeg Bin         ~.       log rate …  -0.288    0.132     2.99e-02   wald            ok     
#> Classes 10/23 [============== 43%               ] Estimated Time Left: 0s[KNeg Bin         ~.       log rate …  -0.288    0.132     2.92e-02   score           ok     
#> Classes 11/23 [============== 47%               ] Estimated Time Left: 0s[KNeg Bin         ~.       log rate …  -0.288    0.132     2.95e-02   lik_ratio       ok     
#> Classes 12/23 [============== 52%               ] Estimated Time Left: 0s[KPoisson         ~.       log rate …  -0.288    0.132     4.33e-02   wald            ok     
#> Classes 13/23 [============== 56%               ] Estimated Time Left: 0s[KPoisson         ~.       log rate …  -0.288    0.132     4.33e-02   score           ok     
#> Classes 14/23 [============== 60% =             ] Estimated Time Left: 0s[KPoisson         ~.       log rate …  -0.288    0.132     4.33e-02   lik_ratio       ok     
#> Classes 15/23 [============== 65% ==            ] Estimated Time Left: 0s[KQuasi Poisson   ~.       log rate …  -0.288    0.127     2.36e-02   wald            ok     
#> Classes 16/23 [============== 69% ===           ] Estimated Time Left: 0s[KRobust Poisson  ~.       log rate …  -0.288    0.119     1.52e-02   wald            ok     
#> Classes 17/23 [============== 73% =====         ] Estimated Time Left: 0s[KZero Infl Neg…  ~.       log rate …  -0.306    NA        NA         wald            ok     
#> Classes 18/23 [============== 78% ======        ] Estimated Time Left: 0s[KZero Infl Neg…  ~.       log rate …  -0.306    NA        7.03e-03   score           ok     
#> Classes 19/23 [============== 82% ========      ] Estimated Time Left: 0s[KZero Infl Neg…  ~.       log rate …  -0.306    NA        2.09e-02   lik_ratio       ok     
#> Classes 20/23 [============== 86% =========     ] Estimated Time Left: 0s[KZero Infl Poi…  ~.       log rate …  -0.306    0.136     2.38e-02   wald            ok     
#> Classes 21/23 [============== 91% ===========   ] Estimated Time Left: 0s[KZero Infl Poi…  ~.       log rate …  -0.306    0.136     2.28e-02   score           ok     
#> Classes 22/23 [============== 95% ============  ] Estimated Time Left: 0s[KZero Infl Poi…  ~.       log rate …  -0.306    0.136     3.26e-02   lik_ratio       ok     
#> Classes 23/23 [============= 100% ==============] Estimated Time Left: 0s[K-------------------------------------------------------------------------------------------
#> Status: Completed in 1s.
#> 
#>   Estimand: HL shift (1 inferences)                : p =     NA
#>   Estimand: log rate ratio cond (5 inferences)     : p = 0.0124
#>   Estimand: log rate ratio marginal (14 inferences): p = 0.0205
#>   Estimand: mean Δ (2 inferences)                  : p = 0.1009
#> 
#> Combined evidence against the sharp null across 4 estimands
#> (22 inferences, weighting = uniform within estimand):
#> p = 0.0268

Sequential design

The one-by-one API is identical across response types — only the generated response changes. Randomization inference is the tool for sequential designs whose assignments depend on earlier subjects.

des_seq = DesignSeqOneByOneKK14$new(n = n, response_type = "count", verbose = FALSE)
for (i in seq_len(n)) {
  w_i = des_seq$add_one_subject_to_experiment_and_assign(X[i, , drop = FALSE])
  mu_i = exp(0.5 + true_log_rr * w_i + 0.15 * X$baseline_rate[i] + 0.3 * X$urban[i])
  des_seq$add_one_subject_response(i, rpois(1, mu_i))
}

inf_seq = InferenceCountPoisson$new(des_seq, verbose = FALSE)
inf_seq$num_cores = 1L
inf_seq$compute_estimate()
#> [1] -0.1474377
inf_seq$set_seed(1)
inf_seq$compute_rand_two_sided_pval(r = 200, show_progress = FALSE)
#> [1] 0.14

Where to go next

Zero-heavy data: InferenceCountHurdlePoisson, InferenceCountHurdleNegBin, InferenceCountZeroInflatedPoisson, InferenceCountZeroInflatedNegBin.
Exposure/offset (rates per unit time) is a planned addition — see ROADMAP.md.
vignette("validation-evidence"): each kernel is checked against stats::glm, MASS::glm.nb, or pscl.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.