How to Choose Statistical Tests for Medical Papers: t-test, ANOVA, Chi-squared, Nonparametric Tests, Correlation and Regression

Decision table

Choosing a test by data type, number of groups, pairing and distribution

Determine the outcome type and the study design first, then look at the distribution. "First choice" is the approach most medical journal reviewers accept; "alternatives" apply when conditions are not met or as sensitivity analyses.

Decide between independent and paired by how the data were generated. The left and right side of the same animal, a patient before and after treatment, and cases and controls matched 1:1 on age and sex are all paired data. Several lesions, teeth or follow-up visits from the same patient are not independent observations either; the unit of analysis must match the unit of randomization or sampling.

When comparing two groups on an ordinal outcome (for example cured, markedly improved, improved, no effect), use the Mann-Whitney U test. A chi-squared test discards the ordering.

Assumptions

Normality and equal variances: how to check them and which practices are mistakes

t-tests and ANOVA require approximate normality within each group (or of regression residuals), not of all data pooled. In large samples they are quite robust to non-normality; the main risk is unequal variances combined with unequal sample sizes.

Tested: two groups

Student's t, Welch's t and Mann-Whitney U on the same data

Tested on this page: R 4.4.3 (car 3.1.5, effectsize 1.0.3, rstatix 1.1.0, afex 1.5.1, survival 3.8.12) and Python 3.12 (SciPy 1.18.1, statsmodels 0.15.0, pingouin 0.7.0, lifelines 0.30.3) on macOS arm64. The data are R built-in datasets exported to CSV so both languages read identical data. The data are ozone concentrations (ppb) from airquality, right-skewed and with clearly different variances between months.

In the June vs August row, the larger group also has the larger variance. Student's t pools the variances toward the larger group, overestimates the standard error and gives a larger P value (0.034 vs 0.0043). If instead the smaller group has the larger variance, Student's t gives too small a P value and more false positives. Welch's t needs no prior judgment in either case.

Effect sizes and confidence intervals: May vs August mean difference −36.3 ppb (Welch 95% CI −54.4 to −18.3), Hedges' g=−1.11 (95% CI −1.69 to −0.53); June vs August mean difference −30.5 ppb (95% CI −50.7 to −10.4), Hedges' g with unpooled SD −0.96 (95% CI −1.61 to −0.30).

Paired data analyzed as independent: the sleep dataset records extra sleep for 10 subjects under each of two drugs. The paired t-test gives t=4.06, df=9, P=0.0028, mean difference 1.58 hours (95% CI 0.70 to 2.46). An independent-samples t-test gives t=1.86, df=18, P=0.079. Between-subject variation is not removed, and a clear difference goes undetected.

Several groups and repeated measures

Three or more groups: ANOVA, post hoc tests and repeated measures

The overall test answers whether all groups are equal; post hoc tests answer which pairs differ. Report both.

The ctrl vs trt2 comparison shows the problem with pairwise t-tests: unadjusted P=0.048 would be written as "statistically significant" at 0.05, while the Tukey-adjusted P=0.198. Three groups need 3 comparisons; five groups need 10, and the inflation grows.

Kruskal-Wallis (airquality ozone, May to September, n=116): H=29.27, df=4, P=6.9×10⁻⁶, η²H=0.23. Post hoc Dunn tests with Holm adjustment show differences for May vs July, May vs August, July vs September and August vs September. Do not follow Kruskal-Wallis with unadjusted pairwise Mann-Whitney tests.

Repeated measures (Indometh, indomethacin plasma concentrations in 6 subjects at 0.5, 1, 2, 4 and 8 hours): Mauchly's W=0.00063, P=0.0059, so sphericity fails; the uncorrected F(4, 20)=106.6; Greenhouse-Geisser ε=0.378, corrected degrees of freedom 1.51 and 7.55, P=4.7×10⁻⁶. R's afex::aov_ez and pingouin.rm_anova agree. Friedman χ²=24.0, P=8.0×10⁻⁵. State which correction was used and the ε value in the paper.

