Cookbook: Count Outcome, End to End

The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
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[@]dogado.de.

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.
Health stats visible at Monitor.