Reproducing the Paper’s Numerical Results

Reproducibility map

The paper contains three numerical components. Here RCT denotes randomized controlled trial.

four Monte Carlo experiments with analytic variance benchmarks;
an analysis of postoperative throat-pain scores from a licorice-gargle RCT;
deterministic power and total-sample-size calculations under two planning laws.

In method labels, DIM denotes the unadjusted difference in means, ANCOVA denotes common-slope analysis of covariance, and ANHECOVA-IF denotes interacted analysis of heterogeneous covariance with joint influence-function inference. CF denotes cross-fitting. HC0 and HC2 denote the heteroskedasticity-consistent type 0 coefficient sandwich and type 2 leverage correction, respectively.

RCTCovAdj separates original inputs, executable code, and derived results. The simulation laws and the aggregate simulation and study-design results are installed with the package. Neither the licorice-gargle trial records nor results derived from them are distributed. A user with the required permission and an authorized local copy can verify that file before running the real-data analysis.

Simulation laws

In every case,

\[ X\sim N(0,1),\quad A\sim\operatorname{Bernoulli}(\pi),\quad \varepsilon\sim N(0,1), \]

independently, and

\[ Y=0.5A+\{b_0+(b_1-b_0)A\}X+c(X^2-1)+\varepsilon. \]

The marginal mean difference is \(0.5\). The installed case table records the four parameter combinations and the scientific role of each design.

data("simulation_cases", package = "RCTCovAdj")
simulation_cases
#>   case                     label  pi b0 b1 quad
#> 1    A        Nonlinear additive 0.5  1  1 0.75
#> 2    B   Heterogeneous, balanced 0.5  2 -1 0.00
#> 3    C Heterogeneous, unbalanced 0.8  2 -1 0.00
#> 4    D           Opposing slopes 0.5 -1  1 0.00

The designs isolate nonlinear additive prognosis, heterogeneous slopes under balanced allocation, the same slopes under unequal allocation, and opposing slopes whose optimally weighted prognostic signal cancels.

For these laws, the analytic variance factors are

\[ \begin{aligned} V_U&=\frac{1+b_1^2+2c^2}{\pi} +\frac{1+b_0^2+2c^2}{1-\pi},\\ V^*&=(b_1-b_0)^2+\frac{1}{\pi}+\frac{1}{1-\pi},\\ V_{\mathrm{ANH}}&=V_U- \frac{\{(1-\pi)b_1+\pi b_0\}^2}{\pi(1-\pi)},\\ V_{\mathrm{ANC}}&=V_{\mathrm{ANH}}+ \frac{(2\pi-1)^2(b_1-b_0)^2}{\pi(1-\pi)}. \end{aligned} \]

Here \(V_U\) is the variance factor for DIM, \(V^*\) is the semiparametric efficiency bound, \(V_{\mathrm{ANH}}\) is the factor for interacted linear adjustment, and \(V_{\mathrm{ANC}}\) is the factor for common-slope ANCOVA. The Oracle method uses the true treatment-specific conditional means and therefore has variance factor \(V^*\). The cross-fitted quadratic method has the same first-order factor in these four laws because its feature span contains the conditional means. rct_variance_factors() evaluates the displayed expressions and also reports the coefficient-only HC2 variance limit.

with(
  simulation_cases[simulation_cases$case == "A", ],
  rct_variance_factors(
    allocation = pi,
    beta_control = b0,
    beta_treatment = b1,
    quadratic = quad,
    noise_variance = 1
  )
)
#>    unadjusted     efficient    interacted        ancova     hc2_limit 
#>          12.5           4.0           8.5           8.5           8.5 
#> optimal_slope  pooled_slope 
#>           1.0           1.0

A short executable check

simulate_rct_case() generates one trial from a named law. Setting the seed makes the example reproducible.

trial_a <- simulate_rct_case(case = "A", n = 200L, seed = 123L)
str(trial_a)
#> 'data.frame':    200 obs. of  3 variables:
#>  $ x        : num  -0.969 0.706 1.489 -1.815 0.33 ...
#>  $ treatment: int  0 0 0 0 1 1 0 1 0 1 ...
#>  $ outcome  : num  0.205 0.534 1.668 1.165 1.558 ...
#>  - attr(*, "case_parameters")='data.frame':  1 obs. of  6 variables:
#>   ..$ case : chr "A"
#>   ..$ label: chr "Nonlinear additive"
#>   ..$ pi   : num 0.5
#>   ..$ b0   : num 1
#>   ..$ b1   : num 1
#>   ..$ quad : num 0.75

