bio-clinical-biostatistics-effect-measures

## Version Compatibility Reference examples tested with: statsmodels 0.14+, numpy 1.26+, pandas 2.1+, matplotlib 3.8+, marginaleffects (Python) 0.0.13+ / (R) 0.20+. R packages cited: ratesci, exact2x2, marginaleffects, riskCommunicator, RobinCar. Before using code patterns, verify installed versions match. If versions differ: - Python: `pip show <package>` then `help(module.function)` to check signatures - R: `packageVersion('<pkg>')` then `?function_name` If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying. # Treatment Effect Measures for Clinical Trials **"Compute treatment effect sizes"** -> Estimate the population-level treatment contrast (OR, RR, RD, HR, NNT) with a confidence interval calibrated to sample size and a clear declaration of whether the estimand is marginal or conditional under ICH E9(R1). ## Algorithmic Taxonomy | Measure | Scale | Collapsible? | Best CI method | When to use | Fails when | |---------|-------|--------------|----------------|-------------|------------| | OR | Log-odds ratio | NO (non-collapsible) | Profile likelihood; Wald acceptable for n>100 per arm | Case-control (only valid measure); logistic regression default | Outcome prevalence > 10% (OR overstates RR); Hauck-Donner pathology near boundary | | RR | Log-risk ratio | YES | Miettinen-Nurminen score; MOVER-R | Cohort, RCT with common outcomes | Sparse strata; one or both p near 0 (Wald log-RR breaks) | | RD (absolute risk difference) | Linear probability | YES | Newcombe-Wilson hybrid; Miettinen-Nurminen | Clinically interpretable absolute scale; FDA-preferred for binary | Predictions outside [0,1] from linear models | | HR (hazard ratio) | Log-hazard ratio | NO | Wald with profile likelihood for small n | Time-to-event with PH | PH violation (see clinical-biostatistics/survival-analysis) | | NNT/NNH | 1/RD | n/a (derived from RD) | Bender 2001 *CCT* 22:102 (Altman 1998 base) | Communicating absolute benefit to clinicians | RD CI crosses zero (NNT becomes NNTB-infinity-NNTH) | | Difference in RMST | Time scale | YES | Wald with delta method; pseudo-observation regression | Time-to-event with PH violation | Different max follow-up across arms (truncation tau ambiguous) | **Postdoc reading:** Permutt 2020 *Stat Biopharm Res* 12:45 ("Do covariates change the estimand?") established that the conditional OR from a multivariable logistic regression is a *different parameter* than the marginal OR -- and the FDA May 2023 Final Guidance "Adjusting for Covariates in RCTs" requires the **marginal** estimand for primary reporting (see clinical-biostatistics/logistic-regression for g-computation/standardisation machinery). ## Decision Tree by Scenario | Scenario | Recommended estimand + CI | Why | |----------|---------------------------|-----| | RCT, binary outcome, prevalence <10% | OR with Wald CI; report as primary | OR ~= RR at low prevalence; standard regulatory currency | | RCT, binary outcome, prevalence >=10% | RR or RD via modified Poisson with HC1/HC3 sandwich SE; OR as secondary | OR substantially overstates RR; modified Poisson directly estimates RR (Zou 2004 *AJE* 159:702) | | RCT, binary outcome, primary endpoint per FDA 2023 | Marginal RD via g-computation with sandwich SE; conditional OR as supportive | FDA 2023 final guidance: marginal estimand for primary; cite Permutt 2020 | | Case-control study | OR only (RR unidentifiable); profile likelihood CI | OR is the only measure estimable from case-control design | | Noninferiority on absolute scale | Miettinen-Nurminen score CI for RD | Regulatory standard for RD CIs; consistent with Pearson chi-square | | Noninferiority on relative scale | Miettinen-Nurminen score CI for RR; Koopman 1984 acceptable | Wald log-RR has poor coverage near boundary; MN is the regulatory expectation | | Stratified design with site/region strata | MH pooled OR with stratified MN CI (R `ratesci::scoreci(..., stratified=TRUE)`) | Preserves stratification; cite Kahan-Morris 2012 | | Reporting NNT for clinicians | Bender 2001 *CCT* 22:102 method; report as NNTB(lower)..infinity..NNTH(upper) when CI crosses zero | Standard in BMJ/Lancet/Cochrane | | Time-to-event with PH violation | RMST difference; cite Royston-Parmar 2013, Uno 2014 | HR is a misleading single-number summary under non-PH (see survival-analysis) | ## Crude Effect Measures from 2x2 Tables **Goal:** Compute unadjusted OR, RR, RD from a contingency table with calibrated CIs. **Approach:** Use Table2x2 for OR/RR with Wald CIs, but switch to score-based methods (Miettinen-Nurminen, Newcombe-Wilson) for regulatory contexts. ```python from statsmodels.stats.contingency_tables import Table2x2 from statsmodels.stats.proportion import confint_proportions_2indep import numpy as np # Table layout: [[treated_event, treated_no_event], [control_event, control_no_event]] table = np.array([[a, b], [c, d]]) t = Table2x2(table) print('OR:', t.oddsratio, t.oddsratio_confint()) # Wald log-OR CI print('RR:', t.riskratio, t.riskratio_confint()) # Wald log-RR CI # Modern CI for RD (regulatory preferred): ci_rd = confint_proportions_2indep(a, a+b, c, c+d, method='newcomb', alpha=0.05) # 'newcomb' = Newcombe-Wilson hybrid; 'agresti-caffo' is the +1+1 +1+1 adjustment # For Miettinen-Nurminen stratified CIs, use R `ratesci::scoreci` ``` **Table orientation is critical:** Table2x2 interprets the first column as "outcome present." `pd.crosstab` orders columns alphabetically; if outcome is coded 0/1, the table will have 0 first and the OR will be the *reciprocal* of intended. Always reorder: `cross = cross[[1, 0]]` or `cross = cross[['Yes', 'No']]`. This is a silent direction-reversing error. **For the underlying significance test of a 2x2:** Pearson chi-square (no Yates) at n>=40 with adequate expected counts; otherwise Boschloo's exact (uniformly more powerful than Fisher's exact at the same Type-I per Mehta-Senchaudhuri 2003; Lydersen-Fagerland-Laake 2009). See clinical-biostatistics/categorical-tests for the algorithmic taxonomy. ## Modern Confidence Intervals -- the Postdoc Toolkit Wald CIs are the textbook default but have well-documented coverage failures (Brown-Cai-DasGupta 2001 *Stat Sci* 16:101 -- "chaotic" coverage near 0 and 1 even with moderate n; Newcombe 1998a *Stat Med* 17:873 documents 11 alternatives for the difference of two proportions). ### Single proportion | Method | When | Status | |--------|------|--------| | Wald | Never preferred | Defunct for serious work; coverage can drop to 0.0 (Brown-Cai-DasGupta) | | Wilson (score) | General default | Brown-Cai-DasGupta recommended; matches Pearson chi-square test | | Jeffreys (Beta(1/2,1/2)) | Small n | Equal-tailed Jeffreys; Bayesian with reference prior; recommended for small n | | Clopper-Pearson | When exact guarantee required | Over-covers by 1-4 percentage points; sometimes regulatory-required | | Agresti-Coull | Simple Wald replacement | +2/+2 adjustment; teaching default | | Mid-p Clopper-Pearson | When CP over-coverage hurts NI margin | Less conservative than CP; slight under-coverage | ```python from statsmodels.stats.proportion import proportion_confint ci = proportion_confint(45, 60, alpha=0.05, method='wilson') # Other methods: 'jeffreys', 'agresti_coull', 'beta' (Clopper-Pearson), 'normal' (Wald) ``` ### Difference of two proportions (RD) | Method | Citation | Verdict | |--------|----------|---------| | Wald | textbook | Poor; can produce limits outside [-1,1] | | Newcombe-Wilson hybrid / MOVER | Newcombe 1998a | Recommended; balanced of computational simplicity and coverage | | Agresti-Caffo (+1+1, +1+1) | Agresti-Caffo 2000 *Am Stat* 54:280 | Add 1 success + 1 failure per arm, then Wald; surprisingly good for small n | | Miettinen-Nurminen score | Miettinen-Nurminen 1985 *Stat Med* 4:213 | **Regulatory standard for RD**; consistent with Pearson chi-square; SAS PROC FREQ `riskdiff(method=mn)` | | Chan-Zhang exact unconditional | Chan-Zhang 1999 *Biometrics* 55:1202 | Guaranteed coverage; needed when MN under-covers near extreme p | ### Ratio of two proportions (RR) | Method | Citation | Verdict | |--------|----------|---------| | Wald on log(RR) ("Katz log") | textbook (Katz 1978) | Defunct for serious work; biased when either p small | | Koopman 1984 score | *Biometrics* 40:513 | Iterative; original score CI for RR | | Miettinen-Nurminen score | Miettinen-Nurminen 1985 | Regulatory standard; R `ratesci::scoreci(contrast='RR')` | | MOVER-R | Donner-Zou 2012 *Stat Methods Med Res* 21:347 | Construct CI on theta1 - R*theta2 not on R directly; allows asymmetric CIs | ### Odds ratio | Method | Citation | Verdict | |--------|----------|---------| | Wald on log(OR) | textbook | Fast; suffers Hauck-Donner effect when cell counts small | | Cornfield exact | Cornfield 1956 | Exact reference standard; often over-conservative | | Mid-p Cornfield | Berry-Armitage 1995 *Statistician* 44:417 | Less conservative; slight under-coverage | | Profile likelihood | Venzon-Moolgavkar 1988 *Appl Stat* 37:87 | **Transformation-invariant; no Hauck-Donner pathology**; R `MASS::confint.glm` default for `glm` | **Hauck-Donner effect (1977 *JASA* 72:851; revived by Yee 2022 *JASA* 117:1763):** the Wald test statistic is *non-monotonic* in the parameter estimate near the boundary -- a large OR can produce a tiny Wald chi-square so the test fails to reject when it should. Use LR or profile-likelihood inference. `VGAM::hdeff()` detects this in fitted models. ## Number Needed to Treat (NNT) -- the Bender 2001 *CCT* 22:102 Way **Goal:** Convert RD to clinically intuitive NNT with a CI that handles the singularity at RD = 0. **Approach:** Compute RD CI first, then transform; when RD CI crosses zero, report NNTB(lower) -> infinity -> NNTH(upper). ```python import numpy as np def nnt_with_ci(treated_events, treated_n, control_events, control_n, alpha=0.05): """Bender 2001 *CCT* 22:102 method; transforms RD CI to NNT CI handling singularity.""" p_t = treated_events / treated_n p_c = control_events / control_n rd = p_c - p_t # positive = treatment helps se_rd = np.sqrt(p_t * (1 - p_t) / treated_n + p_c * (1 - p_c) / control_n) z = 1.96 if alpha == 0.05 else None rd_ci = (rd - z * se_rd, rd + z * se_rd) nnt = 1 / rd if rd != 0 else float('inf') if rd_ci[0] > 0 and rd_ci[1] > 0: return f'NNTB {1/rd_ci[1]:.0f} to {1/rd_ci[0]:.0f}' elif rd_ci[0] < 0 and rd_ci[1] < 0: return f'NNTH {1/-rd_ci[1]:.0f} to {1/-rd_ci[0]:.0f}' else: # CI crosses zero -- Bender convention if rd > 0: return f'NNTB {nnt:.0f} (NNTB {1/rd_ci[1]:.0f} to inf to NNTH {1/-rd_ci[0]:.0f})' return f'NNTH {-nnt:.0f} (similar split)' ``` **Bender 2001 *Controlled Clinical Trials* 22:102-110** resolves the discontinuity at RD = 0: NNT has a singularity there, so a CI that crosses zero produces a disjoint NNT CI ("NNTB(some) -> infinity -> NNTH(some)"). The Cochrane/BMJ convention is to report exactly this -- the infinity in the middle signals non-significance and is more honest than truncating to one side. **NNT from OR + baseline risk** (when only the OR is published): ```python def nnt_from_or(odds_ratio, baseline_risk): baseline_odds = baseline_risk / (1 - baseline_risk) treatment_odds = baseline_odds * odds_ratio treatment_risk = treatment_odds / (1 + treatment_odds) arr = abs(baseline_risk - treatment_risk) return 1 / arr if arr > 0 else float('inf') ``` | OR | Baseline 5% | Baseline 20% | Baseline 50% | |----|-------------|--------------|--------------| | 0.5 | NNT=42 | NNT=12 | NNT=6 | | 0.7 | NNT=70 | NNT=20 | NNT=12 | **Always report baseline risk alongside NNT** -- the same OR produces dramatically different NNTs. ## Marginal RD via G-Computation -- The FDA 2023 Recipe **Goal:** Compute the marginal RD (FDA 2023 primary estimand for binary endpoints) from a fitted logistic regression -- without requiring the analyst to switch model class. **Approach:** Standardise over the observed covariate distribution. Predict per-subject probability under treatment=1 AND under treatment=0; take the mean of each; the difference is the marginal RD. SE via influence function or bootstrap. ```python import statsmodels.formula.api as smf import numpy as np # Step 1: Fit logistic with covariates fit = smf.logit('y ~ z + x1 + x2 + x3', data=df).fit() # Step 2: Predict under each treatment regime (counterfactual prediction) df_z1 = df.assign(z=1) df_z0 = df.assign(z=0) p_z1 = fit.predict(df_z1) p_z0 = fit.predict(df_z0) # Step 3: Marginal estimates marg_p1 = p_z1.mean() marg_p0 = p_z0.mean() marg_rd = marg_p1 - marg_p0 # FDA 2023 primary marg_rr = marg_p1 / marg_p0 marg_or = (marg_p1 / (1 - marg_p1)) / (marg_p0 / (1 - marg_p0)) # Step 4: SE via bootstrap (or analytical influence function in R `marginaleffects`) n_boot = 1000 boot_rds = np.zeros(n_boot) rng = np.random.default_rng(42) for b in range(n_boot): idx = rng.choice(len(df), size=len(df), replace=True) fit_b = smf.logit('y ~ z + x1 + x2 + x3', data=df.iloc[idx]).fit(disp=0) p1_b = fit_b.predict(df_z1.iloc[idx]).mean() p0_b = fit_b.predict(df_z0.iloc[idx]).mean() boot_rds[b] = p1_b - p0_b se_marg_rd = boot_rds.std(ddof=1) ci_marg_rd = (marg_rd - 1.96*se_marg_rd, marg_rd + 1.96*se_marg_rd) ``` **R equivalent (preferred for confirmatory):**