--- title: "testflow: complete pedagogical reference" output: rmarkdown::pdf_document header-includes: - \usepackage{fvextra} - \DefineVerbatimEnvironment{Highlighting}{Verbatim}{breaklines,breakanywhere,commandchars=\\\{\}} vignette: > %\VignetteIndexEntry{testflow: complete pedagogical reference} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, echo = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>", eval = TRUE, fig.width = 6, fig.height = 4) options(width = 80) ``` ```{r setup} library(testflow) cardio <- make_cardio_data() ``` # Introduction This vignette is the complete pedagogical reference for `testflow`: every exported workflow function, organized by study design, alongside the exact mathematical formula each one evaluates. It supersedes and consolidates the previous `effect-size-formulas` and `statistical-test-workflows` vignettes, and folds in the sample-size planning module so that every formula in the package - hypothesis testing, effect sizes, and sample-size planning - lives in one place. Every `test_*()` function returns a `testflow` object: a list containing the cleaned data, descriptive statistics, assumption checks, the recommended test with its rationale, the primary test result (with its null hypothesis), alternative test results for comparison, an effect size, an optional `ggplot`, the original call, and a plain-language interpretation. Five methods work on every object: ```{r} x <- test_two_groups(sbp_3m ~ sex, data = cardio) x # print.testflow(): formatted console report summary(x) # summary.testflow(): compact summary list report(x) # report.testflow(): report-ready sentence as_tibble(x) # as_tibble.testflow(): one-row tidy summary ``` `plot(x)` produces the workflow's `ggplot`; every workflow's own plot is shown alongside its formula in the sections below. `report_test(x)` is a plain-function alias for `report(x)`, provided for callers who prefer not to rely on S3 dispatch. `make_cardio_data(n, seed)` (used throughout this vignette) simulates the small teaching dataset used in every non-survival, non-diagnostic, non-agreement example below - it is a data generator, not a statistical workflow, and carries no formula of its own. Notation used throughout: \(n\) is a sample size, \(\bar x\) and \(s_x\) are a sample mean and standard deviation, \(z_p\) and \(t_{df,p}\) are the \(p\)-th quantiles of the standard normal and \(t\) distributions, and \(\chi^2_{df}\), \(F_{df_1,df_2}\) name a chi-square and \(F\) distribution with the given degrees of freedom. ## How test selection works Most `test_*()` functions run assumption checks first and use them to pick among two or three candidate tests (the *recommended* test), while still computing and reporting every candidate under `alternative_tests` so you can compare. The decision rules (`R/recommend.R`) are simple and deliberately transparent: | Workflow | Normal? & equal variance? | Normal only | Neither | |---|---|---|---| | One sample | One-sample t-test | - | Wilcoxon signed-rank | | Two groups | Student t-test | Welch t-test | Wilcoxon rank-sum | | Paired | Paired t-test | - | Wilcoxon signed-rank | | 3+ groups | One-way ANOVA | Welch ANOVA | Kruskal-Wallis | Normality is checked with Shapiro-Wilk (per group, where applicable); variance homogeneity with Levene's test. `test_categorical()` instead compares chi-square vs. Fisher's exact test based on expected cell counts (`fisher_threshold`, default 5); `test_correlation(method = "auto")` compares Pearson vs. Spearman based on normality and outlier screens; the numeric repeated-measures workflow compares repeated-measures ANOVA vs. Friedman based on within-time normality. # Descriptive statistics `sumtab()` and the internal `descriptives_numeric()`/`descriptives_categorical()` helpers (populating every `testflow` object's `$descriptives` field) report, for numeric variables: \[ n,\quad \bar x = \frac{1}{n}\sum_i x_i,\quad s_x = \sqrt{\frac{1}{n-1}\sum_i (x_i-\bar x)^2},\quad \text{median},\quad IQR = Q_3-Q_1,\quad \min,\ \max \] and, for categorical variables, the count and percentage of each level: ```{r} sumtab(~ age + sbp_3m + ldl | treatment, data = cardio) ``` # Assumption checks `testflow` runs a small set of standardized assumption checks and reports them in every object's `$assumptions` field, each with a `status` (`acceptable`/`not acceptable`/`warning`/`not applicable`/`assumed`/`not checked`) driving the recommended-test logic above. | Check | Method | Statistic | |---|---|---| | Normality | Shapiro-Wilk (`n <= 5000`) | \(W\) | | Symmetry (one-sample/paired) | Cabilio-Masaro | \(z=\sqrt n(\bar x-\text{median})/(s_x\sqrt{\pi/2-1})\) | | Variance homogeneity (2 groups) | `var.test()` (F-test) | \(F=s_1^2/s_2^2\) | | Variance homogeneity (3+ groups) | Levene's test | \(F\) on absolute deviations from group medians | | Variance homogeneity (3+ groups) | Bartlett's test | \(\chi^2\) on group variances | | Sphericity (3+ repeated measures) | Mauchly's test | \(W\) | | Linearity (correlation) | Quadratic-term F-test (RESET, power 2) | \(F\) comparing `y ~ x` to `y ~ x + I(x^2)` | | Monotonicity (correlation) | Spearman correlation of \(x,y\) | \(\rho\) | | Expected cell counts (categorical) | Pearson chi-square expected counts | \(\min(E_{ij}) \geq\) threshold (default 5) | | Multicollinearity (regression) | Variance inflation factor | \(VIF_j = 1/(1-R_j^2)\), flagged at \(\geq 5\) | | Influential points (logistic regression) | Cook's distance | flagged at \(> 4/n\) | | Proportional hazards (Cox) | Schoenfeld residual test (Grambsch-Therneau) | \(\chi^2\) | | Outliers (screen) | IQR rule | \(x < Q_1-1.5\,IQR\) or \(x > Q_3+1.5\,IQR\) | Independence of observations is *assumed* from the study design in every workflow (`testflow` has no way to test it from the data alone) and is always listed explicitly rather than silently taken for granted. # One numerical variable `test_one_sample()` compares a numeric sample against a fixed reference value \(\mu_0\). The recommended test follows normality of the sample: \[ t = \frac{\bar x - \mu_0}{s_x/\sqrt n} \sim t_{n-1} \qquad\text{(One-sample t-test, Gosset 1908)} \] with the Wilcoxon signed-rank test as the non-parametric alternative, and a sign test always reported alongside for reference. The reported effect size is Cohen's \(d\): \[ d = \frac{\bar x - \mu_0}{s_x} \] ```{r} x_one_sample <- test_one_sample(cardio, sbp_3m, mu = 140) x_one_sample plot(x_one_sample) ``` # Two independent groups `test_two_groups()` compares a numeric outcome between two independent groups, selecting among the Student t-test, Welch t-test, and Wilcoxon rank-sum test per the table above: \[ t = \frac{\bar x_1-\bar x_2}{s_p\sqrt{1/n_1+1/n_2}} \sim t_{n_1+n_2-2}, \qquad s_p = \sqrt{\frac{(n_1-1)s_1^2+(n_2-1)s_2^2}{n_1+n_2-2}} \qquad\text{(Student/Gosset 1908)} \] Welch's t-test (1947) instead uses \(s_1^2/n_1+s_2^2/n_2\) as the variance estimate with a Satterthwaite-approximated denominator degrees of freedom - appropriate when variances are unequal. The Wilcoxon rank-sum (Mann-Whitney) test compares rank sums rather than means. Effect size is Cohen's \(d\) (t-test branches) or the rank-biserial correlation \(r_{rb}=1-2W/(n_1n_2)\) (Wilcoxon branch), both with the sign following "first group minus second group." ```{r} x_two_groups <- test_two_groups(sbp_3m ~ sex, data = cardio) x_two_groups plot(x_two_groups) ``` # Paired measurements `test_paired()` compares two measurements on the same subjects via their paired differences \(d_i = \text{after}_i - \text{before}_i\): \[ t = \frac{\bar d}{s_d/\sqrt n} \sim t_{n-1} \qquad\text{(Paired t-test)} \] with the Wilcoxon signed-rank test on the differences as the non-parametric alternative (selected when the differences fail the normality check), and a sign test reported for reference. Effect size is Cohen's \(d_z=\bar d/s_d\). ```{r} x_paired <- test_paired(sbp_3m ~ sbp_baseline, data = cardio) x_paired plot(x_paired) ``` # More than two groups `test_groups()` extends the two-group comparison to \(k\geq 3\) groups: \[ F = \frac{SS_{between}/(k-1)}{SS_{within}/(n-k)} \sim F_{k-1,\,n-k} \qquad\text{(One-way ANOVA, Fisher 1925)} \] Welch's ANOVA relaxes the equal-variance assumption (analogous to Welch's t-test); the Kruskal-Wallis test is the rank-based, distribution-free alternative: \[ H = \frac{12}{n(n+1)}\sum_{j=1}^k \frac{R_j^2}{n_j} - 3(n+1) \sim \chi^2_{k-1} \qquad\text{(Kruskal-Wallis, 1952)} \] Effect size is eta squared \(\eta^2 = SS_{between}/SS_{total}\) for the ANOVA branches, or Kruskal's epsilon squared \(\epsilon^2 = (H-k+1)/(n-k)\) for the Kruskal-Wallis branch. ```{r} x_groups <- test_groups(sbp_3m ~ treatment, data = cardio) x_groups plot(x_groups) ``` ## Post-hoc comparisons When `posthoc = TRUE` (the default), pairwise comparisons follow the selected omnibus test: Tukey's HSD (1949) after one-way ANOVA, BH-adjusted pairwise Welch t-tests after Welch ANOVA, and BH-adjusted pairwise Wilcoxon tests (Dunn-style multiple comparisons, 1964) after Kruskal-Wallis. # Factorial designs `test_factorial()` fits a two- (or more-) way ANOVA with `outcome ~ f1 * f2` and reports sums of squares of the chosen `type`: - **Type I** (`stats::anova()`): sequential SS, order-dependent for unbalanced designs. - **Type II**/**Type III** (`car::Anova()`): each term's SS adjusted for every other term (Type II) or additionally for the intercept under sum-to-zero contrasts (Type III, `testflow` sets `contr.sum` automatically since Type III under the default treatment contrasts tests a different, usually meaningless, hypothesis). Effect size is eta squared per term, computed from whichever SS table (Type I/II/III) was actually used for the reported test - not always recomputed as Type I - so the effect size matches the p-value. Residual normality (Shapiro-Wilk) and variance homogeneity (Levene, on the first factor) are checked as assumptions. ```{r} x_factorial <- test_factorial(sbp_3m ~ sex * treatment, data = cardio) x_factorial plot(x_factorial) ``` # Repeated measurements (numeric) `test_repeated()` (wide data) and `test_repeated_long()` (long data) compare 3 or more repeated numeric measurements on the same subjects, choosing between repeated-measures ANOVA and the Friedman test based on within-time normality: \[ F = \frac{MS_{time}}{MS_{time \times subject}} \sim F_{k-1,\,(k-1)(n-1)} \qquad\text{(Repeated-measures ANOVA, via `Error(id/time)`)} \] \[ \chi_F^2 = \frac{12}{nk(k+1)}\sum_{j=1}^k R_j^2 - 3n(k+1) \sim \chi^2_{k-1} \qquad\text{(Friedman, 1937)} \] Sphericity (Mauchly's test) is checked for the ANOVA branch. Effect size is \(\eta^2_{time}=SS_{time}/SS_{total}\) (ANOVA branch) or Kendall's \(W=\chi_F^2/(n(k-1))\) (Friedman branch). Post-hoc comparisons are BH-adjusted pairwise paired t-tests (ANOVA branch) or pairwise paired Wilcoxon tests (Friedman branch). ```{r} x_repeated <- test_repeated(cardio, c(sbp_baseline, sbp_3m, sbp_6m), id = id) x_repeated plot(x_repeated) ``` `test_repeated_long()` runs the identical workflow (same formulas, same recommendation logic, same plot) starting from long-format data instead of wide: pass the outcome, the within-subject time/condition column, and the subject id column directly, rather than a set of wide columns to pivot. ```{r} long_sbp <- tidyr::pivot_longer( cardio, c(sbp_baseline, sbp_3m, sbp_6m), names_to = "time", values_to = "sbp" ) test_repeated_long(long_sbp, outcome = sbp, within = time, id = id) ``` # Categorical outcomes ## Association between two categorical variables `test_categorical()` tests independence in a contingency table, selecting chi-square or Fisher's exact test based on expected cell counts: \[ \chi^2 = \sum_{i,j}\frac{(O_{ij}-E_{ij})^2}{E_{ij}} \sim \chi^2_{(r-1)(c-1)} \qquad\text{(Pearson 1900)} \] Fisher's exact test (1922) is used instead whenever any expected count falls below `fisher_threshold` (default 5). Effect size is Cramer's V: \[ V = \sqrt{\frac{\chi^2}{n(\min(r,c)-1)}} \] ```{r} x_categorical <- test_categorical(treatment ~ controlled_3m, data = cardio) x_categorical plot(x_categorical) ``` ## One-sample proportion `test_proportion()` tests a single proportion against a reference value \(p_0\), using the exact binomial test (Clopper-Pearson-based) when expected successes or failures fall below 5, or the standard proportion test otherwise: \[ z = \frac{\hat p - p_0}{\sqrt{p_0(1-p_0)/n}} \qquad\text{(one-sample proportion test, normal approximation)} \] ```{r} x_proportion <- test_proportion(cardio, controlled_3m, success = "yes", p = 0.5) x_proportion plot(x_proportion) ``` ## Multinomial goodness of fit `test_multinomial()` compares observed category counts to expected probabilities \(p_1,\ldots,p_k\) (equal by default): \[ \chi^2 = \sum_{i=1}^k \frac{(O_i-E_i)^2}{E_i} \sim \chi^2_{k-1} \] BH-adjusted per-category binomial tests are reported alongside as post-hoc comparisons. ```{r} x_multinomial <- test_multinomial(cardio, treatment) x_multinomial plot(x_multinomial) ``` ## Paired categorical measurements `test_paired_categorical()` compares a binary categorical outcome measured twice on the same subjects via McNemar's test on the discordant cells \(b\) (before-only) and \(c\) (after-only) of the 2x2 table: \[ \chi^2 = \frac{(b-c)^2}{b+c} \sim \chi^2_1 \qquad\text{(McNemar, 1947)} \] Effect size is the number of discordant pairs \(b+c\). ```{r} x_paired_cat <- test_paired_categorical(cardio, controlled_baseline, controlled_3m) x_paired_cat plot(x_paired_cat) ``` ## Repeated categorical measurements (3+ time points) `test_repeated_categorical()` extends McNemar's test to \(k\geq 3\) repeated binary measurements via Cochran's Q test: \[ Q = \frac{(k-1)\left[k\sum_j C_j^2 - \left(\sum_j C_j\right)^2\right]} {k\sum_i R_i - \sum_i R_i^2} \sim \chi^2_{k-1} \qquad\text{(Cochran, 1950)} \] where \(C_j\) is the column (time) total and \(R_i\) the row (subject) total. Effect size is the Cochran-Q analogue of Kendall's \(W\), \(W=Q/(n(k-1))\); pairwise McNemar tests (BH-adjusted) are the post-hoc comparisons. ```{r} x_repeated_cat <- test_repeated_categorical(cardio, c(controlled_baseline, controlled_3m, controlled_6m)) x_repeated_cat plot(x_repeated_cat) ``` # Correlation `test_correlation()` tests association between two numeric variables, with `method = "auto"` choosing Pearson when both variables pass the normality and outlier screens, and Spearman otherwise: \[ r = \frac{\sum_i(x_i-\bar x)(y_i-\bar y)}{\sqrt{\sum_i(x_i-\bar x)^2}\sqrt{\sum_i(y_i-\bar y)^2}} \qquad\text{(Pearson, 1895)} \] \[ \rho = 1-\frac{6\sum_i d_i^2}{n(n^2-1)}, \quad d_i = \text{rank}(x_i)-\text{rank}(y_i) \qquad\text{(Spearman, 1904)} \] Kendall's \(\tau\) (1938), based on the number of concordant vs. discordant pairs, is always computed and reported in `alternative_tests`. Each of the three is a candidate effect size, reported as the estimate of the selected method. Linearity (for Pearson) and monotonicity (for Spearman/Kendall) are checked as assumptions. ```{r} x_correlation <- test_correlation(sbp_3m ~ age, data = cardio) x_correlation plot(x_correlation) ``` `test_correlation_matrix()` extends this to all pairwise correlations among 3+ numeric variables, using pairwise-complete observations per pair and Benjamini-Hochberg-adjusted p-values across all pairs. ```{r} x_corr_matrix <- test_correlation_matrix(cardio, c(age, sbp_baseline, ldl, crp), method = "spearman") x_corr_matrix plot(x_corr_matrix) ``` # Regression ## Linear regression `test_linear_regression()` fits ordinary least squares: \[ y_i = \beta_0+\beta_1x_{1i}+\cdots+\beta_px_{pi}+\varepsilon_i \] and reports the overall model F-test as the primary result: \[ F=\frac{(SS_{total}-SS_{res})/p}{SS_{res}/(n-p-1)} \sim F_{p,\,n-p-1} \] \[ R^2=1-\frac{SS_{res}}{SS_{total}}, \qquad R^2_{adj}=1-(1-R^2)\frac{n-1}{n-p-1} \] Per-term coefficients with confidence intervals are in `alternative_tests$coefficients`. Residual normality (Shapiro-Wilk), homoscedasticity (Breusch-Pagan score test), and multicollinearity (VIF, when \(p>1\)) are checked as assumptions. ```{r} x_linreg <- test_linear_regression(sbp_3m ~ age + ldl, data = cardio) x_linreg plot(x_linreg) ``` ## Logistic regression `test_logistic_regression()` fits: \[ \log\frac{P(y_i=1)}{1-P(y_i=1)}=\beta_0+\beta_1x_{1i}+\cdots+\beta_px_{pi} \] with the overall likelihood-ratio test against the intercept-only model as the primary result: \[ LR = D_{null}-D_{model} \sim \chi^2_p, \qquad D=-2\log L \] \[ R^2_{McFadden}=1-\frac{\log L_{model}}{\log L_{null}} \] Coefficients are reported on both the log-odds and odds-ratio (\(OR_j=e^{\beta_j}\)) scales. Multicollinearity (VIF) and influential observations (Cook's distance, flagged at \(>4/n\)) are checked as assumptions. ```{r} x_logreg <- test_logistic_regression(controlled_3m ~ age + ldl, data = cardio) x_logreg plot(x_logreg) ``` # Survival analysis `make_cardio_data()` has no time-to-event columns, so this section uses a small synthetic example. `Surv()` is re-exported from the `survival` package. ```{r} set.seed(1) n <- 100 survival_dat <- tibble::tibble( time = rexp(n, 0.1), status = rbinom(n, 1, 0.7), arm = rep(c("control", "treatment"), each = n / 2), age = rnorm(n, 60, 10) ) ``` ## Log-rank test `test_survival()` compares two Kaplan-Meier survival curves via the log-rank test: \[ \chi^2 = \sum_{g\in\{A,B\}}\frac{(O_g-E_g)^2}{V_g} \sim \chi^2_1 \qquad\text{(Mantel, 1966)} \] with a companion univariate Cox hazard ratio \(HR=e^{\hat\beta}\) reported as the effect size; the hazard ratio is not itself part of the log-rank test and can disagree with it under strongly non-proportional hazards. ```{r} x_survival <- test_survival(Surv(time, status) ~ arm, data = survival_dat) x_survival plot(x_survival) ``` ## Cox proportional hazards regression `test_cox()` leaves the baseline hazard \(h_0(t)\) unspecified: \[ h(t\mid x) = h_0(t)\exp(\beta_1x_1+\cdots+\beta_px_p) \] Coefficients are reported as hazard ratios \(HR_j=e^{\beta_j}\); the primary result is the overall likelihood-ratio test on \(p\) degrees of freedom. Effect size is Harrell's concordance index. The proportional-hazards assumption is checked via the correlation between scaled Schoenfeld residuals and time (Grambsch & Therneau, 1994). ```{r} x_cox <- test_cox(Surv(time, status) ~ age + arm, data = survival_dat) x_cox plot(x_cox) ``` # Diagnostic test accuracy and agreement `make_cardio_data()` has no diagnostic-test or multi-rater columns, so this section uses small synthetic examples. ## Diagnostic accuracy With true/false positives/negatives \(TP,FP,FN,TN\) from a 2x2 test-vs. -reference table: \[ Sensitivity=\frac{TP}{TP+FN}, \quad Specificity=\frac{TN}{TN+FP}, \quad PPV=\frac{TP}{TP+FP}, \quad NPV=\frac{TN}{TN+FN} \] \[ LR^+=\frac{Sensitivity}{1-Specificity}, \qquad LR^-=\frac{1-Sensitivity}{Specificity} \] Sensitivity, specificity, PPV, NPV, and accuracy each get an exact (Clopper-Pearson) confidence interval; likelihood ratios use the closed-form log-scale interval of Simel, Samsa & Matchar (1991), \(SE[\log LR^+]=\sqrt{1/TP-1/(TP+FN)+1/FP-1/(FP+TN)}\). `test_diagnostic()` tests overall accuracy against the no-information rate (the larger reference-class proportion) as the primary result, matching `caret::confusionMatrix()`'s convention. ```{r} set.seed(1) diag_dat <- tibble::tibble( test = c(rep("positive", 55), rep("negative", 98)), reference = c(rep("positive", 45), rep("negative", 10), rep("positive", 8), rep("negative", 90)) ) x_diagnostic <- test_diagnostic(diag_dat, test, reference) x_diagnostic plot(x_diagnostic) ``` ## ROC / AUC `test_roc()` computes the AUC from the Mann-Whitney U relationship (the probability a random positive-class observation outranks a random negative-class one): \[ AUC = \frac{U}{n_1n_0} \] with the closed-form Hanley & McNeil (1982) standard error: \[ SE(AUC)=\sqrt{\frac{AUC(1-AUC)+(n_1-1)(Q_1-AUC^2)+(n_0-1)(Q_2-AUC^2)}{n_1n_0}}, \quad Q_1=\frac{AUC}{2-AUC}, \ Q_2=\frac{2AUC^2}{1+AUC} \] The optimal threshold maximizes Youden's J (1950), \(J=Sensitivity+Specificity-1\). ```{r} roc_dat <- tibble::tibble( marker = c(rnorm(60, 2, 1), rnorm(50, 0, 1)), disease = c(rep("yes", 60), rep("no", 50)) ) x_roc <- test_roc(roc_dat, marker, disease) x_roc plot(x_roc) ``` ## Inter-rater agreement (Cohen's kappa) `test_agreement()` reports Cohen's kappa from observed agreement \(p_o\) and chance-expected agreement \(p_e\): \[ \kappa = \frac{p_o-p_e}{1-p_e} \] Two *different* standard errors are used for two different purposes: the confidence interval uses the large-sample SE of Fleiss, Cohen & Everitt (1969), evaluated at \(\hat\kappa\); the z-test of \(H_0:\kappa=0\) uses the smaller-magnitude null-hypothesis SE of Fleiss (1981), evaluated at \(\kappa=0\). Reusing the CI SE for the null test (an error `testflow` previously made and has since fixed) understates the null SE whenever agreement is materially above chance, inflating the z-statistic and anti-conservatively shrinking the p-value. ```{r} agree_dat <- tibble::tibble( rater1 = sample(c("mild", "moderate", "severe"), 100, replace = TRUE), rater2 = sample(c("mild", "moderate", "severe"), 100, replace = TRUE) ) x_agreement <- test_agreement(agree_dat, rater1, rater2) x_agreement plot(x_agreement) ``` ## Intraclass correlation coefficient `test_icc()` reports three ICC forms from one-way and two-way ANOVA variance decompositions (Shrout & Fleiss, 1979) of \(n\) subjects by \(k\) raters: \[ ICC(1,1)=\frac{BMS-WMS}{BMS+(k-1)WMS} \] \[ ICC(2,1)=\frac{BMS-EMS}{BMS+(k-1)EMS+\frac{k}{n}(JMS-EMS)} \qquad\text{(two-way random, absolute agreement)} \] \[ ICC(3,1)=\frac{BMS-EMS}{BMS+(k-1)EMS} \qquad\text{(two-way fixed, consistency)} \] ICC(2,1) is reported as the primary/effect-size estimate, following Koo & Li's (2016) recommendation for typical reliability studies. Confidence intervals use the F-distribution formulas of McGraw & Wong (1996), with a Satterthwaite-approximated denominator degrees of freedom for ICC(2,1). ```{r} icc_dat <- tibble::tibble( rater1 = rnorm(30, 50, 10), rater2 = rnorm(30, 50, 10), rater3 = rnorm(30, 50, 10) ) x_icc <- test_icc(icc_dat, c(rater1, rater2, rater3)) x_icc plot(x_icc) ``` # Outlier screening `test_outliers()` is a screening workflow, not a hypothesis test: it flags rows via the IQR rule (univariate, per variable), \[ x < Q_1-1.5\,IQR \quad\text{or}\quad x > Q_3+1.5\,IQR \] and/or Mahalanobis distance (multivariate, `method = "mahalanobis"`), flagged against a \(\chi^2\) cutoff at the 97.5th percentile with degrees of freedom equal to the number of variables: \[ D_i^2 = (x_i-\bar x)^\top S^{-1}(x_i-\bar x), \qquad \text{flag if } D_i^2 > \chi^2_{0.975,\,p} \] ```{r} x_outliers <- test_outliers(c(sbp_3m, ldl, crp), data = cardio) x_outliers plot(x_outliers) ``` # Effect-size reference The table below consolidates every effect-size formula reported across the workflows above, since some functions (Kendall's W, discordant pairs, observed proportion) are only summarized inline elsewhere. | Workflow | Effect size | Formula | |---|---|---| | One sample | Cohen's d | \(d=(\bar x-\mu_0)/s_x\) | | Two groups (t-test) | Cohen's d | \(d=(\bar x_1-\bar x_2)/s_p\) | | Two groups (Wilcoxon) | Rank-biserial r | \(r_{rb}=1-2W/(n_1n_2)\) | | Paired | Cohen's dz | \(d_z=\bar d/s_d\) | | ANOVA (one-way, factorial) | Eta squared | \(\eta_j^2=SS_j/\sum SS\) | | Kruskal-Wallis | Epsilon squared | \(\epsilon^2=(H-k+1)/(n-k)\) | | Repeated ANOVA | Eta squared (time) | \(\eta^2_{time}=SS_{time}/SS_{total}\) | | Friedman | Kendall's W | \(W=\chi_F^2/(n(k-1))\) | | Cochran Q | Kendall's-W analogue | \(W=Q/(n(k-1))\) | | Categorical association | Cramer's V | \(V=\sqrt{\chi^2/(n(\min(r,c)-1))}\) | | One-sample proportion | Observed proportion | \(\hat p=x/n\) | | Multinomial | Chi-square GOF statistic | \(\chi^2=\sum_i(O_i-E_i)^2/E_i\) | | Paired categorical | Discordant pairs | \(b+c\) | | Correlation | r / rho / tau | see "Correlation" above | | Linear regression | R squared, adjusted R squared | \(R^2=1-SS_{res}/SS_{total}\) | | Logistic regression | McFadden's pseudo R squared | \(1-\log L_{model}/\log L_{null}\) | | Survival (log-rank) | Hazard ratio | \(HR=e^{\hat\beta}\) (no magnitude label) | | Cox regression | Concordance index | Harrell's C (no magnitude label) | | Diagnostic accuracy | Accuracy | \((TP+TN)/n\) (no magnitude label) | | ROC | AUC | \(U/(n_1n_0)\) | | Agreement | Cohen's kappa | \((p_o-p_e)/(1-p_e)\) | | ICC | ICC(2,1) | see "Intraclass correlation" above | ## Magnitude labels - Cohen's \(d\): negligible \(<0.2\), small \(<0.5\), moderate \(<0.8\), otherwise large. - Eta/epsilon/R-squared family: negligible \(<0.01\), small \(<0.06\), moderate \(<0.14\), otherwise large. (McFadden's pseudo \(R^2\) is not on the OLS \(R^2\) scale; values above ~0.2-0.4 are already a good fit in practice.) - Cramer's V, rank-biserial \(|r|\), Kendall's W: negligible \(<0.1\), small \(<0.3\), moderate \(<0.5\), otherwise large. - AUC: poor if \(\max(AUC,1-AUC)<0.7\), acceptable \(<0.8\), excellent \(<0.9\), otherwise outstanding (Hosmer, Lemeshow & Sturdivant 2013, p. 177). - Cohen's kappa: poor \(<0\), slight \(<0.21\), fair \(<0.41\), moderate \(<0.61\), substantial \(<0.81\), otherwise almost perfect (Landis & Koch, 1977). - ICC: poor \(<0.5\), moderate \(<0.75\), good \(<0.9\), otherwise excellent (Koo & Li, 2016). - Hazard ratio, concordance index, accuracy: no magnitude label - there is no widely agreed small/moderate/large convention, and inventing one would overstate how standardized that judgment is. The estimate itself is always reported. # Sample-size planning `sample_size()` and the endpoint-specific `sample_size_*()` functions plan continuous, binary, survival, ordinal, bioequivalence, and precision-based designs, each under superiority, non-inferiority, or equivalence objectives where applicable. Every result is a `sample_size` object with a raw and dropout-adjusted \(n\), the assumptions used, a report-ready sentence, and (for some designs) a power curve. Notation: \(\alpha\) is the type I error rate, \(1-\beta\) is power, \(z_{1-\alpha}\)/\(z_{1-\alpha/2}\)/\(z_{1-\beta}\) are standard normal quantiles, \(r=n_B/n_A\) is the allocation ratio. Superiority uses two-sided \(\alpha\); non-inferiority/equivalence use one-sided \(\alpha\) (pass `alpha = 0.025` explicitly - `testflow` does not rescale it). Dropout adjustment is always: \[ n_{recruit} = \left\lceil \frac{n_{eval}}{1-d} \right\rceil \] ```{r, eval = FALSE} sample_size( endpoint = c("continuous", "binary", "survival", "ordinal"), design = c("parallel", "paired", "repeated"), objective = c("superiority", "noninferiority", "equivalence"), ... ) ``` ## Continuous endpoints `sample_size_continuous()` covers parallel, paired, and repeated designs. For a parallel-group superiority comparison with common variance \(\sigma^2\) and target difference \(\Delta\): \[ n_A = \frac{(r+1)\sigma^2(z_{1-\beta}+z_{1-\alpha/2})^2}{r\Delta^2}, \qquad n_B = rn_A \] ```{r} ss_parallel <- sample_size_continuous(design = "parallel", objective = "superiority", delta = 5, sd = 10, alpha = 0.05, power = 0.90) ss_parallel plot(ss_parallel, type = "summary") ``` Paired designs replace \(\sigma\) with the paired-difference SD and drop the allocation ratio: \[ n_{pairs} = \frac{(z_{1-\alpha/2}+z_{1-\beta})^2\sigma_{diff}^2}{\Delta^2} \] `curve_data` (and therefore `plot(type = "curve"/"both")`) is populated only for `design = "paired"`/`"repeated"`, `objective = "superiority"`; other combinations fall back to the summary bar chart shown above. ```{r} ss_paired <- sample_size_continuous(design = "paired", objective = "superiority", delta = 5, sd_diff = 10, alpha = 0.05, power = 0.90) ss_paired plot(ss_paired, type = "both") ``` For `design = "repeated"` with 3+ time points, \(\sigma_{diff}\) is rescaled under a compound-symmetry working correlation before the paired formula is applied (`testflow` design-effect extension, not a standalone textbook formula): \[ \sigma_{eff} = \sigma_{diff}\sqrt{\frac{1+(m-1)\rho}{m}} \] Non-inferiority uses \(d_{NI}=\Delta+\delta\); equivalence uses \(d_{EQ}=\delta-|\Delta|\) (or, when \(\Delta=0\) exactly, \(z_{1-\beta/2}\) in place of \(z_{1-\beta}\)), both plugged into the same form as superiority with a one-sided \(\alpha\). ## Binary endpoints `sample_size_binary()` covers parallel risk-difference and paired (discordant-pairs) designs. Parallel superiority (`method = "pooled"`, with \(\bar p=(p_A+rp_B)/(1+r)\)): \[ n_A = \frac{\left[z_{1-\alpha/2}\sqrt{(1+1/r)\bar p(1-\bar p)} +z_{1-\beta}\sqrt{p_A(1-p_A)+p_B(1-p_B)/r}\right]^2}{(p_A-p_B)^2} \] ```{r} sample_size_binary(design = "parallel", objective = "superiority", p1 = 0.4, p2 = 0.25, method = "pooled", alpha = 0.05, power = 0.90) ``` Paired binary designs plan on the discordant cells directly, with discordant odds ratio \(OR_D=p_{10}/p_{01}\) and discordance rate \(\lambda_D=p_{10}+p_{01}\): \[ n_{total} = \left\lceil\frac{(z_{1-\alpha/2}+z_{1-\beta})^2(OR_D+1)^2}{(OR_D-1)^2\lambda_D}\right\rceil \] Non-inferiority/equivalence use the same asymmetric variance form at the non-inferiority/equivalence distance instead of \(p_A-p_B\). ## Survival endpoints `sample_size_survival()` plans required events, then converts to total N. Under an assumed exponential hazard: \[ D_{total} = \frac{4(z_{1-\alpha/2}+z_{1-\beta})^2}{[\log(HR)]^2} \] or, proportional-hazards-only (`method = "ph_only"`): \[ D_{total} = \frac{2(HR+1)^2(z_{1-\alpha/2}+z_{1-\beta})^2}{(HR-1)^2} \] Given planned survival probabilities \(S_A(t)\), \(S_B(t)\), equal allocation: \[ N_{total} = \frac{2D_{total}}{(1-S_A(t))+(1-S_B(t))} \] ```{r} sample_size_survival(hr = 0.7, survival_a = 0.8, survival_b = 0.7, alpha = 0.05, power = 0.90) ``` With staggered uniform accrual (`accrual_duration`, `follow_up`), the flat \(1-S_g(t)\) is replaced by the accrual-averaged event probability: \[ \bar q_g = 1-\frac{e^{-\lambda_g(T-R)}-e^{-\lambda_gT}}{\lambda_gR}, \qquad T=R+\text{follow\_up} \] which always requires more subjects than instantaneous accrual for the same event target, and reduces to the flat conversion as \(R\to0\). Non-inferiority/equivalence work on the log-hazard-ratio scale with \(d_{NI}=\log(HR_M)-\log(HR)\) or \(d_{EQ}=\min(\theta-\theta_L,\theta_U-\theta)\). ## Ordinal endpoints `sample_size_ordinal()` implements Noether's method for a two-group Mann-Whitney/Wilcoxon comparison, with \(P=P(A>B)+\tfrac12P(A=B)\): \[ n_{per\ group} = \frac{(z_{1-\alpha/2}+z_{1-\beta})^2}{6(P-0.5)^2} \] ```{r} sample_size_ordinal(p_superiority = 0.65, alpha = 0.05, power = 0.90) ``` ## Bioequivalence `sample_size_bioequivalence()` plans a two one-sided test (TOST) study on the log scale, with \(\theta_0=\log(GMR)\), bounds \(\theta_L=\log(L)\)/\(\theta_U=\log(U)\) (usually 0.80/1.25), and \(d_{BE}=\min(\theta_0-\theta_L,\theta_U-\theta_0)\). The default `method = "iterative_tost"` searches for the smallest \(n\) achieving the *exact* TOST power via the noncentral t distribution (Phillips 1990): \[ Power(n) = T_{df,\delta_U}(-t_{1-\alpha,df}) - T_{df,\delta_L}(t_{1-\alpha,df}) \] where \(T_{df,\delta}\) is the noncentral t CDF, \(\delta_U=(\theta_0-\theta_U)/SE\), \(\delta_L=(\theta_0-\theta_L)/SE\), \(\sigma=\sqrt{\log(1+CV^2)}\). Verified in package tests to match `PowerTOST::power.TOST()` (the regulatory-standard, Owen's-Q-based calculation) for balanced designs; `method = "normal_approx"` gives the closed-form approximation \(n_{total}\approx2\sigma_w^2(z_{1-\beta}+z_{1-\alpha})^2/d_{BE}^2\), which agrees closely near \(GMR=1\) but can diverge further away from it (an earlier version of `testflow` used only the normal approximation here, under-sizing studies by up to ~20% at small \(n\)). ```{r} sample_size_bioequivalence(design = "crossover", gmr = 0.95, cv_within = 0.30, alpha = 0.05, power = 0.90) ``` ## Precision-based sample size `sample_size_precision()` sizes a study to a target confidence-interval half-width \(w\) rather than power against an effect - there is no `power` argument. ```{r} sample_size_precision(endpoint = "continuous", design = "one_sample", width = 2, sd = 10, alpha = 0.05) sample_size_precision(endpoint = "binary", design = "two_sample", width = 0.08, p1 = 0.4, p2 = 0.3, allocation = 2) ``` For `endpoint = "binary"`, `design = "one_sample"`, `method` chooses the interval used for planning: `"wald"` (default, closed-form normal approximation) or `"wilson"`/`"exact"` (Clopper-Pearson), which run a deterministic search for the minimum integer \(n\) whose interval has maximum one-sided half-width \(\leq w\) - important for rare-prevalence studies (newborn screening, rare adverse events) where the Wald approximation is unreliable: ```{r} sample_size_precision(endpoint = "binary", design = "one_sample", width = 0.02, p = 0.01, method = "wilson") ``` `criterion = "expected"` (default) checks precision at the anticipated event count `round(n*p)`; `criterion = "worst_case"` additionally requires the target to hold across nearby plausible event counts. `x$diagnostics` reports rare-event planning quantities a CI-width target alone does not guarantee (`expected_events`, `probability_zero_events`, `complete_case_n` vs. `adjusted_n`), with an `approximation_warning` when the Wald approximation is used or expected counts are small. For binary `design = "two_sample"`, with \(r=n_B/n_A\): \[ n_A = \frac{z_{1-\alpha/2}^2}{w^2}\left[p_1(1-p_1)+\frac{p_2(1-p_2)}{r}\right], \qquad n_B = rn_A \] ## Cluster-randomized adjustment and dropout `sample_size_cluster_adjust()` inflates an individually randomized \(n\) by the design effect \(DE=1+(m-1)\rho\) (equal cluster size \(m\), intracluster correlation \(\rho\)), or \(DE\approx1+((1+CV_m^2)m-1)\rho\) for unevenly sized clusters. `sample_size_adjust_dropout()` exposes the dropout formula above directly. Both return a plain integer, not a `sample_size` object - they are post-processing adjustments, not full planning calculations. ```{r} sample_size_cluster_adjust(100, m = 20, rho = 0.02) sample_size_adjust_dropout(100, dropout = 0.15) ``` ## Summary of supported settings | Endpoint | Design(s) | Objective(s) | Formula family | |---|---|---|---| | Continuous | parallel, paired, repeated (2+ time points) | superiority, non-inferiority, equivalence | normal-approximation mean-difference formulas, compound-symmetry repeated extension | | Binary | parallel, paired, repeated (2 time points) | superiority, non-inferiority, equivalence | risk-difference and discordant-pair planning | | Survival | parallel | superiority, non-inferiority, equivalence | event-based proportional-hazards planning, optional uniform-accrual adjustment | | Ordinal | parallel | superiority | Noether approximation | | Bioequivalence | crossover, parallel | equivalence (TOST) | iterative TOST or normal approximation, log scale | | Precision | one-sample, two-sample, log-OR | CI half-width (no power/objective) | continuous and binary CI-width formulas | Every `sample_size` object also has `print()`, `summary()`, `report()`, `as_tibble()`, and `plot(type = "summary"/"curve"/"both")` methods. See `?sample_size_continuous`, `?sample_size_binary`, `?sample_size_survival`, `?sample_size_ordinal`, `?sample_size_bioequivalence`, `?sample_size_precision`, and `?sample_size_cluster_adjust` for the complete argument reference and every `design`/`objective`/`method` combination. # References Gosset WS (Student). The probable error of a mean. *Biometrika*. 1908;6(1):1-25. Welch BL. The generalization of "Student's" problem when several different population variances are involved. *Biometrika*. 1947;34(1-2):28-35. Wilcoxon F. Individual comparisons by ranking methods. *Biometrics Bulletin*. 1945;1(6):80-83. Mann HB, Whitney DR. On a test of whether one of two random variables is stochastically larger than the other. *Annals of Mathematical Statistics*. 1947;18(1):50-60. Fisher RA. *Statistical Methods for Research Workers*. Oliver and Boyd; 1925. Levene H. Robust tests for equality of variances. In: *Contributions to Probability and Statistics*. Stanford University Press; 1960. Kruskal WH, Wallis WA. Use of ranks in one-criterion variance analysis. *Journal of the American Statistical Association*. 1952;47(260):583-621. Tukey JW. Comparing individual means in the analysis of variance. *Biometrics*. 1949;5(2):99-114. Dunn OJ. Multiple comparisons using rank sums. *Technometrics*. 1964;6(3):241-252. Friedman M. The use of ranks to avoid the assumption of normality implicit in the analysis of variance. *Journal of the American Statistical Association*. 1937;32(200):675-701. Girden ER. *ANOVA: Repeated Measures*. Sage; 1992. Cochran WG. The comparison of percentages in matched samples. *Biometrika*. 1950;37(3/4):256-266. McNemar Q. Note on the sampling error of the difference between correlated proportions or percentages. *Psychometrika*. 1947;12(2):153-157. Pearson K. Notes on regression and inheritance in the case of two parents. *Proceedings of the Royal Society of London*. 1895;58:240-242. Pearson K. On the criterion that a given system of deviations from the probable in the case of a correlated system of variables is such that it can be reasonably supposed to have arisen from random sampling. *Philosophical Magazine*. 1900;50(302):157-175. Fisher RA. On the interpretation of chi-square from contingency tables, and the calculation of P. *Journal of the Royal Statistical Society*. 1922;85(1):87-94. Spearman C. The proof and measurement of association between two things. *The American Journal of Psychology*. 1904;15(1):72-101. Kendall MG. A new measure of rank correlation. *Biometrika*. 1938;30(1/2):81-93. Cramer H. *Mathematical Methods of Statistics*. Princeton University Press; 1946. Clopper CJ, Pearson ES. The use of confidence or fiducial limits illustrated in the case of the binomial. *Biometrika*. 1934;26(4):404-413. Cohen J. *Statistical Power Analysis for the Behavioral Sciences* (2nd ed.). Lawrence Erlbaum; 1988. Draper NR, Smith H. *Applied Regression Analysis* (3rd ed.). Wiley; 1998. Hosmer DW, Lemeshow S, Sturdivant RX. *Applied Logistic Regression* (3rd ed.). Wiley; 2013. Kaplan EL, Meier P. Nonparametric estimation from incomplete observations. *Journal of the American Statistical Association*. 1958;53(282):457-481. Mantel N. Evaluation of survival data and two new rank order statistics arising in its consideration. *Cancer Chemotherapy Reports*. 1966;50(3):163-170. Cox DR. Regression models and life-tables. *Journal of the Royal Statistical Society: Series B*. 1972;34(2):187-202. Grambsch PM, Therneau TM. Proportional hazards tests and diagnostics based on weighted residuals. *Biometrika*. 1994;81(3):515-526. Simel DL, Samsa GP, Matchar DB. Likelihood ratios with confidence: sample size estimation for diagnostic test studies. *Journal of Clinical Epidemiology*. 1991;44(8):763-770. Altman DG, Bland JM. Diagnostic tests 1: sensitivity and specificity. *BMJ*. 1994;308(6943):1552. Hanley JA, McNeil BJ. The meaning and use of the area under a receiver operating characteristic (ROC) curve. *Radiology*. 1982;143(1):29-36. Youden WJ. Index for rating diagnostic tests. *Cancer*. 1950;3(1):32-35. Cohen J. A coefficient of agreement for nominal scales. *Educational and Psychological Measurement*. 1960;20(1):37-46. Fleiss JL, Cohen J, Everitt BS. Large sample standard errors of kappa and weighted kappa. *Psychological Bulletin*. 1969;72(5):323-327. Fleiss JL. *Statistical Methods for Rates and Proportions* (2nd ed.). Wiley; 1981. Landis JR, Koch GG. The measurement of observer agreement for categorical data. *Biometrics*. 1977;33(1):159-174. Shrout PE, Fleiss JL. Intraclass correlations: uses in assessing rater reliability. *Psychological Bulletin*. 1979;86(2):420-428. McGraw KO, Wong SP. Forming inferences about some intraclass correlation coefficients. *Psychological Methods*. 1996;1(1):30-46. Koo TK, Li MY. A guideline of selecting and reporting intraclass correlation coefficients for reliability research. *Journal of Chiropractic Medicine*. 2016;15(2):155-163. Julious SA. *Sample Sizes for Clinical Trials*. Chapman & Hall/CRC; 2010. Phillips KF. Power of the two one-sided tests procedure in bioequivalence. *Journal of Pharmacokinetics and Biopharmaceutics*. 1990;18(2):137-144. Donner A, Klar N. *Design and Analysis of Cluster Randomization Trials in Health Research*. Arnold; 2000.