rct_adjust(
  outcome = trial_a$outcome,
  treatment = trial_a$treatment,
  covariates = trial_a["x"],
  allocation = 0.5,
  methods = c("unadjusted", "ancova", "interacted")
)
#> Covariate-adjusted marginal mean differences
#>                      method estimate std.error conf.low conf.high
#>                  Unadjusted   0.5951    0.2304   0.1435    1.0468
#>         Common-slope ANCOVA   0.5544    0.1941   0.1741    0.9348
#>  Interacted standardization   0.5545    0.1940   0.1743    0.9347
#>                                 inference   n n.control n.treated allocation
#>              empirical influence function 200       106        94        0.5
#>                  HC0 coefficient sandwich 200       106        94        0.5
#>  joint standardization influence function 200       106        94        0.5

The quadratic term is absent from the analysis of covariance (ANCOVA) and interacted fits, so their linear span cannot attain the semiparametric bound in Case A. A cross-fitted quadratic learner includes the missing prognostic direction.

rct_crossfit(
  outcome = trial_a$outcome,
  treatment = trial_a$treatment,
  covariates = trial_a["x"],
  allocation = 0.5,
  learner = "quadratic",
  folds = 2L,
  seed = 456L
)
#> Covariate-adjusted marginal mean differences
#>                             method estimate std.error conf.low conf.high
#>  Cross-fitted quadratic adjustment   0.4511    0.1503   0.1566    0.7457
#>                        inference   n n.control n.treated allocation
#>  cross-fitted influence function 200       106        94        0.5

Full Monte Carlo experiment

The paper used sample sizes 200 and 800, 4,000 independent replications in each case–sample-size cell, and master seed 20260903. Distinct L’Ecuyer combined multiple-recursive random-number streams make the cells reproducible without sharing streams. The full experiment creates 32,000 trial datasets and is deliberately excluded from vignette evaluation.

full_output_dir <- tempfile("rctcovadj-full-")
full_files <- reproduce_paper_simulations(
  output_dir = full_output_dir,
  reps = 4000L,
  sample_sizes = c(200L, 800L),
  seed = 20260903L
)
unlink(full_output_dir, recursive = TRUE)

This example removes its temporary output directory after the run. To retain the files, choose an explicit permanent destination and omit the cleanup call.

For a fast installation check, reduce the number of replications. Such a run checks code paths and file structure; it is too small for scientific Monte Carlo conclusions.

quick <- run_paper_simulations(
  reps = 32L,
  sample_sizes = 200L,
  seed = 20260903L
)

Complete aggregate summaries from the 4,000-replication run are included for direct inspection. Replicate-level output is created when the full experiment is rerun with keep_replicates = TRUE.

data("paper_population_benchmarks", package = "RCTCovAdj")
data("paper_simulation_results", package = "RCTCovAdj")

head(paper_population_benchmarks)
#>   case                     label  pi b0 b1 quad   VU Vstar Vlinear Vancova
#> 1    A        Nonlinear additive 0.5  1  1 0.75 12.5  4.00    8.50     8.5
#> 2    B   Heterogeneous, balanced 0.5  2 -1 0.00 14.0 13.00   13.00    13.0
#> 3    C Heterogeneous, unbalanced 0.8  2 -1 0.00 27.5 15.25   15.25    35.5
#> 4    D           Opposing slopes 0.5 -1  1 0.00  8.0  8.00    8.00     8.0
#>   Vhc2_limit beta_opt beta_pool
#> 1       8.50      1.0       1.0
#> 2       4.00      0.5       0.5
#> 3       6.25      1.4      -0.4
#> 4       4.00      0.0       0.0
head(paper_simulation_results)
#>   case   n  pi       method reps          bias   bias_mcse    emp_sd
#> 1    A 200 0.5          DIM 4000 -0.0007503803 0.003922982 0.2481112
#> 2    A 200 0.5       ANCOVA 4000 -0.0007117228 0.003225973 0.2040285
#> 3    A 200 0.5  ANHECOVA-IF 4000 -0.0006166324 0.003243016 0.2051063
#> 4    A 200 0.5 ANHECOVA-HC2 4000 -0.0006166324 0.003243016 0.2051063
#> 5    A 200 0.5    CF-linear 4000 -0.0006726284 0.003290503 0.2081097
#> 6    A 200 0.5 CF-quadratic 4000 -0.0003160130 0.002254944 0.1426152
#>   emp_sd_mcse   mean_se mean_se_mcse coverage coverage_mcse empirical_nvar
#> 1 0.002726944 0.2487015 0.0003426429  0.95025   0.003437839      12.311831
#> 2 0.002271743 0.2035131 0.0002516665  0.94875   0.003486525       8.325522
#> 3 0.002281916 0.2035495 0.0002518441  0.94875   0.003486525       8.413720
#> 4 0.002281916 0.2073628 0.0002745164  0.95400   0.003312250       8.413720
#> 5 0.002316432 0.2094763 0.0002868159  0.95150   0.003396607       8.661930
#> 6 0.001564719 0.1439481 0.0001210376  0.95300   0.003346304       4.067819
#>    nvar_mcse point_estimator_nvar_limit reported_se_nvar_limit
#> 1 0.27063409                       12.5                   12.5
#> 2 0.18540005                        8.5                    8.5
#> 3 0.18721411                        8.5                    8.5
#> 4 0.18721411                        8.5                    8.5
#> 5 0.19282880                        8.5                    8.5
#> 6 0.08926111                        4.0                    4.0
#>   relative_to_bound training_fallbacks full_data_fallbacks empty_arm_fallbacks
#> 1          3.077958                  0                   0                   0
#> 2          2.081380                  0                   0                   0
#> 3          2.103430                  0                   0                   0
#> 4          2.103430                  0                   0                   0
#> 5          2.165482                  0                   0                   0
#> 6          1.016955                  0                   0                   0
#>   hc2_leverage_fallbacks hc2_leverage_floor_candidates
#> 1                      0                             0
#> 2                      0                             0
#> 3                      0                             0
#> 4                      0                             0
#> 5                      0                             0
#> 6                      0                             0

