Cookbook: Continuous Outcome, End to End

This cookbook is part of EDI (Experimental Design and Inference), an R package for randomized experimental designs with inference matched to the design and response type.

One complete, runnable script: design a two-arm experiment with a continuous outcome, assign treatment, record responses, and run every matched inference procedure. Every chunk executes when the vignette is built; copy the whole thing into a session and it works. The other cookbooks (vignette("cookbook-incidence"), -count, -proportion, -survival, -ordinal) follow exactly this shape, differing only in the response and the inference class.

The pattern is always the same four steps:

construct a Design for n subjects and a response_type;
assign treatment — all at once (fixed designs) or one subject at a time as they arrive (sequential designs);
record responses;
construct an Inference class on the completed design and call its methods — or hand the design to InferenceSuite to run every applicable procedure at once.

Setup and a simulated population

EDI is on CRAN. Not evaluated here: the vignette builds inside an already-installed package.

install.packages("EDI")
# or the development version: remotes::install_github("kapelner/EDI", subdir = "R/EDI")

library(EDI)
#> Welcome to EDI v1.0.3
set.seed(20260916)

n = 60
X = data.frame(
  age  = round(rnorm(n, 45, 12)),
  bmi  = round(rnorm(n, 27, 4), 1),
  male = rbinom(n, 1, 0.5)
)
true_effect = 2.5

A fixed design: everyone is known up front

DesignFixedBernoulli assigns each subject by an independent coin flip — the simplest fixed design, and the one under which every inference method is valid without caveat. Subjects are added all at once, assigned all at once, and responses recorded all at once.

des = DesignFixedBernoulli$new(n = n, response_type = "continuous", verbose = FALSE)
des$add_all_subjects_to_experiment(X)
des$assign_w_to_all_subjects()
w = des$get_w()                                  # 0/1 treatment vector

y = 10 + true_effect * w + 0.1 * X$age - 0.3 * X$bmi + rnorm(n, sd = 3)
des$add_all_subject_responses(y)

Inference: covariate-adjusted OLS

InferenceContinOLS regresses the response on treatment and the covariates. Every inference class exposes the same core verbs: a point estimate, an asymptotic (Wald) interval and p-value, a randomization test that re-randomizes under the design’s own mechanism, and — for fixed designs — a nonparametric bootstrap.

inf = InferenceContinOLS$new(des, verbose = FALSE)
inf$num_cores = 1L

inf$compute_estimate()
#> [1] 3.296457
inf$compute_asymp_confidence_interval(alpha = 0.05)
#>     2.5%    97.5% 
#> 1.517762 5.075152
inf$compute_asymp_two_sided_pval()
#> [1] 0.000478099

The randomization test uses no distributional assumption: it re-draws treatment r times from the same design and compares the observed statistic to that reference distribution. set_seed() makes the result reproducible (see vignette("reproducibility")).

inf$set_seed(1)
inf$compute_rand_two_sided_pval(r = 200, show_progress = FALSE)
#> [1] 0.01

Bootstrap intervals resample subjects (the design’s own resampling structure — rows here; matched pairs for matched designs):

inf$set_seed(1)
inf$compute_bootstrap_confidence_interval(alpha = 0.05, B = 200, show_progress = FALSE)
#>  2.5% 97.5% 
#>    NA    NA

Everything at once: `InferenceSuite`

InferenceSuite discovers every inference class compatible with this design and response type, runs each applicable procedure, and reports a single Cauchy-combined p-value across them. With screen = TRUE it prints the full results table — the one-call answer to “what does every valid analysis say?” — and returns the same results as an object (res) for programmatic use.