• The Tukey test works with unequal group sizes (Tukey-Kramer method); R's TukeyHSD and SPSS's Tukey handle this automatically.
• When only comparing each treatment with one control, use Dunnett's test: fewer comparisons and more power than Tukey.
• The LSD method controls the familywise error rate only with exactly 3 groups and a significant overall F test; do not use it with more groups.
• The SNK method does not control the familywise error rate when only some group means are equal; reviewers often ask for Tukey or Holm instead.
• With dropouts, repeated-measures ANOVA deletes the whole subject; a linear mixed model uses incomplete observations.

Multiple comparisons

Bonferroni, Holm, Tukey and FDR: which to use when

With 10 independent tests whose null hypotheses are all true, the probability of at least one P<0.05 is 1−0.95¹⁰=40%. The adjustment method depends on the structure of the comparisons and the nature of the study.

Six of these raw P values are below 0.05; only the first remains after Bonferroni or Holm adjustment, and the first three after BH adjustment. In the paper, state the adjustment method and the family of tests adjusted (for example, "Holm adjustment for the 6 between-group comparisons of 3 secondary outcomes"), and report adjusted P values or the adjusted significance level.

raw P    Bonferroni  Holm    BH(FDR)
0.001    0.010       0.010   0.0100
0.008    0.080       0.072   0.0400
0.012    0.120       0.096   0.0400
0.021    0.210       0.147   0.0525
0.035    0.350       0.210   0.0700
0.048    0.480       0.240   0.0800
0.090    0.900       0.360   0.1286
0.200    1.000       0.600   0.2500
0.410    1.000       0.820   0.4556
0.740    1.000       0.820   0.7400

Categorical data

Chi-squared, Fisher's exact and McNemar's tests

For 2×2 tables, check the total sample size and the smallest expected count before choosing a test. R's chisq.test and SciPy's chi2_contingency both apply the Yates continuity correction to 2×2 tables by default; SPSS prints both the uncorrected and corrected rows.

The rule in Chinese medical statistics textbooks for 2×2 tables: with n≥40 and all expected counts T≥5, use Pearson's chi-squared; with n≥40 and some 1≤T<5, use the continuity-corrected chi-squared; with n<40 or any T<1, use Fisher's exact test. R×C tables need no cell with T<1 and no more than 1/5 of cells with 1≤T<5.

Campbell (2007, Stat Med) compared methods by simulation and found the Yates correction too conservative. He recommended the N−1 chi-squared test (Pearson χ² multiplied by (N−1)/N) when all expected counts are at least 1, and the Fisher test otherwise. For Chinese journals, follow the textbook rule; for international journals, Fisher's exact test or the uncorrected chi-squared can be reported directly with the reason stated.

When reporting an OR from the Fisher test, note that R gives the conditional maximum likelihood estimate and SciPy gives the sample OR ad/bc; with small samples they differ noticeably (8.15 vs 9.33). State how the OR and its 95% CI were computed, or estimate it with logistic regression throughout.

McNemar's test (the 1600-person repeated poll from the R help page, discordant pairs b=150, c=86): continuity-corrected χ²=16.82, P=4.1×10⁻⁵; uncorrected χ²=17.36, P=3.1×10⁻⁵; exact binomial test on discordant pairs P=3.7×10⁻⁵. R's mcnemar.test applies the correction by default, while statsmodels' mcnemar defaults to exact=True, so the default results differ.

Correlation and regression

Correlation analysis, regression analysis and common misreadings

A correlation coefficient describes the strength of a linear or monotonic association; a regression coefficient describes an effect size adjusted for other variables. Neither alone supports a causal conclusion.

The number of predictors in logistic and Cox regression is limited by the number of events. Peduzzi et al. (1996) proposed at least 10 events per variable (EPV≥10) from simulation; later work treats this only as a rough lower bound, and sample sizes for prediction models should be calculated with the method of Riley et al. Putting 10 predictors into a model with only 30 events gives severely biased coefficients.