The simulation figures can be rebuilt from the installed aggregate summaries without rerunning the Monte Carlo experiment. In the efficiency figure, \(S_B^2\) is the empirical variance of the estimates across replications, so the vertical axis is \(nS_B^2/V^*\). Filled circles and vertical bars show the \(n=800\) empirical ratios and their 95% Monte Carlo intervals; open diamonds show the corresponding first-order ratios. The short axis label ANH-IF denotes ANHECOVA-IF, and Oracle denotes adjustment with the true conditional means.

plot_simulation_results(
  results = paper_simulation_results,
  benchmarks = paper_population_benchmarks,
  metric = "efficiency"
)

Four panels show empirical root-n variance relative to the semiparametric efficiency bound for six estimators in Cases A through D. Filled circles have vertical 95 percent Monte Carlo intervals, and open diamonds show first-order theory.

Empirical variance relative to the semiparametric bound at n=800. Bars are 95% Monte Carlo intervals for the empirical ratios; open diamonds are the first-order ratios.

The empirical points fluctuate around the analytic limits. Equality is not expected in a finite Monte Carlo run. Comparisons should use the stored Monte Carlo standard errors, especially when two methods have similar variance.

Postoperative throat-pain analysis

The package distributes no trial record, aggregate table, estimate, or figure derived from the licorice-gargle data. The source portal states that permission from the data contributor or corresponding author is required before the data are used in a new publication. Installing RCTCovAdj grants no right to analyze, publish, or redistribute the records or derived output.

After obtaining the required permission, an authorized source file can be verified and analyzed as follows:

source_data <- read_licorice_data(file.choose())
trial_results <- analyze_licorice(source_data)
trial_results$estimates
plot_licorice_results(trial_results)

read_licorice_data() rejects a file whose fingerprint or structure does not match the documented source. This guard prevents an unrelated object with the same file name from silently entering the analysis. Permission to obtain and use the data and its derived output remains the user’s responsibility.

The optional medicaldata package supplies a source-format copy through a separate distribution. Its availability does not establish permission for a proposed use. A user whose permission covers that copy may construct the reviewed analysis object and run the same calculation:

source_data <- prepare_licorice_data()
trial_results <- analyze_licorice(source_data)

This preparation route does not grant permission to analyze, publish, or redistribute the participant records or derived output. The returned object contains the aggregate estimates, working standard errors, bounded-score completion analysis, diagnostics, and model coefficients generated from the authorized local file. These outputs remain local unless the user’s permission allows their publication or redistribution.

Power and sample-size application

The design application is deterministic. It uses analytic variance factors, a marginal planning effect of 0.5, and a two-sided level-0.05 normal-reference test. No Monte Carlo result enters the calculation.