By default it runs every method each class supports, including the resampling ones (randomization, bootstrap, Bayesian bootstrap, jackknife) at their full default replicate counts — a few minutes per class, which is exactly what you want for a real analysis but not inside a vignette that rebuilds on every R CMD check. So this chunk restricts methods to the asymptotic procedures ("wald", "score", "lik_ratio"; classes that don’t support one simply skip it) and sets a per-class time guard. Drop the methods argument to get the complete report.

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/11  [                0%                 ] Status: Estimating...Avg Δ                    mean Δ      3.31      0.820     2.34e-04   wald            ok     
#> Classes 1/11  [===            9%                ] Estimated Time Left: 1sAvg Δ Pooled …           mean Δ      3.31      0.895     4.76e-04   wald            ok     
#> Classes 2/11  [======         18%               ] Estimated Time Left: 1sWilcox                   HL shift    3.24      0.930     9.17e-04   wald            ok     
#> Classes 3/11  [=========      27%               ] Estimated Time Left: 1sLin             ~.       mean Δ      3.35      0.846     2.32e-04   wald            ok     
#> Classes 4/11  [============   36%               ] Estimated Time Left: 1sLin             ~.       mean Δ      3.35      0.846     2.35e-04   score           ok     
#> Classes 5/11  [============== 45%               ] Estimated Time Left: 1sLin             ~.       mean Δ      3.35      0.846     2.35e-04   lik_ratio       ok     
#> Classes 6/11  [============== 54%               ] Estimated Time Left: 1sOLS             ~.       mean Δ      3.30      0.888     4.78e-04   wald            ok     
#> Classes 7/11  [============== 63% ==            ] Estimated Time Left: 0sOLS             ~.       mean Δ      3.30      0.888     2.04e-04   score           ok     
#> Classes 8/11  [============== 72% =====         ] Estimated Time Left: 0sOLS             ~.       mean Δ      3.30      0.888     2.04e-04   lik_ratio       ok     
#> Classes 9/11  [============== 81% ========      ] Estimated Time Left: 0sMedian Regr     ~.       median ef…  2.95      1.20      1.68e-02   wald            ok     
#> Classes 10/11 [============== 90% ===========   ] Estimated Time Left: 0sRobust Regr     ~.       mean Δ      3.29      0.919     7.19e-04   wald            ok     
#> Classes 11/11 [============= 100% ==============] Estimated Time Left: 0s-------------------------------------------------------------------------------------------
#> Status: Completed in 1s.
#> 
#>   Estimand: HL shift (1 inferences)     : p =       NA
#>   Estimand: mean Δ (9 inferences)       : p = 0.000277
#>   Estimand: median effect (1 inferences): p =       NA
#> 
#> Combined evidence against the sharp null across 3 estimands
#> (11 inferences, weighting = uniform within estimand):
#> p = 0.00063

A sequential design: subjects arrive one at a time

Most real trials enroll sequentially. DesignSeqOneByOneKK21 — the Kapelner–Krieger matching-on-the-fly design — decides each arrival’s treatment by trying to match it to an earlier unmatched subject on the covariates (using the reservoir of unmatched subjects otherwise), so covariate balance is built as the trial runs. The API changes only at step 2: assign and record per subject.

des_seq = DesignSeqOneByOneKK21$new(n = n, response_type = "continuous", verbose = FALSE)
for (i in seq_len(n)) {
  w_i = des_seq$add_one_subject_to_experiment_and_assign(X[i, , drop = FALSE])
  y_i = 10 + true_effect * w_i + 0.1 * X$age[i] - 0.3 * X$bmi[i] + rnorm(1, sd = 3)
  des_seq$add_one_subject_response(i, y_i)
}

The matched inference class for this design, InferenceContinKKOLSOneLik, fits one stacked OLS regression over the matched-pair differences and the reservoir subjects, so both contribute to a single likelihood. Note the nonparametric bootstrap is not offered on sequential designs whose assignment depends on earlier subjects (row resampling would not replicate the design); the randomization test, which replays the design’s actual mechanism, is the right tool here.

inf_seq = InferenceContinKKOLSOneLik$new(des_seq, verbose = FALSE)
inf_seq$num_cores = 1L
inf_seq$compute_estimate()
#> [1] 2.738078
inf_seq$compute_asymp_confidence_interval(alpha = 0.05)
#>     2.5%    97.5% 
#> 1.494033 3.982123
inf_seq$set_seed(1)
inf_seq$compute_rand_two_sided_pval(r = 200, show_progress = FALSE)
#> [1] 0.01

Where to go next