Common errors

The statistical errors reviewers flag most often

Each of these can be checked before submission; each item gives the consequence and the fix.

Reporting

Writing effect sizes, 95% CIs and exact P values

Compiled from the SAMPL guidelines (Lang and Altman). Where the target journal has its own statistical requirements, follow the journal.

Choosing an effect size: use Hedges' g for two-group means (less biased than Cohen's d in small samples); for paired designs, dz (mean difference divided by the SD of differences) is common. Software defines "paired Cohen's d" differently: on the sleep data dz=1.28, while pingouin's ttest(paired=True) reports cohen_d=0.83 (using the average SD of the two measurements as the denominator). State the formula in the paper. Cohen's 0.2, 0.5 and 0.8 benchmarks come from the behavioral sciences; in medical research, interpret differences in original units against the minimal clinically important difference.

• Report P values as equalities: P=0.03, P=0.22, not P<0.05 or "NS". The smallest P value that needs reporting is P<0.001, except in genetic association studies.
• For primary outcomes, report effect estimates with 95% CIs, such as mean differences, risk differences, ORs and HRs.
• State which method was used for each analysis instead of listing all method names at the end of the statistics paragraph.
• State whether tests were one- or two-sided (justify one-sided tests), the alpha level, and whether and how multiple comparisons were adjusted.
• State how assumptions were checked: skewed data were analyzed with nonparametric methods, paired data with paired methods, and linear regression checked for linearity and residuals.
• Give numerators and denominators for percentages and the sample size for each analysis.
• Name the statistical software and version, for example R 4.4.3 (car 3.1) or SPSS 27.0.
• Label post hoc and unplanned subgroup analyses as exploratory.

Software map

R, Python and SPSS equivalents and their default differences

The three packages use different defaults for the same method, so the same data can give different P values. R and Python results in the table were run on this page; SPSS menus were checked against IBM documentation and not run locally.

R aov, cylinders entered first (type I)      P = 0.056
R aov, transmission entered first (type I)   P = 4.9e-07
car::Anova type = 2 / pingouin               P = 0.056
car::Anova type = 3 (contr.sum)              P = 0.083
statsmodels anova_lm(typ=3)                  P = 0.083   # matches SPSS default type III

Runnable code

Complete R and Python code

The code below ran in this page's test environment. The R code uses built-in datasets directly; the Python code reads CSV files exported from R so both sides use identical data.

If rstatix and effectsize are both loaded in R, both provide a cohens_d function; the one loaded later masks the other, and incompatible arguments raise unused argument (pooled_sd = FALSE). Writing effectsize::cohens_d avoids this.

pingouin is slow to import the first time (about 80 seconds on the test machine while it builds its cache) and normal afterwards.

# R 4.4：install.packages(c("car", "effectsize", "rstatix", "afex", "survival"))
library(car); library(effectsize); library(rstatix); library(afex); library(survival)

# Two independent groups: airquality ozone, June (n=9) vs August (n=26)
aq <- subset(airquality, Month %in% c(6, 8) & !is.na(Ozone)); aq$Month <- factor(aq$Month)
t.test(Ozone ~ Month, data = aq)                    # R defaults to Welch
t.test(Ozone ~ Month, data = aq, var.equal = TRUE)  # Student t, the first row in SPSS
wilcox.test(Ozone ~ Month, data = aq, conf.int = TRUE)
effectsize::hedges_g(Ozone ~ Month, data = aq, pooled_sd = FALSE)

# Paired: sleep data, same subjects under two drugs
w <- reshape(sleep, idvar = "ID", timevar = "group", direction = "wide")
t.test(w$extra.2, w$extra.1, paired = TRUE)
wilcox.test(w$extra.2, w$extra.1, paired = TRUE, exact = FALSE)

