Changelog
Source:NEWS.md
sampleSizeCalc 3.1.1 (Development)
Corrections (P0 series: wrong formulas / wrong citations)
-
steppedWedgehadsigma^2andtau^2swapped in Hussey & Hughes Eq 8. The printed equation (p. 187) isVar = I*sigma^2*(sigma^2 + T*tau^2) / [(IU-W)*sigma^2 + (U^2+ITU-TW-IV)*tau^2], wheresigma^2is the cluster-mean variance (sigma2_ind/m) andtau^2the between-cluster variance. The previous code had the two factors interchanged. For the paper’s own Washington EPT example (I=24, T=5, m=100, p=0.05, CV=0.3) it returned SE=0.00486 and power=0.959 at RR=0.64, whereas Eq 8 gives SE=0.006638 and power=0.774, consistent with the paper’s “about 80% power to detect a decrease in prevalence of roughly 36%”. The spurious “80% crossing at RR=0.73” claim is withdrawn. -
twoOrdinalwas missing the allocation factor(1+phi)^2/phi. Machin et al. (2018) Eq 4.3 (p. 57) isN = Gamma (1+phi)^2/phi (z_{1-a/2}+z_{1-b})^2/(log OR)^2withGamma = 3/(1-sum pbar^3)(Eq 4.4); Whitehead (1993) Eq 10 (p. 2260) isn = 3(A+1)^2 (u_{a/2}+u_b)^2 / [A theta^2 (1-sum pbar^3)]. They are algebraically identical withA = phi, and atphi = 1the factor is 4. The previous code omitted it and returned one quarter of the required size (e.g. 16 instead of 61 for Machin Example 4.1). Newallocationargument (default 1) exposes unequal allocation. Verified against Machin Example 4.1 (N4 = 60.09 -> 60) and Whitehead Example 1 (n = 187, 94 per group). -
kappa_Q()now implements Cantor (1996) Eq 3 verbatim. The previousp_o(1-p_o)/(1-p_e)^2coincides with Cantor’s Q only when both marginal rating proportions are 0.5. Cantor’s expression is `Q = (1-pi_e)^(-4) * { sum_i pi_ii[(1-pi_e) - (pi_i. + pi_.i)(1-pi_o)]^2- (1-pi_o)^2 * sum_{i!=j} pi_ij(pi_.i + pi_j.)^2 - (pi_o pi_e - 2pi_e + pi_o)^2 }
with the Eq 4 identities. The implementation now reproduces **Cantor Table 1** at all 18 spot-checked cells and **Table 2** at all 5 maxima ((0.3,0.3) 1.07003 at kappa=0.187,(0.4,0.4) 1.00558 at 0.067,(0.5,0.5) 1.00000 at 0,(0.7,0.3) 0.52438 at 0,(0.9,0.1) 0.04819 at 0`), plus both worked examples (165 and 214). The earlier claim that “pi = 0.3 gives Q = 1.0” was arithmetically impossible.
- (1-pi_o)^2 * sum_{i!=j} pi_ij(pi_.i + pi_j.)^2 - (pi_o pi_e - 2pi_e + pi_o)^2 }
-
AUC variance default is now binormal (Zhou Eq 6.5). The exponential form (Eq 6.4) systematically under-estimates variance for ordinal results (Zhou p. 197), and the binormal form is what Zhou’s own CTC example (p. 198, n = 46) and Hajian-Tilaki’s Table 7 use.
var_methodremains selectable. This makesdiagAUCreproduce the book’sn = 13CTC superiority anchor (the previous one-sidedZ_alphaconvention gave 10) anddiagAUCComparereproduce Hajian’sn = 211example (was 185). -
singleMeannow iterates the Student t quantile (Machin 2018 Eq 9.18, p. 152), instead of using the normal z once. Machin explicitly requires the iteration when N0 < 40. sigma=3, d=0.5 now gives 141 (was 139; Bhardwaj’s z-only arithmetic prints 138.29 -> 138). Thehintsfield reports the starting value when it is below 40. -
cohort+rateswaps the two rates whenrate1 < rate0and says so inhints, so the exposed/unexposed person-time columns stay meaningful.ratiois now documented consistently as the exposed:unexposed person-time ratio (Machin Eq 6.2, p. 84). -
singleArmPhaseIIminimax now matches Simon/clinfun exactly. The tie-break is “minimise N, then minimise EN(p0)”, reproducing Machin 2018 Example 17.2(1/15, 5/25)and clinfun::ph2simon for both (0.1, 0.3) and (0.2, 0.4). -
clusterRCTupgraded to Machin (2018) Eq 12.10 (Student t correction). The previous implementation wrapped the individually-randomised answer inDE = 1+(m-1)ICCand rounded twice (once on subjects, once on clusters), which over-stated the cluster count. The new implementation follows Eq 12.10 (p. 199): `N = DE * [ (1+phi)^2/phi * (z_{1-a/2}+z_{1-b})2/(delta/sigma)2- z^2_{1-a/2}/2 ]
, rounds **once at the cluster level**, and reproduces **Machin Example 12.1 (p. 206) exactly**: delta=3.5, sigma=6.94, m=2, ICC=0.025, 2-sided alpha=0.05, power=0.90 -> **N=172, K=86, 43 per arm** (was 176/44). The ICC sensitivity in the same example is also reproduced (rho=0.1 -> 184, rho=0.2 -> 204). Newcluster_size_inflation()exposes Eq 12.5 (constant m) and Eq 12.7 (variable m viacv_m), and Eq 12.12/12.13 give the binarysigma_Planandrho_Binary`.
- z^2_{1-a/2}/2 ]
-
crossoverreported a double-counted total. Julious (2004) Eq 6.6 (Section 6.3 “Cross-over trials”) writesnas the total sample size, son_requiredis now the per-sequence count andtotal_n = 2 * n_required. The previous code used the total asn_requiredand then doubled it again. Example: delta=2, sigma_w=3 now gives 18 per sequence / 36 total (was 36/72). Senn (2002) Section 9.5.1 (Eq 9.18-9.20, pp. 278-280) confirms the parameterisation (sigma^2 = 2 sigma_w^2,Var = sigma^2/n,df = n-2,n= total withn/2per sequence) and refines it with a Student t correction; thehintsfield now says so. -
mrmc_power()now implements Zhou Eq 6.32 verbatim. The previous code used(1 - rho1)instead of(1 - rho_b), divided the patient-sampling variancesigma_c2byn, and added a spuriousrho2*sigma_b2*((1-rho1)+(J-1)(rho2-rho3))term. With the book’s own mammography parameters it returned 0.979 (the book reports 89%); the old “0.888” anchor came from a different parameter set. The function now has arho_bargument (Table 6.3p_b), which defaults torho1for backward compatibility, and reproduces the book exactly:lambda = 14.0672,power = 0.893. -
Legacy AUC and likelihood-ratio paths no longer return
Inf.calculate_auc_single(),calculate_auc_superiority()andcalculate_auc_comparison()used a mis-expanded Zhou Eq 6.4 variance (A(1-A) + (1/K)(Q1-A^2) + (Q2-A^2)) and divided byK = (1-P)/P, which is zero at the historical defaultprevalence = 1. They now use the printed Eq 6.4 variance functionQ1/R + Q2 - A^2(1/R + 1); whenprevalenceis omitted the non-diseased:diseased ratioRfalls back to 1 (a 1:1 case:control quota) and the deviation is reported inhints.calculate_likelihood_ratio_test()had the sameP/(1-P)division and now falls back identically. This changes previously published numbers:auc_single(auc=0.8, d=0.05, prevalence=0.2)is 120 (was 366),aucSuperiority26 (was 66),aucComparison65 (was 122),likelihoodRatioTest35 (was 120). -
A’Hern (2001) citation corrected in
calculate_ahern(),calculate_diagnostic_superiority_single()and the Chinese guide: the paper is Statistics in Medicine 20(6), 859-866 (doi:10.1002/sim.721), not British Journal of Cancer 85(2), 223-226. - Two more citation corrections verified against the source PDFs: Obuchowski (1994) is Investigative Radiology 29(2), 238-243 (not Medical Decision Making 14(2), 154-162); Simel et al. (1991) is Journal of Clinical Epidemiology 44(8), 763-770 (not 44(11), 1123-1130). Julious (2004) bioequivalence/crossover is Section 6.3 with Eq 6.1-6.9 (the non-central t power expression is Eq 6.8), not “Eq 13, p. 1941”.
-
alternativeis now a real control argument ofcalculate(). It reproduces either Akoglu (2022) Table 2 Eq 2 (one-sided, the default that preserves existing numbers) or Hajian-Tilaki (2014) Eq 6.8 (two-sided). It is currently consumed bydiagAccuracy+test+proportiononly. The unimplementedsamplingparameter has been removed frominst/reference/DESIGN_MATRIX.md, and the partial scope ofdropoutis now stated explicitly there and in?calculate.
Output & Reporting Enhancements
-
Comprehensive Output Display: Default
print()andsummary()methods now consistently output:- Plain-text mathematical formula (
Formula (plain text)) - LaTeX mathematical formula (
Formula (LaTeX)) - Step-by-step calculation process (
Calculation process), documenting hypotheses, input parameters, critical values, intermediate components, and rounding - Literature references (
Reference)
- Plain-text mathematical formula (
- Calculation reports now respect the selected AUC variance model and allocation, retain non-default power for diagnostic CI-width planning, and show the rounding used before cohort dropout inflation.
-
Calculation Steps In S3 Object: All
sample_size_resultobjects now carrycalculationandcalculation_stepsfields documenting the intermediate steps.
Numerical Enhancements & Cross-Package Alignments
diagLR: Corrected the log-LR variance to use sensitivity and false-positive probability for LR+, and false-negative probability and specificity for LR-. The published Se=0.90, Sp=0.50, upper LR-=0.40 example now requires 80 cases and 80 controls, rather than 9 per group. Both formula representations and the clinical guides use the corrected calculation.singleCorrelation: Fixed implementation to strictly match Machin et al. (2018, Chapter 19, Equations 19.5-19.8). Now supports bothprecision(absolute CI half-width) andepsilon(relative precision). Produces exact agreement withpresize::prec_cor(e.g., r=0.5, conf.width=0.10 -> 867; r=0.5, conf.width=0.20 -> 219).bioequivalence: Upgraded 2x2 crossover TOST from crude normal approximation to exact noncentral t distribution (Phillips 1990; Julious 2004 Equation 13). Added support fortheta0(expected geometric mean ratio, default 1.0) and custom equivalence bounds (margin_low,margin_high). Results matchPowerTOST::sampleN.TOSTbit-for-bit (e.g., CV=0.25, theta0=1.0 -> N=24, 12/sequence; theta0=0.95 -> N=28, 14/sequence).diagAccuracy(test): Addedmethod = c("asymptotic", "exact"). Settingmethod = "exact"executes an exact binomial search per A’Hern (2001), returning critical successescritical_successesalongside sample size (e.g., prob=0.85, p0=0.70 gives n=49, critical cutoff r=40).-
Methodological hints and parameter role clarifications added:
-
steppedWedge: Explicit warning regarding the absence of calendar-time fixed effects in Hussey & Hughes (2007) Equation 8 and potential power loss. -
singleArmPhaseII: Documented that ties in minimax total N are broken by minimizing Type I error alpha. -
caseControlMatched: Clarified thatpsiis the discordant pair proportion (pi10 + pi01), not the discordant odds ratio. -
cohortRate: Clarified that required subjects = person-time / follow-up duration. -
diagLR: Clarified that planning targets the approximate CI limit at the anticipated sensitivity and specificity, distinct from fixed-width CI estimation or assurance about every realized interval.
-
sampleSizeCalc 3.1.0
3.1 — new designs (mrmcPower, steppedWedge, kappa, Riley prediction models)
-
predictiveModel+ task=“plan” + outcome=“binary”: Riley 2020 BMJ method for prediction model development sample size (original paper obtained and implemented). Criteria B1 (outcome proportion precision), B3 (shrinkage target via Cox-Snell R2), and B2 (van Smeden MAPE, P <= 30) with the entropy-based max R2cs helper. Anchors exact: B3 P=30/R2cs=0.05/S=0.9 -> 5249 (matches pmsampsize), B1 phi=0.5/delta=0.05 -> 385, B2 MAPE=0.05/ phi=0.3/P=10 -> 461 (Fig 2 worked example). -
predictiveModel+ task=“plan” + outcome=“rate”: Riley 2020 criterion T1 for time-to-event outcomes — person-years of follow-up needed for precise cumulative incidence under an exponential model (exact binomial-CI search; anchor: rate=0.10/yr, t=10y, delta=0.05 -> 2366 person-years, exactly the paper’s Figure 1 worked example). -
mrmc_power(): power calculation for MRMC studies (paired-patient, paired-reader) using the Obuchowski-Rockette mixed model / non-central F framework (Zhou et al. 2011 Eq 6.32, Tables 6.3-6.4). Verified against the book’s mammography worked example (power ~= 0.89 at n=50, J=8, delta=0.10). -
steppedWedge+ task=“test”: Hussey-Hughes Eq 8 WLS Wald power for full stepped wedge designs (original paper obtained and Eq 8 implemented directly; the EPT figure-read anchor discrepancy is documented — Eq 8 gives power 0.959 at RR=0.64 and the 80% crossing at RR=0.73). -
kappa+ task=“estimate”/“test”: Cohen’s kappa sample size per Cantor (1996), using the Fleiss-Cohen-Everitt variance factor Q = po(1-po)/(1-pe)^2. Anchors exact: Q(0.5,0.5,0.7) = 0.510 and the two-questionnaire test N = 214. (Table 1/2 QMAX values require the full FCE lambda1 term; estimate mode is conditional on the planning kappa.)
Breaking / API
- New unified three-axis entry point
calculate(design, task, outcome, ...). Combinations are validated against a master table; illegal combinations error with the list of legal ones. Seeinst/reference/DESIGN_MATRIX.md. (sample_size_calculator()remains fully and silently compatible: all 17 legacy keywords now dispatch through the same shared code path, and golden-master tests assert identical results between both APIs.) -
calculate()resolves arguments by exact name only. This fixes a latent partial-matching bug where design-specific arguments such asdcould be silently swallowed into thedesignformal.
New features (Phase 5)
-
clusterRCT+ task=“test” + outcome=“mean”/“proportion”: design-effect wrapper (DE = 1+(m-1)ICC) around the individually-randomised kernels, with cluster-level rounding (n_clusters per arm = ceil(nDE/m)) (Machin 2018 Ch.12; Donner & Klar 2000). -
singleArmPhaseII: exact Simon two-stage designs (optimal = min E[N|p0], minimax = min N) via vectorised exact binomial search (Simon 1989), plus A’Hern single-stage (A’Hern 2001). Anchors verified against Machin 2018 Chapter 17: optimal (1/10, 5/29) and A’Hern 6/25, 12/35 reproduced exactly; the book’s minimax (2/15, 5/25) has exact power 0.7696 < 0.80 under the strict definition and is excluded (documented discrepancy). -
rct+ task=“noninferiority” + outcome=“proportion”: fills the legacy margin gap (nonInferiorityProportions). -
twoOrdinal: ordered categorical two-group comparison via Whitehead’s formula N = 3(Za/2+Zb)2/((logOR)2(1-sum(p^3))) (Whitehead 1993; Machin 2018 Ch.4). -
predictiveModel+ task=“plan”: events-per-variable rule (EPV >= 10, Peduzzi 1996) with a pointer to Riley et al. (2019). -
its+ task=“test”: interrupted time series level-change power with an AR(1) effective-sample-size correction - marked EXPERIMENTAL in output hints (framework: Bernal et al. 2017; transparent approximation, not a validated source formula). - Deferred with documented reasons:
kappaandsteppedWedge(Machin Ch.20/Ch.13 formulas are OCR-garbled in the available PDF; require paper derivation before implementation).
New features (Phase 4)
- Diagnostic non-inferiority/equivalence (margin designs):
-
diagAccuracyCompare+ task=“noninferiority”/“equivalence”: two- proportion margin framework with one-sided alpha for NI and Z_{beta/2} for equivalence, plus prevalence screening conversion. -
diagAUCCompare+ task=“noninferiority”/“equivalence”: margin framework on the AUC difference with correlation r and all three variance functions (framework per Chow et al. 2017; application per Zhou 2011 Section 6.3).
-
-
diagAccuracyCompare+ paired=TRUE via the three-axis API now uses the Connor kernel with psi = discordant-pair PROPORTION (the legacy keyword keeps the old sqrt-scale parameterization, locked by test). -
crossover: 2x2 crossover with within-subject SD parameterization (n = 2(Za/2+Zb)^2 sigma_w^2 / delta^2; Senn 2002). -
bioequivalence: 2x2 TOST on the log scale, 80-125% criterion (n = 2(Z_{1-a}+Z_{1-b})^2 log(1+CV2)/log(1.25)2; Schuirmann 1987; Julious 2004), with a replicate-design hint for CV > 30%. -
diagAccuracy+ dual=TRUE: dual-endpoint strategy computing sensitivity and specificity requirements each at dual_power (default 0.90 -> ~80% joint) and taking the max (Akoglu 2022 p.184). - Dispatcher gains a
pathflag (“legacy” vs “three-axis”) so upgrades can diverge cleanly from locked legacy behavior.
New features (Phase 3)
- Observational study designs:
-
cohort+ task=“test” + outcome=“proportion”: exposed vs unexposed cumulative incidence with 1:m allocation (ratio), optionaldropoutinflation, and RR parameter roles (Bhardwaj 2024 p.2559 anchor; exact-z value 294, paper’s rounded-z 293, ±1 documented). -
cohort+ outcome=“rate”: person-time incidence-rate comparison with unequal allocation (Machin 2018 Ch.6; IRR roles). -
caseControlMatched: individually matched 1:1 case-control via the discordant-pairs formulation with the independence-within-pair approximation (Dupont 1988; Schlesselman 1982 p.174). -
crosstab: cross-sectional prevalence comparison (PR/OR roles), sharing the two-proportion kernel.
-
-
paired+ outcome=“proportion”: McNemar/Connor paired binary kernel with explicit discordant-pair proportionpsi(Connor 1987; Akoglu 2022 Table 6A uncorrected anchor 155 at psi = delta = 0.05, 80% power). -
rctSingleArm+ task=“test”: single-arm vs fixed value for mean (Machin Ch.5) and proportion outcomes (Hajian-Tilaki Eq 6.8 anchor 153). -
inflate_n()gainsmethod = "additive"(n*(1+d), matching the arithmetic of the Bhardwaj worked examples) alongside the canonical “inflate” (n/(1-d)). - Result objects promote new metric fields:
n_exposed,n_unexposed,person_time_exposed,person_time_unexposed,expected_events,n_pairs.
New features
- Three estimation designs (task = “estimate”):
-
singleProportion: precision-based n = Z2*p(1-p)/d2 with absolute or relative precision and extreme-prevalence convention hints (Bhardwaj 2024 p.2558; anchor 323). -
singleMean: n = (Z*sigma/d)^2 (Machin 2018 Ch.9). -
singleCorrelation: Fisher-z precision n = ((1+r2)/(r*eps))2 + 3 (Machin 2018 Ch.19).
-
- Diagnostic estimation upgrades (three-axis calls only; legacy keywords unchanged):
-
diagAccuracy+ci_method: “precision” (Zhou Eq 6.2) or “ci_width_power” (Zhou Eq 6.1, Flahault 2005; anchors 246/503). -
diagAUC+var_method: “exponential” (Zhou Eq 6.4), “binormal” (Obuchowski 1994, Zhou Eq 6.5; CTC anchor 46/92), “conservative” (distribution-free, Zhou Eq 6.6; anchor 62/123); plus explicit case: controlratioparameter. -
diagLR+lr_type/target_lr: log-LR CI method (Hajian-Tilaki 2014 Sections 6.4, Eqs 6.13/6.14) replacing the simplified Se/Sp-max formula for new-style calls.
-
- New exported utility
screen_required(n, prevalence, beta): total screening enrollment with an exact binomial probability guarantee of capturing n cases (Zhou 2011 Eq 6.3; anchor 451 at 350 cases, 80% prevalence, 90% guarantee). - Result objects gain
hints(methodological warnings surfaced to the user) andn_required(single-sample estimation results).
New features (infrastructure)
- New
sample_size_resultinternals:new_sample_size_result(),validate_sample_size_result(),assert_sample_size_result()(internal three-piece constructor family) with canonical result fields (n_1,events,cases,subjects_to_screen, …) derived from legacy storage fields, preparing for new design modules in 3.0. - Internal numeric kernels
z_alpha(),z_power(),inflate_n(),design_effect()shared by all future design modules (dropout inflation n/(1-d); cluster design effect 1+(m-1)*ICC). -
as_sample_size_result()assembler so each design module can stay a plain calculation list.
Documentation
-
inst/reference/DESIGN_MATRIX.md: full design x task x outcome combination table with formulas, sources, and the 17 legacy-keyword mappings. -
inst/reference/NUMERIC_ANCHORS.md: published worked examples used as numeric test anchors, including documented discrepancies (Bhardwaj 2024 dropout formula vs its own worked example).
Tests
- Golden master suite: all 17 legacy keywords asserted identical between
sample_size_calculator()andcalculate()(228 assertions in the calculate suite alone).