data("paper_power_design_results", package = "RCTCovAdj")
paper_power_design_results
#>   case         analysis      curve variance_factor relative_variance
#> 1    A Unadjusted (DIM) unadjusted           12.50         1.0000000
#> 2    A      ANHECOVA-IF interacted            8.50         0.6800000
#> 3    A     CF-quadratic  efficient            4.00         0.3200000
#> 4    C Unadjusted (DIM) unadjusted           27.50         1.0000000
#> 5    C      ANHECOVA-IF interacted           15.25         0.5545455
#> 6    C    Pooled ANCOVA     pooled           35.50         1.2909091
#>   required_n_80 achieved_power_80 required_n_90 achieved_power_90 power_n_300
#> 1           393         0.8005559           526         0.9003400   0.6877704
#> 2           267         0.8002038           358         0.9005937   0.8438674
#> 3           126         0.8013024           169         0.9014810   0.9911099
#> 4           864         0.8002839          1156         0.9000452   0.3789962
#> 5           479         0.8001797           641         0.9000210   0.6016954
#> 6          1115         0.8001624          1493         0.9001803   0.3065879
#>   power_n_500 power_n_800
#> 1   0.8853791   0.9793266
#> 2   0.9695941   0.9980784
#> 3   0.9998584   0.9999998
#> 4   0.5683196   0.7693904
#> 5   0.8167451   0.9516901
#> 6   0.4667904   0.6604236

The same table and power curves can be regenerated from the package functions.

recomputed_design <- paper_power_design()
recomputed_design$sample_sizes
#>   case         analysis      curve variance_factor relative_variance
#> 1    A Unadjusted (DIM) unadjusted           12.50         1.0000000
#> 2    A      ANHECOVA-IF interacted            8.50         0.6800000
#> 3    A     CF-quadratic  efficient            4.00         0.3200000
#> 4    C Unadjusted (DIM) unadjusted           27.50         1.0000000
#> 5    C      ANHECOVA-IF interacted           15.25         0.5545455
#> 6    C    Pooled ANCOVA     pooled           35.50         1.2909091
#>   required_n_80 achieved_power_80 required_n_90 achieved_power_90 power_n_300
#> 1           393         0.8005559           526         0.9003400   0.6877704
#> 2           267         0.8002038           358         0.9005937   0.8438674
#> 3           126         0.8013024           169         0.9014810   0.9911099
#> 4           864         0.8002839          1156         0.9000452   0.3789962
#> 5           479         0.8001797           641         0.9000210   0.6016954
#> 6          1115         0.8001624          1493         0.9001803   0.3065879
#>   power_n_500 power_n_800
#> 1   0.8853791   0.9793266
#> 2   0.9695941   0.9980784
#> 3   0.9998584   0.9999998
#> 4   0.5683196   0.7693904
#> 5   0.8167451   0.9516901
#> 6   0.4667904   0.6604236
plot_power_design(recomputed_design)

Two-panel normal-approximation power curves under balanced nonlinear Case A and unequal-allocation Case C, comparing unadjusted, interacted, and efficient quadratic or pooled ANCOVA analyses; symbols mark 80 percent power.

Normal-approximation power under planning Cases A and C. Symbols mark the smallest integer total sample sizes attaining 80% power for each analysis.

Integer sample sizes are rounded upward after the continuous normal approximation is inverted. They are total randomized sample sizes before any inflation for attrition. Reproduction of the arithmetic does not validate the planning law for a new trial; its variance inputs still require scientific justification.

Standalone replication scripts

The package installs three command-line scripts for rebuilding the simulation, power, and authorized-data analyses. Their installed location can be found without relying on a machine-specific path.

replication_dir <- system.file("replication", package = "RCTCovAdj")
list.files(replication_dir, pattern = "[.]R$")
#> [1] "run_licorice_application.R" "run_power_application.R"   
#> [3] "run_simulations.R"

Each script requires an explicit output directory, and the licorice analysis also requires the path to an authorized source file. The README.md in the same directory documents the command-line arguments. Every workflow also writes session information and SHA-256 and MD5 checksums. For example, a short simulation check can be launched from a terminal with

Rscript --vanilla run_simulations.R --quick --output=/path/to/results

after changing to the installed replication directory. The full simulation omits --quick and uses the article’s 4,000 replications per cell.

Recommended audit sequence

A compact audit proceeds in the following order:

inspect simulation_cases and recompute the population variance factors;
compare the installed aggregate summaries in paper_simulation_results with those analytic benchmarks;
rerun a small simulation to check the local software environment;
run the 4,000-replication experiment only when the full archive must be rebuilt;
after obtaining permission, verify and analyze the authorized local trial file and inspect the generated aggregate output; and
rebuild the deterministic power and sample-size results separately from the Monte Carlo study.

This sequence distinguishes exact population calculations, stochastic Monte Carlo evidence, restricted source data, and deterministic study-design calculations. That distinction is necessary for a reproducibility audit because the four objects have different sources of uncertainty.