# Several groups: one-way ANOVA, Welch ANOVA, post hoc tests, Kruskal-Wallis + Dunn
fit <- aov(weight ~ group, data = PlantGrowth); summary(fit)
oneway.test(weight ~ group, data = PlantGrowth)     # Welch ANOVA
TukeyHSD(fit)
games_howell_test(PlantGrowth, weight ~ group)      # post hoc test for unequal variances
effectsize::omega_squared(fit)
aq_all <- subset(airquality, !is.na(Ozone)); aq_all$Month <- factor(aq_all$Month)
kruskal.test(Ozone ~ Month, data = aq_all)
dunn_test(aq_all, Ozone ~ Month, p.adjust.method = "holm")

# Repeated measures: Indometh, 6 subjects at 5 time points
ind <- subset(as.data.frame(Indometh), time %in% c(0.5, 1, 2, 4, 8))
ind$Subject <- factor(as.character(ind$Subject)); ind$time <- factor(ind$time)
rm <- aov_ez(id = "Subject", dv = "conc", data = ind, within = "time")
summary(rm)                       # Mauchly test and GG/HF corrections
friedman.test(conc ~ time | Subject, data = ind)

# Categorical: chi-squared, Fisher, McNemar
tab <- table(infert$spontaneous > 0, infert$case)
chisq.test(tab)$expected          # check expected counts first
chisq.test(tab, correct = FALSE); fisher.test(tab)
mcnemar.test(matrix(c(794, 86, 150, 570), 2))

# Correlation and regression
cor.test(aq_all$Ozone, aq_all$Temp)                                   # Pearson, with 95% CI
cor.test(aq_all$Ozone, aq_all$Temp, method = "spearman", exact = FALSE)
m <- lm(log(Ozone) ~ Temp + Wind, data = aq_all); summary(m); confint(m)
g <- glm(case ~ spontaneous + induced + age, family = binomial, data = infert)
exp(cbind(OR = coef(g), confint.default(g)))
cx <- coxph(Surv(time, status) ~ age + sex, data = lung); summary(cx); cox.zph(cx)

# Multiple comparison adjustment
p <- c(0.001, 0.008, 0.012, 0.021, 0.035, 0.048, 0.09, 0.2, 0.41, 0.74)
round(cbind(p, bonf = p.adjust(p, "bonferroni"), holm = p.adjust(p, "holm"), BH = p.adjust(p, "BH")), 4)

library(survival)
write.csv(airquality, "airquality.csv", row.names = FALSE); write.csv(sleep, "sleep.csv", row.names = FALSE)
write.csv(PlantGrowth, "PlantGrowth.csv", row.names = FALSE); write.csv(as.data.frame(Indometh), "Indometh.csv", row.names = FALSE)
write.csv(infert, "infert.csv", row.names = FALSE); write.csv(lung, "lung.csv", row.names = FALSE)

# Python 3.12：pip install scipy statsmodels pingouin pandas lifelines
# Export the data from R first: write.csv(airquality, "airquality.csv", row.names = FALSE), etc.
import numpy as np, pandas as pd, pingouin as pg
from scipy import stats
import statsmodels.formula.api as smf
from statsmodels.stats.multitest import multipletests
from statsmodels.stats.contingency_tables import mcnemar

aq = pd.read_csv("airquality.csv").dropna(subset=["Ozone"])
x6, x8 = aq.Ozone[aq.Month == 6], aq.Ozone[aq.Month == 8]
print(stats.ttest_ind(x6, x8, equal_var=False))   # SciPy defaults to equal_var=True; set False explicitly for Welch
print(pg.ttest(x6, x8, correction=True))           # pingouin default 'auto' uses Welch only when sample sizes differ
print(stats.mannwhitneyu(x6, x8))

sleep = pd.read_csv("sleep.csv").pivot(index="ID", columns="group", values="extra")
print(stats.ttest_rel(sleep[2], sleep[1]))
print(stats.wilcoxon(sleep[2], sleep[1], method="approx", correction=True))  # matches the R default

pg_df = pd.read_csv("PlantGrowth.csv")
print(pg.anova(pg_df, dv="weight", between="group", detailed=True))
print(pg.welch_anova(pg_df, dv="weight", between="group"))
print(pg.pairwise_tukey(pg_df, dv="weight", between="group"))
print(pg.pairwise_gameshowell(pg_df, dv="weight", between="group"))
print(pg.kruskal(aq, dv="Ozone", between="Month"))

ind = pd.read_csv("Indometh.csv"); ind = ind[ind.time.isin([0.5, 1, 2, 4, 8])]
print(pg.rm_anova(ind, dv="conc", within="time", subject="Subject", correction=True))
print(pg.friedman(ind, dv="conc", within="time", subject="Subject"))

inf = pd.read_csv("infert.csv")
tab = pd.crosstab(inf.spontaneous > 0, inf.case).to_numpy()
chi2, p, dof, expected = stats.chi2_contingency(tab, correction=False); print(chi2, p, expected.round(1))
print(stats.fisher_exact(tab))
print(mcnemar(np.array([[794, 150], [86, 570]]), exact=False, correction=True))

print(pg.corr(aq.Ozone, aq.Temp), pg.corr(aq.Ozone, aq.Temp, method="spearman"))
print(smf.ols("np.log(Ozone) ~ Temp + Wind", data=aq).fit().summary())
logit = smf.logit("case ~ spontaneous + induced + age", data=inf).fit(disp=0)
print(np.exp(pd.concat([logit.params, logit.conf_int()], axis=1)))
from lifelines import CoxPHFitter
lung = pd.read_csv("lung.csv")[["time", "status", "age", "sex"]].dropna(); lung["status"] -= 1
print(CoxPHFitter().fit(lung, "time", "status").summary[["exp(coef)", "exp(coef) lower 95%", "exp(coef) upper 95%", "p"]])

p = [0.001, 0.008, 0.012, 0.021, 0.035, 0.048, 0.09, 0.2, 0.41, 0.74]
for m in ["bonferroni", "holm", "fdr_bh"]: print(m, multipletests(p, method=m)[1].round(4))

Principal component analysis

Principal component analysis: when to use it and how to read loadings, variance explained and the number of components

Principal component analysis (PCA) combines several correlated continuous variables into a few uncorrelated composite variables, for dimension reduction, composite indices or handling collinearity before regression. It describes the variance structure of the variables and tests no hypothesis.

R's prcomp returns eigenvectors in rotation, while the SPSS "component matrix" from the Factor menu with principal components extraction contains loadings. To get the coefficients of the component equations from the SPSS component matrix, divide each column by the square root of its eigenvalue. sklearn's explained_variance_ uses n−1 as the divisor while StandardScaler standardizes with n, so here sklearn's first eigenvalue is 2.622 versus the correlation-matrix eigenvalue 2.613; the proportions of variance explained are unaffected.

library(survival)
labs <- c("bili", "chol", "albumin", "copper", "alk.phos", "ast", "trig", "platelet", "protime")
pb <- na.omit(pbc[, labs])                      # 276 complete cases
pca <- prcomp(pb, scale. = TRUE)                # standardize: the variables have different units
summary(pca)                                    # proportion and cumulative proportion of variance
eig <- pca$sdev^2                               # eigenvalues
loadings <- sweep(pca$rotation, 2, pca$sdev, "*")   # loadings = eigenvectors x sqrt(eigenvalues), i.e. the SPSS component matrix
round(loadings[, 1:2], 3)
# Parallel analysis: keep components whose eigenvalue exceeds the 95th percentile of random data
set.seed(2026)
sim <- replicate(500, eigen(cor(matrix(rnorm(nrow(pb) * ncol(pb)), nrow(pb))))$values)
rbind(observed = round(eig, 3), random95 = round(apply(sim, 1, quantile, 0.95), 3))

import numpy as np, pandas as pd
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
pb = pd.read_csv("pbc_labs.csv")                # exported in R with write.csv(na.omit(pbc[, labs]), ...)
pca = PCA().fit(StandardScaler().fit_transform(pb))
print(pca.explained_variance_ratio_.round(4))   # same as Proportion of Variance in R summary(pca)
eig = np.linalg.eigvalsh(np.corrcoef(pb.T.to_numpy()))[::-1]   # eigenvalues of the correlation matrix, same as R and SPSS
loadings = pca.components_.T * np.sqrt(eig)     # signs may be flipped relative to R; interpret relative directions
print(pd.DataFrame(loadings[:, :2], index=pb.columns, columns=["PC1", "PC2"]).round(3))

1. 01

   Check that PCA applies

   Variables are continuous and moderately correlated (most |r|>0.3 in the correlation matrix). A common rule of thumb is a sample size of at least 5 to 10 times the number of variables. The KMO measure and Bartlett's test of sphericity come from factor analysis; they can be consulted for PCA but are not required.

2. 02

   Standardize first

   Variables with different units must be standardized (scale. = TRUE in R, i.e. a correlation-matrix PCA). On the PBC data (9 laboratory variables from survival::pbc, 276 complete cases): without standardization, the first component explains 98.3% of the variance, almost all from alkaline phosphatase, which has the largest values and a loading of −1.000; after standardization the first component explains 29.0%.

3. 03

   Decide how many components to keep

   The standardized eigenvalues are 2.613, 1.483, 0.991, 0.876 and so on. The Kaiser rule (eigenvalue >1) keeps 2; parallel analysis (comparison with the 95th percentiles of eigenvalues from random data of the same size, 1.362, 1.249, 1.165) also keeps 2. The third eigenvalue, 0.991, is close to 1; the Kaiser rule is unstable in such borderline cases and parallel analysis is more reliable. A cumulative 85% of variance is a common criterion in Chinese tutorials; here the first 2 components explain only 45.5% cumulatively, so these variables cannot be summarized by a few components.

4. 04

   Read the loadings

   A loading is the correlation between a variable and a component, equal to the eigenvector multiplied by the square root of its eigenvalue. Here bilirubin (0.836), copper (0.661), AST (0.618), triglycerides (0.561) and cholesterol (0.527) load strongly on the first component and albumin at −0.476, which can be read as severity of liver damage; the second component is mainly platelets (0.756), prothrombin time (−0.533) and cholesterol (0.513). Component signs are arbitrary: the first component has opposite signs in R and Python, so interpret the relative directions of variables.

5. 05

   Report

   Report the standardization, the number of components kept and why, the eigenvalue and variance explained for each component, and the loading matrix (usually listing variables with |loading|≥0.4). When component scores enter a regression, state how the scores were computed.

Chinese practice

Statistical advice from Chinese sources: what to follow and what to correct

The items below come from the "Xiaobai Learns Statistics" column by Feng Guoshuang, study notes on a Chinese medical statistics textbook and SPSS PCA tutorials on CSDN, checked against this page's tests and the English literature.

Hand off to an agent

Hand this workflow to an agent

Data cleaning, test selection, double computation in R and Python and writing results in SAMPL format can be handed to the science agent; judgment and interpretation remain the researcher's job.

Example instruction: "Here is my clinical data, data.xlsx. The primary outcome is HbA1c at 6 months, comparing three treatment groups, with repeated measurements at baseline and 3 months. Choose tests by data type, compute them in both R and Python and cross-check, and write the results paragraph and tables following the SAMPL guidelines."

The agent works in an isolated cloud computer: it reads the data and checks missing values, outliers and variable types; draws histograms and Q-Q plots by group and computes skewness and kurtosis; chooses the overall and post hoc tests from this page's decision table; computes them in R and Python and compares P values, effect sizes and confidence intervals item by item, tracing any difference to a default setting; and produces result tables, figures and a methods paragraph. The workspace keeps the scripts, software versions, logs and intermediate results, and the agent reviews its results adversarially. The task keeps running after you shut down your computer.

You still need to check: whether the unit of analysis matches the sampling unit; whether the primary and secondary outcomes match the protocol; the family of tests used for multiple comparison adjustment; and the clinical interpretation of the results.