V The V Lab

Common Statistical Tests in Medical Research

A Reproducible Tutorial from Study Design and Method Selection to R Implementation and Rigorous Reporting

The V Lab Home Biostatistics Foundations Medical Statistical Tests Linear Regression Logistic Regression Multilevel Modelling Causal Inference Propensity Score Matching Survival Analysis AFT and Cox PH 中文版

AudienceLearners in medicine, nursing, public health, and clinical research

Study timeAbout 150–210 minutes

PrerequisitesDescriptive statistics, confidence intervals, and basic R

About the data and their intended use All individual records, effects, and test results in this tutorial are simulated with a fixed random seed solely for teaching statistics. They contain no real patient information and must not be treated as clinical evidence that any treatment is effective or safe.

How to use this tutorial

This tutorial is not a menu for mechanically looking up tests by variable name. Every section follows the same path: define the research question and target estimand → identify the study design and data structure → check data quality and method assumptions → estimate the effect and its uncertainty → test when needed → report the result in the context of clinical importance.

Begin with the “Overview of test selection,” then study the sections that match your data structure. Code is shown by default and can be run one chunk at a time. Each example presents an effect size, a 95% confidence interval, and a p-value whenever possible.

Learning objectives

By the end of this tutorial, you should be able to:

• select common statistical methods according to outcome type, number of groups, independent or paired structure, and study design;
• distinguish among one-sample, two-independent-sample, paired-sample, multiple-independent-group, and repeated-measures questions;
• use and interpret t-tests, Wilcoxon tests, analysis of variance, and the Kruskal–Wallis test correctly;
• use chi-squared tests, Fisher’s exact test, McNemar’s test, and tests for trend correctly;
• choose Pearson, Spearman, or Kendall correlation to match the question;
• compare survival curves with the log-rank test and understand why it does not replace effect estimation;
• check key assumptions without treating a preliminary normality test as an automatic decision rule;
• address multiple comparisons and distinguish statistical significance, clinical importance, equivalence, and non-inferiority;
• report effect size, direction, 95% confidence interval, p-value, sample size, assumptions, and limitations.

1 Select a method by starting with the research question

1.1 A statistical test cannot repair the study design

A statistical test measures how compatible the data are with a particular null hypothesis under a set of model assumptions. It cannot correct selection bias, loss to follow-up, measurement error, uncontrolled confounding, incorrect temporal ordering, or a poorly framed research question.

Before calculating any p-value, write down:

1. Population and unit of analysis: Does each row represent a patient, an eye, a hospitalization, a lesion, or a visit?
2. Outcome: Is it continuous, binary, nominal multicategory, ordinal, a count, or time to event?
3. Comparison structure: Is one group compared with a reference value, are there two independent groups, are the same patients measured before and after, are there multiple groups, or are there repeated measurements?
4. Target estimand: Is it a mean difference, risk difference, risk ratio, odds ratio, correlation coefficient, or difference in survival probabilities?
5. Dependence structure: Are observations paired, clustered, repeated, censored, or survey-weighted?
6. Clinical threshold: How large must a difference be to matter clinically or in public health?

The most common mistake at the starting line “Is my variable normally distributed, and therefore which test should I use?” is usually not the first question. First identify the comparison, independence, pairing, outcome scale, and target estimand. These features are often more important than whether a marginal distribution is perfectly normal.

1.2 Quick selection matrix

Research question and structure	Common primary method	Main effect or estimand	Important reminder
One continuous sample versus a reference value	One-sample t-test	Difference between the mean and reference value	Inference targets the mean; check outliers and the sampling mechanism
Continuous outcome in two independent groups	Welch two-sample t-test	Mean difference	Does not require equal group variances by default
Continuous outcome before and after in the same individuals	Paired t-test	Mean of individual differences	Analyze the differences; do not treat observations as independent groups
Ordinal or markedly skewed outcome in two independent groups	Wilcoxon rank-sum test	Distributional location / probability of superiority	Usually not a “test of medians”
Paired ordinal outcome or skewed paired differences	Wilcoxon signed-rank test	Location of the difference distribution	Differences should be approximately symmetric; describe how zeros were handled
Continuous outcome in three or more independent groups	One-way / Welch ANOVA	Differences among group means	Follow the omnibus test with prespecified or multiplicity-adjusted comparisons
Ordinal or skewed outcome in three or more independent groups	Kruskal–Wallis test	Differences in rank distributions	A significant result requires adjusted pairwise comparisons
Three or more repeated measurements	Repeated-measures model / Friedman test	Differences across time or conditions	Preserve within-patient correlation
Two independent categorical variables	Pearson chi-squared test	Association between proportions	Consider Fisher’s exact test when expected counts are small
Sparse 2 x 2 table from two independent groups	Fisher’s exact test	Odds ratio and exact inference	The odds ratio is not the risk ratio; understand the design assumptions about fixed margins
Binary outcome before and after in the same individuals	McNemar’s test	Directional difference among discordant pairs	Uses information only from discordant pairs
Linear association between two continuous variables	Pearson correlation	Correlation coefficient $r$	Inspect the scatterplot, outliers, and independence
Monotonic nonlinear or ordinal association	Spearman or Kendall correlation	Rank correlation	Not a causal effect and not a measure of agreement
Censored event times in two groups	Kaplan–Meier + log-rank	Difference between survival curves	Also report survival probabilities, time points, and numbers at risk

Core principle: A test asks, “Are the data compatible with the null hypothesis?” An effect estimate asks, “How large is the difference, and in which direction?” Medical reports usually need the latter and its confidence interval, not merely an asterisk for statistical significance.

1.3 Distinguish independent, paired, and clustered data first

• Independent samples: The groups consist of different patients, and each patient contributes one independent observation.
• Paired samples: Observations correspond explicitly, as in before-and-after measurements from the same patient, measurements of the left and right eyes, or matched cases and controls.
• Repeated measures: Each patient contributes observations at three or more time points.
• Clustered data: Patients are nested within hospitals, teeth within patients, or multiple lesions within one person.

Ordinary t-tests, chi-squared tests, and rank tests generally assume that units of analysis are independent. Treating multiple observations from the same patient as independent underestimates standard errors and exaggerates the evidence. More complex repeated or clustered data generally require mixed-effects models, generalized estimating equations, cluster-robust standard errors, or randomization methods matched to the design.

2 Continuous outcomes: one- and two-group comparisons

2.1 One-sample t-test

The one-sample t-test compares a population mean with a prespecified reference value $\mu_0$:

$$H_0:\mu=\mu_0,\qquad t=\frac{\bar x-\mu_0}{s/\sqrt{n}}.$$

Example question: Is the mean reduction in systolic blood pressure at 12 weeks in these patients different from the prespecified value of 8 mmHg? The reference value should come from the protocol, a historical standard, or the clinical question; it must not be chosen after looking at the data.

new_treatment_reduction <- subset(
  trial_data,
  treatment == "New treatment"
)$sbp_reduction

one_sample_result <- t.test(
  new_treatment_reduction,
  mu = 8,
  conf.level = 0.95
)

data.frame(
  `Sample size` = length(new_treatment_reduction),
  `Sample mean` = mean(new_treatment_reduction),
  `Difference from 8` = mean(new_treatment_reduction) - 8,
  `Lower 95% CI` = unname(one_sample_result$conf.int[1] - 8),
  `Upper 95% CI` = unname(one_sample_result$conf.int[2] - 8),
  `p-value` = one_sample_result$p.value,
  check.names = FALSE
) |>
  knitr::kable(
    digits = 3,
    caption = "Mean systolic blood pressure reduction in the new-treatment group versus 8 mmHg"
  )

Mean systolic blood pressure reduction in the new-treatment group versus 8 mmHg

Sample size	Sample mean	Difference from 8	Lower 95% CI	Upper 95% CI	p-value
160	11.6	3.62	2.23	5.01	0

Key conditions include independent observations; a meaningful mean for the outcome; no serious erroneous values that dominate the mean and standard error; and a sampling process consistent with the target of inference. With small samples, the outcome distribution should also be approximately normal. With larger samples, the sampling distribution of the mean is generally more robust, but extreme skewness and outliers still require careful attention.

2.2 Two independent samples: Prefer the Welch t-test

When the target is the mean difference between two independent groups, the Welch t-test is usually a sensible default. It does not force the two population variances to be equal and typically sacrifices little efficiency when they truly are equal.

welch_result <- t.test(
  sbp_reduction ~ treatment,
  data = trial_data,
  var.equal = FALSE
)

group_summary <- do.call(
  rbind,
  lapply(split(trial_data$sbp_reduction, trial_data$treatment), function(x) {
    c(n = length(x), mean = mean(x), sd = sd(x))
  })
)
group_summary <- data.frame(
  treatment = rownames(group_summary),
  group_summary,
  row.names = NULL
)

knitr::kable(
  group_summary,
  digits = 2,
  col.names = c("Treatment group", "n", "Mean", "Standard deviation"),
  caption = "Descriptive statistics for systolic blood pressure reduction by group"
)

Descriptive statistics for systolic blood pressure reduction by group

Treatment group	n	Mean	Standard deviation
Standard care	160	7.07	8.24
New treatment	160	11.62	8.90

data.frame(
  Comparison = "Standard care − New treatment",
  `Mean difference` = unname(welch_result$estimate[1] - welch_result$estimate[2]),
  `Lower 95% CI` = unname(welch_result$conf.int[1]),
  `Upper 95% CI` = unname(welch_result$conf.int[2]),
  `Degrees of freedom` = unname(welch_result$parameter),
  `p-value` = welch_result$p.value,
  check.names = FALSE
) |>
  knitr::kable(
    digits = 3,
    caption = "Welch two-sample t-test"
  )

Welch two-sample t-test

Comparison	Mean difference	Lower 95% CI	Upper 95% CI	Degrees of freedom	p-value
Standard care − New treatment	-4.55	-6.43	-2.66	316	0

In t.test(), the first factor level is subtracted from the second, so the difference here is “Standard care − New treatment.” If a positive value should represent a larger reduction under the new treatment, explicitly reverse the direction in the reporting table; never copy output without checking the coding.

2.2.1 Effect size: Prioritize the mean difference and supplement it with a standardized difference

The mean difference in the original measurement units is usually the easiest quantity for judging clinical importance. A standardized mean difference can support comparisons across scales, but it is affected by the amount of variability in the study population.

standard_values <- subset(
  trial_data,
  treatment == "Standard care"
)$sbp_reduction
new_values <- subset(
  trial_data,
  treatment == "New treatment"
)$sbp_reduction

n0 <- length(standard_values)
n1 <- length(new_values)
pooled_sd <- sqrt(
  ((n0 - 1) * var(standard_values) + (n1 - 1) * var(new_values)) /
    (n0 + n1 - 2)
)
cohens_d <- (mean(new_values) - mean(standard_values)) / pooled_sd
small_sample_correction <- 1 - 3 / (4 * (n0 + n1) - 9)
hedges_g <- small_sample_correction * cohens_d

data.frame(
  Measure = c("Mean difference: new treatment − standard care (mmHg)", "Hedges' g"),
  Estimate = c(mean(new_values) - mean(standard_values), hedges_g)
) |>
  knitr::kable(digits = 3, caption = "Effect sizes in original and standardized units")

Effect sizes in original and standardized units

Measure	Estimate
Mean difference: new treatment − standard care (mmHg)	4.548
Hedges’ g	0.529

Recommended reporting template State how many participants were included in the new-treatment and standard-care groups; report the mean (SD) in each group; then report the mean difference in the prespecified direction, its 95% CI, and the p-value from the Welch t-test; finally explain what the interval means relative to the minimum clinically important difference.

2.3 Paired t-test

The paired t-test first converts each patient’s measurements into one difference, “after treatment − before treatment,” and then tests whether the mean difference is zero. It assumes that differences from distinct patients are independent and focuses on the distribution of the differences; it does not separately require the before- and after-treatment measurements to be normal.

paired_result <- t.test(
  trial_data$followup_sbp,
  trial_data$baseline_sbp,
  paired = TRUE
)

within_person_change <- trial_data$followup_sbp - trial_data$baseline_sbp

data.frame(
  `Sample size` = length(within_person_change),
  `Mean change` = mean(within_person_change),
  `SD of change` = sd(within_person_change),
  `Lower 95% CI` = unname(paired_result$conf.int[1]),
  `Upper 95% CI` = unname(paired_result$conf.int[2]),
  `p-value` = paired_result$p.value,
  check.names = FALSE
) |>
  knitr::kable(
    digits = 3,
    caption = "Paired t-test of follow-up minus baseline values among all participants"
  )

Paired t-test of follow-up minus baseline values among all participants

Sample size	Mean change	SD of change	Lower 95% CI	Upper 95% CI	p-value
320	-9.35	8.86	-10.3	-8.37	0

plot_ids <- seq_len(50)
matplot(
  x = c(1, 2),
  y = t(as.matrix(trial_data[plot_ids, c("baseline_sbp", "followup_sbp")])),
  type = "l",
  lty = 1,
  col = grDevices::adjustcolor(palette_test[["gray"]], alpha.f = 0.32),
  xaxt = "n",
  xlab = "Time point",
  ylab = "Systolic blood pressure (mmHg)"
)
axis(1, at = c(1, 2), labels = c("Baseline", "Week 12"))
points(
  c(1, 2),
  c(mean(trial_data$baseline_sbp), mean(trial_data$followup_sbp)),
  type = "b",
  pch = 19,
  lwd = 3,
  col = palette_test[["vermillion"]]
)

Baseline and 12-week systolic blood pressure for each participant. Thin lines connect measurements from the same person.

Different within-group significance does not imply a difference in change between groups A statistically significant before-and-after comparison in one treatment group and a nonsignificant comparison in the other do not show that the groups changed differently. Directly compare individual changes between groups or, more commonly, compare follow-up outcomes with a baseline-adjusted model while retaining randomized treatment assignment.

2.4 Wilcoxon rank-sum test: Two independent groups

The Wilcoxon rank-sum test (Mann–Whitney U test) can be used when the outcome is ordinal or when a continuous outcome is strongly skewed and a rank-based comparison matches the question. It uses the ranks after pooling observations from the two groups.

rank_sum_result <- wilcox.test(
  crp_week12 ~ treatment,
  data = trial_data,
  exact = FALSE,
  conf.int = TRUE
)

crp_summary <- do.call(
  rbind,
  lapply(split(trial_data$crp_week12, trial_data$treatment), function(x) {
    c(
      n = length(x),
      median = median(x),
      q1 = unname(quantile(x, 0.25)),
      q3 = unname(quantile(x, 0.75))
    )
  })
)

knitr::kable(
  data.frame(Treatment = rownames(crp_summary), crp_summary, row.names = NULL),
  digits = 2,
  col.names = c("Treatment group", "n", "Median", "Q1", "Q3"),
  caption = "Median and quartiles of CRP at 12 weeks"
)

Median and quartiles of CRP at 12 weeks

Treatment group	n	Median	Q1	Q3
Standard care	160	3.09	1.96	5.53
New treatment	160	2.88	1.62	4.69

data.frame(
  `W statistic` = unname(rank_sum_result$statistic),
  `Location-shift estimate` = unname(rank_sum_result$estimate),
  `Lower 95% CI` = unname(rank_sum_result$conf.int[1]),
  `Upper 95% CI` = unname(rank_sum_result$conf.int[2]),
  `p-value` = rank_sum_result$p.value,
  check.names = FALSE
) |>
  knitr::kable(
    digits = 3,
    caption = "Wilcoxon rank-sum test and Hodges–Lehmann-type location shift"
  )

Wilcoxon rank-sum test and Hodges–Lehmann-type location shift

W statistic	Location-shift estimate	Lower 95% CI	Upper 95% CI	p-value
14232	0.36	-0.06	0.81	0.084

If the two distributions have the same shape and their main difference can be described as an overall shift, the result can be interpreted as a location shift. More generally, the test evaluates whether the groups have the same rank distribution; it must not automatically be described as a test of equal medians. Even when using a rank test, plot and report the distribution in each group.

boxplot(
  crp_week12 ~ treatment,
  data = trial_data,
  col = c(palette_test[["sky"]], palette_test[["orange"]]),
  border = palette_test[["navy"]],
  ylab = "CRP (mg/L)",
  xlab = "Treatment group",
  outline = FALSE
)
stripchart(
  crp_week12 ~ treatment,
  data = trial_data,
  vertical = TRUE,
  method = "jitter",
  pch = 16,
  cex = 0.55,
  col = grDevices::adjustcolor(palette_test[["navy"]], alpha.f = 0.32),
  add = TRUE
)

Right-skewed distributions of 12-week CRP in the two groups.

2.5 Wilcoxon signed-rank and sign tests: Paired data

The Wilcoxon signed-rank test is a rank-based analysis of paired differences. In addition to independent pairs, it generally requires the difference distribution to be approximately symmetric around its center. If symmetry is difficult to justify, an exact sign test that uses only the signs can be used, but it discards information about the magnitudes of the differences.

signed_rank_result <- wilcox.test(
  trial_data$followup_sbp,
  trial_data$baseline_sbp,
  paired = TRUE,
  exact = FALSE,
  conf.int = TRUE
)

nonzero_change <- within_person_change[within_person_change != 0]
negative_changes <- sum(nonzero_change < 0)
sign_test_result <- binom.test(
  negative_changes,
  length(nonzero_change),
  p = 0.5
)

data.frame(
  Method = c("Wilcoxon signed-rank test", "Exact sign test"),
  `Statistic or count` = c(
    unname(signed_rank_result$statistic),
    negative_changes
  ),
  `p-value` = c(signed_rank_result$p.value, sign_test_result$p.value),
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "Two nonparametric tests for paired data")

Two nonparametric tests for paired data

Method	Statistic or count	p-value
Wilcoxon signed-rank test	3260	0
Exact sign test	276	0

Check your understanding: t-test or rank test?

A randomized trial compares hospital length of stay between two groups. Many patients stay for 1–3 days, while a few stay for more than 60 days. The researchers care most about mean use of inpatient resources. Should they automatically switch to the Wilcoxon test because the distribution is skewed?

Answer: No. If the target estimand is mean hospital length of stay, the Wilcoxon test does not answer the same question. First verify whether the extreme values are genuine, then consider robust or bootstrap intervals for the mean difference, permutation methods, models suited to right-skewed outcomes, and reporting both the mean and the distribution. The method must match the target estimand.

3 Continuous or Ordinal Outcomes Across Three or More Groups

3.1 One-Way ANOVA and Welch ANOVA

The omnibus null hypothesis for a one-way analysis of variance is that all population group means are equal. It answers whether at least one group mean differs, but it does not automatically identify which groups differ. Conventional ANOVA assumes independent errors, a reasonable mean model, approximately normal within-group residuals, and equal population variances; Welch ANOVA relaxes the equal-variance assumption.

dose_summary <- do.call(
  rbind,
  lapply(split(dose_data$hemoglobin_change, dose_data$dose), function(x) {
    c(n = length(x), mean = mean(x), sd = sd(x))
  })
)

knitr::kable(
  data.frame(Group = rownames(dose_summary), dose_summary, row.names = NULL),
  digits = 3,
  col.names = c("Group", "n", "Mean", "Standard deviation"),
  caption = "Descriptive statistics for hemoglobin change in three groups"
)

Descriptive statistics for hemoglobin change in three groups

Group	n	Mean	Standard deviation
Placebo	60	0.075	0.773
Low dose	60	0.695	0.847
High dose	60	1.219	1.146

traditional_anova <- aov(hemoglobin_change ~ dose, data = dose_data)
welch_anova <- oneway.test(
  hemoglobin_change ~ dose,
  data = dose_data,
  var.equal = FALSE
)

anova_table <- summary(traditional_anova)[[1]]
ss_between <- anova_table["dose", "Sum Sq"]
ss_total <- sum(anova_table[, "Sum Sq"])
ms_error <- anova_table["Residuals", "Mean Sq"]
k_groups <- nlevels(dose_data$dose)
omega_squared <- max(
  0,
  (ss_between - (k_groups - 1) * ms_error) / (ss_total + ms_error)
)

data.frame(
  Method = c("Conventional one-way ANOVA", "Welch ANOVA"),
  Statistic = c(
    anova_table["dose", "F value"],
    unname(welch_anova$statistic)
  ),
  `Numerator degrees of freedom` = c(
    anova_table["dose", "Df"],
    unname(welch_anova$parameter[1])
  ),
  `Denominator degrees of freedom` = c(
    anova_table["Residuals", "Df"],
    unname(welch_anova$parameter[2])
  ),
  `p value` = c(
    anova_table["dose", "Pr(>F)"],
    welch_anova$p.value
  ),
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "Conventional ANOVA and Welch ANOVA")

Conventional ANOVA and Welch ANOVA

Method	Statistic	Numerator degrees of freedom	Denominator degrees of freedom	p value
Conventional one-way ANOVA	22.5	2	177	0
Welch ANOVA	22.3	2	115	0

data.frame(
  `Effect size` = "omega-squared",
  Estimate = omega_squared,
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "Omnibus effect size for conventional one-way ANOVA")

Omnibus effect size for conventional one-way ANOVA

Effect size	Estimate
omega-squared	0.193

$\omega^2$ approximately describes the proportion of total outcome variation attributable to between-group differences, but fixed “small,” “medium,” and “large” thresholds should not be applied without considering the research context. It also cannot replace estimates of specific between-group mean differences.

3.1.1 Comparisons After the Omnibus Test

If the protocol prespecifies “Low dose versus Placebo” and “High dose versus Placebo,” these contrasts align more closely with the research question than first performing an omnibus test and then exhaustively examining every combination only if it is significant. If all pairwise comparisons are genuinely of interest, the Tukey method can be used; when variances are unequal, pairwise Welch t-tests with a multiplicity adjustment are an option.

tukey_result <- TukeyHSD(traditional_anova, "dose")$dose
knitr::kable(
  data.frame(Comparison = rownames(tukey_result), tukey_result, row.names = NULL),
  digits = 3,
  col.names = c(
    "Comparison",
    "Mean difference",
    "Lower simultaneous CI",
    "Upper simultaneous CI",
    "Adjusted p value"
  ),
  caption = "Tukey comparisons for all pairs"
)

Tukey comparisons for all pairs

Comparison	Mean difference	Lower simultaneous CI	Upper simultaneous CI	Adjusted p value
Low dose-Placebo	0.620	0.216	1.024	0.001
High dose-Placebo	1.144	0.740	1.548	0.000
High dose-Low dose	0.524	0.120	0.928	0.007

pairwise_welch <- pairwise.t.test(
  dose_data$hemoglobin_change,
  dose_data$dose,
  p.adjust.method = "holm",
  pool.sd = FALSE
)
pairwise_welch$p.value |>
  knitr::kable(
    digits = 3,
    caption = "Holm-adjusted p values for pairwise Welch t-tests"
  )

Holm-adjusted p values for pairwise Welch t-tests

	Placebo	Low dose
Low dose	0	NA
High dose	0	0.005

Note that the pairwise.t.test() table above provides adjusted p values but not simultaneous confidence intervals. If multiple contrasts are central to the conclusions of a formal report, use a post-model contrast tool that produces simultaneous intervals consistent with the chosen adjustment method, rather than presenting adjusted p values alongside unadjusted intervals.

3.2 Kruskal–Wallis Test

The Kruskal–Wallis test extends the Wilcoxon rank-sum method to three or more independent groups. Its null hypothesis concerns identical rank distributions across groups; it is appropriate to simplify the result as a comparison of locations only when the distributions have similar shapes and their differences can reasonably be viewed as location shifts.

kruskal_result <- kruskal.test(
  inflammation_score ~ dose,
  data = dose_data
)

h_statistic <- unname(kruskal_result$statistic)
epsilon_squared <- max(
  0,
  (h_statistic - k_groups + 1) / (nrow(dose_data) - k_groups)
)

rank_summary <- do.call(
  rbind,
  lapply(split(dose_data$inflammation_score, dose_data$dose), function(x) {
    c(
      n = length(x),
      median = median(x),
      q1 = unname(quantile(x, 0.25)),
      q3 = unname(quantile(x, 0.75))
    )
  })
)

knitr::kable(
  data.frame(Group = rownames(rank_summary), rank_summary, row.names = NULL),
  digits = 2,
  col.names = c("Group", "n", "Median", "Q1", "Q3"),
  caption = "Descriptive statistics for inflammation scores in three groups"
)

Descriptive statistics for inflammation scores in three groups

Group	n	Median	Q1	Q3
Placebo	60	2.96	1.40	5.43
Low dose	60	2.51	1.65	3.70
High dose	60	2.58	1.60	4.14

data.frame(
  `H statistic` = h_statistic,
  `Degrees of freedom` = unname(kruskal_result$parameter),
  `p value` = kruskal_result$p.value,
  `Rank epsilon-squared` = epsilon_squared,
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "Kruskal–Wallis test")

Kruskal–Wallis test

H statistic	Degrees of freedom	p value	Rank epsilon-squared
1.46	2	0.483	0

pairwise_rank <- pairwise.wilcox.test(
  dose_data$inflammation_score,
  dose_data$dose,
  p.adjust.method = "holm",
  exact = FALSE
)

pairwise_rank$p.value |>
  knitr::kable(
    digits = 3,
    caption = "Holm-adjusted p values for pairwise Wilcoxon rank-sum tests"
  )

Holm-adjusted p values for pairwise Wilcoxon rank-sum tests

	Placebo	Low dose
Low dose	0.839	NA
High dose	0.839	0.967

An omnibus result does not mean that every pair differs A small p value from ANOVA or the Kruskal–Wallis test only indicates evidence against its respective omnibus null hypothesis; it does not show that every pair differs. Specific between-group differences must be estimated using prespecified contrasts or comparisons with appropriate multiplicity control.

3.3 Three or More Repeated Measurements

Observations from the same patient at multiple time points are usually correlated. Repeated-measures ANOVA can be used for complete, balanced continuous data; the Friedman test is a rank-based method for complete blocks. In real longitudinal studies with unequally spaced times, missing follow-up, time-varying covariates, or between-group comparisons, mixed-effects models or generalized estimating equations are usually more appropriate.

repeated_summary <- do.call(
  rbind,
  lapply(split(repeated_data$pain_score, repeated_data$time), function(x) {
    c(n = length(x), mean = mean(x), sd = sd(x), median = median(x))
  })
)

knitr::kable(
  data.frame(`Time point` = rownames(repeated_summary), repeated_summary, row.names = NULL),
  digits = 2,
  col.names = c("Time point", "n", "Mean", "Standard deviation", "Median"),
  caption = "Descriptive statistics for pain scores at three time points"
)

Descriptive statistics for pain scores at three time points

Time point	n	Mean	Standard deviation	Median
Baseline	90	6.53	1.59	6.55
Week 4	90	4.97	1.59	5.05
Week 12	90	4.08	1.48	4.10

repeated_anova <- aov(
  pain_score ~ time + Error(patient_id / time),
  data = repeated_data
)
summary(repeated_anova)

## 
## Error: patient_id
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 89    541    6.08               
## 
## Error: patient_id:time
##            Df Sum Sq Mean Sq F value Pr(>F)    
## time        2    276   137.9     235 <2e-16 ***
## Residuals 178    104     0.6                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

For more than two time points, the conventional repeated-measures ANOVA above also involves the sphericity assumption. The simplified output from aov() does not automatically provide Mauchly’s test or a Greenhouse–Geisser correction; when formal inference is needed, use a method that explicitly handles sphericity or the covariance structure.

friedman_result <- friedman.test(
  pain_score ~ time | patient_id,
  data = repeated_data
)

data.frame(
  `Chi-square statistic` = unname(friedman_result$statistic),
  `Degrees of freedom` = unname(friedman_result$parameter),
  `p value` = friedman_result$p.value,
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "Friedman rank test for repeated measurements")

Friedman rank test for repeated measurements

Chi-square statistic	Degrees of freedom	p value
134	2	0

To explore which time points differ, adjust the p values from paired comparisons:

time_pairs <- combn(colnames(pain_matrix), 2, simplify = FALSE)
raw_p <- vapply(time_pairs, function(pair) {
  wilcox.test(
    pain_matrix[, pair[1]],
    pain_matrix[, pair[2]],
    paired = TRUE,
    exact = FALSE
  )$p.value
}, numeric(1))

repeated_pairwise <- data.frame(
  Comparison = vapply(time_pairs, paste, collapse = " vs ", FUN.VALUE = character(1)),
  `Unadjusted p value` = raw_p,
  `Holm-adjusted p value` = p.adjust(raw_p, method = "holm"),
  check.names = FALSE
)

knitr::kable(
  repeated_pairwise,
  digits = 3,
  caption = "Multiplicity adjustment for paired rank tests across three time points"
)

Multiplicity adjustment for paired rank tests across three time points

Comparison	Unadjusted p value	Holm-adjusted p value
Baseline vs Week 4	0	0
Baseline vs Week 12	0	0
Week 4 vs Week 12	0	0

When a randomized trial includes a baseline measurement A two-group randomized trial generally should not compare only the within-group pre–post p values. For a continuous follow-up outcome, a prespecified baseline-adjusted ANCOVA can directly compare follow-up means between groups, improve precision, and account for chance baseline imbalance. If there are multiple follow-up time points, consider a longitudinal model.

4 Categorical Outcomes and Proportions

4.1 One Proportion: Exact Binomial and Score Approximation

If 62 of 80 patients achieve a treatment response, we can test whether the true response proportion equals a prespecified value of 0.70. binom.test() uses the exact binomial distribution, whereas prop.test() provides a large-sample score-type approximation; their intervals and p values may differ slightly.

x_response <- 62
n_response <- 80
exact_proportion <- binom.test(x_response, n_response, p = 0.70)
score_proportion <- prop.test(
  x_response,
  n_response,
  p = 0.70,
  correct = FALSE
)

data.frame(
  Method = c("Exact binomial", "Score approximation"),
  `Observed proportion` = x_response / n_response,
  `Lower CI` = c(exact_proportion$conf.int[1], score_proportion$conf.int[1]),
  `Upper CI` = c(exact_proportion$conf.int[2], score_proportion$conf.int[2]),
  `p value` = c(exact_proportion$p.value, score_proportion$p.value),
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "Two inferential methods for a single response proportion")

Two inferential methods for a single response proportion

Method	Observed proportion	Lower CI	Upper CI	p value
Exact binomial	0.775	0.668	0.861	0.179
Score approximation	0.775	0.672	0.853	0.143

4.2 Two Independent Proportions: Chi-Square, Proportion Test, and Effect Measures

For a binary outcome, a chi-square p value alone is not enough. In randomized trials or cohort studies, the risk difference (RD) and risk ratio (RR) are usually more direct; in case–control studies, where sampling is based on outcome status, the odds ratio (OR) is usually the primary measure of association.

response_table <- table(trial_data$treatment, trial_data$response)
response_table |>
  addmargins() |>
  knitr::kable(caption = "2×2 table of treatment group and response at 12 weeks")

2×2 table of treatment group and response at 12 weeks

	No	Yes	Sum
Standard care	67	93	160
New treatment	49	111	160
Sum	116	204	320

response_risks <- prop.table(response_table, margin = 1)[, "Yes"]
data.frame(
  `Treatment group` = names(response_risks),
  `Response risk` = unname(response_risks),
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "Observed response risk in each group")

Observed response risk in each group

Treatment group	Response risk
Standard care	0.581
New treatment	0.694

chi_result <- chisq.test(response_table, correct = FALSE)
prop_result <- prop.test(
  x = response_table[, "Yes"],
  n = rowSums(response_table),
  correct = FALSE
)

data.frame(
  Method = c("Pearson chi-square test", "Two-sample proportion test"),
  Statistic = c(unname(chi_result$statistic), unname(prop_result$statistic)),
  `Degrees of freedom` = c(unname(chi_result$parameter), unname(prop_result$parameter)),
  `p value` = c(chi_result$p.value, prop_result$p.value),
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "Omnibus tests for two independent proportions")

Omnibus tests for two independent proportions

Method	Statistic	Degrees of freedom	p value
Pearson chi-square test	4.38	1	0.036
Two-sample proportion test	4.38	1	0.036

risk_measures(response_table) |>
  knitr::kable(
    digits = 3,
    caption = "Absolute and relative effects for response (large-sample Wald intervals)"
  )

Absolute and relative effects for response (large-sample Wald intervals)

Measure	Estimate	Lower 95% CI	Upper 95% CI
Risk under new treatment	0.694	NA	NA
Risk under standard care	0.581	NA	NA
Risk difference	0.112	0.008	0.217
Risk ratio	1.194	1.010	1.411
Odds ratio	1.632	1.030	2.585

The null value is 0 for the risk difference and 1 for both the risk ratio and odds ratio. The manually calculated intervals above are for illustration; sparse data, clustered data, or adjusted analyses require methods matched to the design and model.

4.2.1 Expected Counts and Cramér’s V

The chi-square approximation requires expected counts large enough to support the asymptotic distribution; it should not be assessed using a fixed slogan such as “the total sample size is less than 30.” First inspect the expected counts and the extent of sparsity:

knitr::kable(
  chi_result$expected,
  digits = 2,
  caption = "Expected counts under the null hypothesis of independence"
)

Expected counts under the null hypothesis of independence

	No	Yes
Standard care	58	102
New treatment	58	102

cramers_v <- sqrt(
  unname(chi_result$statistic) /
    (sum(response_table) * min(nrow(response_table) - 1, ncol(response_table) - 1))
)
data.frame(`Effect size` = "Cramér's V", Estimate = cramers_v, check.names = FALSE) |>
  knitr::kable(digits = 3, caption = "Standardized effect size for a categorical association")

Standardized effect size for a categorical association

Effect size	Estimate
Cramér’s V	0.117

4.3 Fisher’s Exact Test: A Sparse 2×2 Table

Rare serious adverse events often produce small cell counts. Fisher’s exact test avoids relying on the large-sample chi-square approximation, but the original denominators, absolute risks, and effect intervals should still be reported.

rare_event_table <- matrix(
  c(1, 39, 7, 33),
  nrow = 2,
  byrow = TRUE,
  dimnames = list(
    Treatment = c("New treatment", "Standard care"),
    `Serious adverse event` = c("Yes", "No")
  )
)

fisher_result <- fisher.test(rare_event_table)

knitr::kable(
  addmargins(rare_event_table),
  caption = "Sparse 2×2 table of serious adverse events"
)

Sparse 2×2 table of serious adverse events

	Yes	No	Sum
New treatment	1	39	40
Standard care	7	33	40
Sum	8	72	80

data.frame(
  `Conditional odds ratio estimate` = unname(fisher_result$estimate),
  `Lower CI` = unname(fisher_result$conf.int[1]),
  `Upper CI` = unname(fisher_result$conf.int[2]),
  `p value` = fisher_result$p.value,
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "Fisher's exact test")

Fisher’s exact test

Conditional odds ratio estimate	Lower CI	Upper CI	p value
0.124	0.003	1.04	0.057

The conditional maximum-likelihood OR estimate returned by fisher.test() may differ slightly from the simple cross-product $ad/bc$. A zero cell can also make the simple OR equal to 0 or infinity. Do not add 0.5 arbitrarily merely to obtain a finite value unless that analytical method was prespecified and has a sound justification.

4.4 Paired Binary Outcomes: McNemar’s Test

Whether the same patient has symptoms before and after an intervention is a paired binary outcome. McNemar’s test uses only the two types of discordant pairs: symptomatic before and asymptomatic after, and asymptomatic before and symptomatic after.

knitr::kable(
  addmargins(paired_symptom),
  caption = "Symptom status before and after intervention in the same patients"
)

Symptom status before and after intervention in the same patients

	Asymptomatic	Symptomatic	Sum
Asymptomatic	49	6	55
Symptomatic	45	40	85
Sum	94	46	140

mcnemar_result <- mcnemar.test(paired_symptom, correct = TRUE)
b <- paired_symptom["Asymptomatic", "Symptomatic"]
c <- paired_symptom["Symptomatic", "Asymptomatic"]
exact_mcnemar <- binom.test(b, b + c, p = 0.5)

data.frame(
  Method = c(
    "McNemar chi-square approximation (continuity corrected)",
    "Exact binomial test of discordant pairs"
  ),
  `Statistic or count` = c(unname(mcnemar_result$statistic), b),
  `p value` = c(mcnemar_result$p.value, exact_mcnemar$p.value),
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "Tests for a paired binary outcome")

Tests for a paired binary outcome

Method	Statistic or count	p value
McNemar chi-square approximation (continuity corrected)	28.3	0
Exact binomial test of discordant pairs	6.0	0

before_risk <- mean(symptom_before)
after_risk <- mean(symptom_after)
paired_risk_difference <- after_risk - before_risk
matched_or <- b / c

data.frame(
  Measure = c(
    "Symptom proportion before intervention",
    "Symptom proportion after intervention",
    "After − before risk difference",
    "Discordant-pair odds ratio"
  ),
  Estimate = c(before_risk, after_risk, paired_risk_difference, matched_or)
) |>
  knitr::kable(digits = 3, caption = "Descriptive effects for paired binary data")

Descriptive effects for paired binary data

Measure	Estimate
Symptom proportion before intervention	0.607
Symptom proportion after intervention	0.329
After − before risk difference	-0.279
Discordant-pair odds ratio	0.133

When discordant pairs are few, prioritize the exact result and transparently present the wide interval. The paired OR describes only the direction of the discordant pairs and should not be confused with an OR from an independent 2×2 table.

4.5 Trend in Proportions Across an Ordered Exposure

When dose levels have a scientifically prespecified order, prop.trend.test() can test whether the proportions follow a linear trend. It is more focused than an unordered chi-square test, but it cannot detect every nonlinear pattern.

responders_by_dose <- c(21, 31, 43)
total_by_dose <- c(60, 60, 60)
trend_result <- prop.trend.test(
  responders_by_dose,
  total_by_dose,
  score = c(0, 1, 2)
)

data.frame(
  Dose = c("Placebo", "Low dose", "High dose"),
  Responders = responders_by_dose,
  Total = total_by_dose,
  `Response proportion` = responders_by_dose / total_by_dose,
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "Response proportions by ordered dose group")

Response proportions by ordered dose group

Dose	Responders	Total	Response proportion
Placebo	21	60	0.350
Low dose	31	60	0.517
High dose	43	60	0.717

data.frame(
  `Trend chi-square` = unname(trend_result$statistic),
  `Degrees of freedom` = unname(trend_result$parameter),
  `p value` = trend_result$p.value,
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "Cochran–Armitage-type test for a trend in proportions")

Cochran–Armitage-type test for a trend in proportions

Trend chi-square	Degrees of freedom	p value
16.2	1	0

4.6 Stratified 2×2 Tables: Mantel–Haenszel Test

If the 2×2 table of treatment and outcome is stratified by center or by a prespecified confounder, the Mantel–Haenszel method can estimate a common OR and test the conditional association. It assumes that the stratum-specific ORs are sufficiently similar; if effects clearly differ across strata, a single common OR will conceal effect modification.

mh_array <- array(
  c(
    28, 22, 19, 31,
    34, 16, 25, 25,
    19, 31, 13, 37
  ),
  dim = c(2, 2, 3),
  dimnames = list(
    Treatment = c("New treatment", "Standard care"),
    Response = c("Yes", "No"),
    Center = c("Center A", "Center B", "Center C")
  )
)

mh_result <- mantelhaen.test(mh_array)

data.frame(
  `Common odds ratio` = unname(mh_result$estimate),
  `Lower CI` = unname(mh_result$conf.int[1]),
  `Upper CI` = unname(mh_result$conf.int[2]),
  `p value` = mh_result$p.value,
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "Mantel–Haenszel analysis stratified by study center")

Mantel–Haenszel analysis stratified by study center

Common odds ratio	Lower CI	Upper CI	p value
1.98	1.24	3.18	0.007

Check your understanding: chi-square or McNemar?

A study gives every patient two rapid tests in sequence and compares their positivity rates. Should the analysis use Pearson’s chi-square test or McNemar’s test?

Answer: McNemar’s test. The two test results come from the same patient and are paired. An ordinary Pearson chi-square test would incorrectly treat the two sets of results as independent. If the research objective is to compare sensitivity, the analysis must also be restricted to patients confirmed by the reference standard while retaining the paired structure.

5 Event Counts and Incidence Rates

5.1 Poisson Rate Test

When participants have different follow-up times, simply comparing the proportion who ever experienced an event may discard useful information. When the event rate is the target, report the number of events, total person-time at risk, and the incidence rate in clearly specified units.

Suppose 18 events occurred over 410 person-years in the new treatment group and 31 events occurred over 395 person-years in the standard care group:

event_counts <- c("New treatment" = 18, "Standard care" = 31)
person_years <- c("New treatment" = 410, "Standard care" = 395)
rate_result <- poisson.test(event_counts, person_years)

rate_table <- data.frame(
  Group = names(event_counts),
  Events = unname(event_counts),
  `Person-years` = unname(person_years),
  `Rate per 100 person-years` = 100 * unname(event_counts / person_years),
  check.names = FALSE
)

knitr::kable(
  rate_table,
  digits = 2,
  caption = "Event counts, person-time, and rates in the two groups"
)

Event counts, person-time, and rates in the two groups

Group	Events	Person-years	Rate per 100 person-years
New treatment	18	410	4.39
Standard care	31	395	7.85

data.frame(
  `Rate ratio` = unname(rate_result$estimate),
  `Lower 95% CI` = unname(rate_result$conf.int[1]),
  `Upper 95% CI` = unname(rate_result$conf.int[2]),
  `p value` = rate_result$p.value,
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "Exact two-sample Poisson rate test")

Exact two-sample Poisson rate test

	Rate ratio	Lower 95% CI	Upper 95% CI	p value
New treatment	0.559	0.295	1.03	0.062

Verify that the direction of the rate ratio returned by poisson.test() matches the input order. A simple Poisson method assumes that events arise through a process compatible with the person-time and that there is no unaddressed overdispersion, within-patient dependence from recurrent events, or clustering. If the same patient can experience repeated events, or if adjustment for age, center, and follow-up characteristics is needed, consider Poisson or negative binomial regression, robust variance estimation, or a dedicated recurrent-event model.

6 Correlation Tests

6.1 Pearson, Spearman, and Kendall

Correlation analysis describes the degree to which two variables vary together:

• Pearson $r$: Measures linear association. It is sensitive to outliers, and its interval and test commonly rely on independent paired observations and conditions such as approximate bivariate normality.
• Spearman $\rho$: Measures monotonic association after replacing values with ranks; it is useful for ordinal data or monotonic but nonlinear relationships.
• Kendall $\tau$: Is based on concordant and discordant pairs and often has a clear interpretation when the sample is small or tied ranks are common.

plot(
  trial_data$age,
  trial_data$baseline_sbp,
  pch = 16,
  cex = 0.7,
  col = grDevices::adjustcolor(palette_test[["blue"]], alpha.f = 0.45),
  xlab = "Age (years)",
  ylab = "Baseline systolic blood pressure (mmHg)"
)
abline(
  lm(baseline_sbp ~ age, data = trial_data),
  col = palette_test[["vermillion"]],
  lwd = 2.5
)

Scatterplot of baseline age and systolic blood pressure with a fitted linear trend.

pearson_result <- cor.test(
  trial_data$age,
  trial_data$baseline_sbp,
  method = "pearson"
)
spearman_result <- cor.test(
  trial_data$age,
  trial_data$baseline_sbp,
  method = "spearman",
  exact = FALSE
)
kendall_result <- cor.test(
  trial_data$age,
  trial_data$baseline_sbp,
  method = "kendall",
  exact = FALSE
)

data.frame(
  Method = c("Pearson", "Spearman", "Kendall"),
  Estimate = c(
    unname(pearson_result$estimate),
    unname(spearman_result$estimate),
    unname(kendall_result$estimate)
  ),
  `Lower 95% CI` = c(pearson_result$conf.int[1], NA, NA),
  `Upper 95% CI` = c(pearson_result$conf.int[2], NA, NA),
  `p value` = c(
    pearson_result$p.value,
    spearman_result$p.value,
    kendall_result$p.value
  ),
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "Three commonly used correlation analyses")

Three commonly used correlation analyses

Method	Estimate	Lower 95% CI	Upper 95% CI	p value
Pearson	0.001	-0.109	0.111	0.985
Spearman	-0.006	NA	NA	0.918
Kendall	-0.003	NA	NA	0.935

Base R directly provides the usual confidence interval for Pearson correlation. If an interval for a rank correlation is needed, use a bootstrap method that respects the sampling design or dedicated software. A correlation coefficient near 0 rules out only the corresponding type of association, not every possible relationship.

Correlation is neither agreement nor causation Two instruments can have a Pearson correlation near 1 while still exhibiting a fixed bias or a bias that changes with the measurement level. Agreement between continuous measurements should be evaluated using a difference plot, Bland–Altman limits, or an appropriate ICC; agreement between categorical classifications may be assessed with kappa. Correlation can also arise from confounding, range restriction, selection mechanisms, or a shared time trend, so a p value alone cannot justify a causal interpretation.

7 Time-to-Event Outcomes: Kaplan–Meier and Log-Rank

7.1 Why Event Times Require Specialized Methods

Survival data contain whether an event occurred, the time of the event or censoring, and a risk set that changes over time. Treating censored participants as “event-free” or comparing mean times only among participants with observed events misuses the available information.

The log-rank test compares the complete survival curves of two groups and is generally efficient for proportional-hazards-type differences. Its null hypothesis can be understood as asking whether the observed event count at each event time is compatible with the count expected if the two groups had the same survival experience throughout follow-up.

survival_object <- survival::Surv(
  survival_data$time_months,
  survival_data$event
)
km_fit <- survival::survfit(
  survival_object ~ group,
  data = survival_data
)
logrank_result <- survival::survdiff(
  survival_object ~ group,
  data = survival_data
)
logrank_p <- pchisq(
  logrank_result$chisq,
  df = length(logrank_result$n) - 1,
  lower.tail = FALSE
)

data.frame(
  Group = names(logrank_result$n),
  `Sample size` = unname(logrank_result$n),
  `Observed events` = unname(logrank_result$obs),
  `Expected events under the null` = unname(logrank_result$exp),
  check.names = FALSE
) |>
  knitr::kable(digits = 2, caption = "Observed and expected event counts for the log-rank test")

Observed and expected event counts for the log-rank test

Group	Sample size.Var1	Sample size.Freq	Observed events	Expected events under the null
group=Standard care	A	130	86	63
group=New treatment	B	130	61	84

data.frame(
  `Chi-squared statistic` = unname(logrank_result$chisq),
  df = length(logrank_result$n) - 1,
  `p value` = logrank_p,
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "Log-rank test comparing the two survival curves")

Log-rank test comparing the two survival curves

Chi-squared statistic	df	p value
15.1	1	0

plot(
  km_fit,
  col = c(palette_test[["navy"]], palette_test[["vermillion"]]),
  lwd = 2.5,
  mark.time = TRUE,
  xlab = "Follow-up time (months)",
  ylab = "Event-free probability",
  conf.int = FALSE
)
legend(
  "bottomleft",
  legend = levels(survival_data$group),
  col = c(palette_test[["navy"]], palette_test[["vermillion"]]),
  lwd = 2.5,
  bty = "n"
)

Kaplan–Meier survival curves for the two groups; tick marks indicate censoring.

km_12 <- summary(km_fit, times = 12, extend = TRUE)
km_12_table <- data.frame(
  Group = sub("group=", "", km_12$strata),
  `Time (months)` = km_12$time,
  `Number at risk` = km_12$n.risk,
  `Event-free probability` = km_12$surv,
  `Lower 95% CI` = km_12$lower,
  `Upper 95% CI` = km_12$upper,
  check.names = FALSE
)

knitr::kable(
  km_12_table,
  digits = 3,
  caption = "Kaplan–Meier event-free probabilities at 12 months"
)

Kaplan–Meier event-free probabilities at 12 months

Group	Time (months)	Number at risk	Event-free probability	Lower 95% CI	Upper 95% CI
Standard care	12	40	0.423	0.342	0.522
New treatment	12	57	0.656	0.576	0.745

cox_model <- survival::coxph(
  survival::Surv(time_months, event) ~ group,
  data = survival_data
)
cox_ci <- exp(confint(cox_model))
data.frame(
  Comparison = "New treatment vs standard care",
  `Hazard ratio` = exp(coef(cox_model)),
  `Lower 95% CI` = cox_ci[1, 1],
  `Upper 95% CI` = cox_ci[1, 2],
  `p value` = summary(cox_model)$coefficients[1, "Pr(>|z|)"],
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "Hazard ratio from an unadjusted Cox model")

Hazard ratio from an unadjusted Cox model

	Comparison	Hazard ratio	Lower 95% CI	Upper 95% CI	p value
groupNew treatment	New treatment vs standard care	0.523	0.375	0.73	0

The log-rank test itself does not provide an effect size, so also report absolute Kaplan–Meier survival probabilities, the time point, the risk set, and confidence intervals. If a constant HR is reported, assess the proportional hazards assumption. When curves cross substantially, a single log-rank p value or HR may not summarize the difference well; consider an estimand aligned with the question, such as restricted mean survival time.

For more on time zero, censoring, proportional hazards diagnostics, Cox PH models, and AFT models, continue with the project’s dedicated survival-analysis module.

8 Diagnostic Testing and Agreement

8.1 Sensitivity, Specificity, and Predictive Values

Do not confuse a “statistical test” with a “diagnostic test.” Evaluating a diagnostic test requires prespecifying the target population, reference standard, threshold, temporal relationship, and approach to uncertainty.

diagnostic_table <- matrix(
  c(86, 14, 27, 173),
  nrow = 2,
  byrow = TRUE,
  dimnames = list(
    `Reference standard` = c("Disease present", "Disease absent"),
    `New test` = c("Positive", "Negative")
  )
)

tp <- diagnostic_table["Disease present", "Positive"]
fn <- diagnostic_table["Disease present", "Negative"]
fp <- diagnostic_table["Disease absent", "Positive"]
tn <- diagnostic_table["Disease absent", "Negative"]

sensitivity_ci <- binom.test(tp, tp + fn)$conf.int
specificity_ci <- binom.test(tn, tn + fp)$conf.int
ppv_ci <- binom.test(tp, tp + fp)$conf.int
npv_ci <- binom.test(tn, tn + fn)$conf.int

knitr::kable(
  addmargins(diagnostic_table),
  caption = "Cross-tabulation of the new test and reference standard"
)

Cross-tabulation of the new test and reference standard

	Positive	Negative	Sum
Disease present	86	14	100
Disease absent	27	173	200
Sum	113	187	300

diagnostic_metrics <- data.frame(
  Metric = c(
    "Sensitivity",
    "Specificity",
    "Positive predictive value",
    "Negative predictive value"
  ),
  Estimate = c(
    tp / (tp + fn),
    tn / (tn + fp),
    tp / (tp + fp),
    tn / (tn + fn)
  ),
  `Lower 95% CI` = c(
    sensitivity_ci[1], specificity_ci[1], ppv_ci[1], npv_ci[1]
  ),
  `Upper 95% CI` = c(
    sensitivity_ci[2], specificity_ci[2], ppv_ci[2], npv_ci[2]
  ),
  check.names = FALSE
)

knitr::kable(
  diagnostic_metrics,
  digits = 3,
  caption = "Diagnostic performance with exact binomial 95% confidence intervals"
)

Diagnostic performance with exact binomial 95% confidence intervals

Metric	Estimate	Lower 95% CI	Upper 95% CI
Sensitivity	0.860	0.776	0.921
Specificity	0.865	0.810	0.909
Positive predictive value	0.761	0.672	0.836
Negative predictive value	0.925	0.878	0.958

Sensitivity and specificity condition on the reference standard. Positive and negative predictive values instead condition on the test result and change with disease prevalence in the target population. Case-control-type sampling can estimate sensitivity and specificity, but the artificially constructed case fraction in the sample usually cannot be used directly to estimate PPV and NPV in a clinical setting.

If the same patient receives two tests, comparisons of positivity, sensitivity, or specificity have a paired design. For example, to compare sensitivity, apply a McNemar-type analysis to the paired results of the two tests among patients with disease according to the reference standard, rather than treating the tests as two independent samples.

9 Pre-Test Data Checks and Method Assumptions

9.1 Check Denominators, Coding, Missingness, and Plots First

At a minimum, a data audit should include:

• the number enrolled, number actually analyzed, and number missing in each group;
• whether any patient was duplicated unintentionally and whether paired identifiers are unique;
• whether units, factor levels, event coding, and reference groups are correct;
• impossible values, data-entry errors, and measurement floors or ceilings;
• raw distributions, outliers, scatterplots, and time trajectories by group;
• whether complete-case selection changed the target population represented by the analysis sample.

data.frame(
  Variable = names(trial_data),
  Class = vapply(trial_data, function(x) class(x)[1], character(1)),
  Missing = vapply(trial_data, function(x) sum(is.na(x)), integer(1)),
  `Unique values` = vapply(
    trial_data,
    function(x) length(unique(x)),
    integer(1)
  ),
  check.names = FALSE
) |>
  knitr::kable(caption = "Basic integrity checks for the simulated trial data")

Basic integrity checks for the simulated trial data

	Variable	Class	Missing	Unique values
participant_id	participant_id	character	0	320
treatment	treatment	factor	0	2
age	age	numeric	0	56
sex	sex	factor	0	2
baseline_sbp	baseline_sbp	numeric	0	243
followup_sbp	followup_sbp	numeric	0	251
sbp_reduction	sbp_reduction	numeric	0	212
crp_week12	crp_week12	numeric	0	261
response	response	factor	0	2
adverse_event	adverse_event	factor	0	2

stopifnot(!anyDuplicated(trial_data$participant_id))

An outlier must not be removed simply because deleting it produces statistical significance. First determine whether it is a data-entry error, a measurement error, a rare but protocol-consistent true value, or an observation outside the target population. Rules for retaining or excluding observations should follow data quality considerations and the study protocol, and sensitivity analyses should show their impact.

9.2 Normality: Plots and the Scientific Question Matter More Than One p Value

The relevant normality assumption for a t test or ANOVA concerns model errors, within-group distributions, or paired differences, not the outcome after all groups have been pooled. Judge it using Q–Q plots, sample size, degree of skewness, outliers, and the target estimand.

anova_residuals <- residuals(traditional_anova)
qqnorm(
  anova_residuals,
  pch = 16,
  col = grDevices::adjustcolor(palette_test[["blue"]], alpha.f = 0.55),
  main = "Q–Q Plot of ANOVA Residuals"
)
qqline(anova_residuals, col = palette_test[["vermillion"]], lwd = 2.5)

Normal Q–Q plot of residuals from the three-group ANOVA.

shapiro_result <- shapiro.test(anova_residuals)
fligner_result <- fligner.test(
  hemoglobin_change ~ dose,
  data = dose_data
)

data.frame(
  Check = c(
    "Shapiro–Wilk test of residual normality",
    "Fligner–Killeen test of equal variances"
  ),
  Statistic = c(
    unname(shapiro_result$statistic),
    unname(fligner_result$statistic)
  ),
  `p value` = c(shapiro_result$p.value, fligner_result$p.value),
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "Two common tests used to check assumptions")

Two common tests used to check assumptions

Check	Statistic	p value
Shapiro–Wilk test of residual normality	0.99	0.235
Fligner–Killeen test of equal variances	9.05	0.011

Do not use Shapiro–Wilk as an automatic switch With a small sample, the test may fail to detect an important departure; with a large sample, it may return a very small p value for a minor departure with little practical impact. Do not automatically choose a method using a two-stage rule such as “test normality first; use a t test if nonsignificant and Wilcoxon otherwise,” because the two methods do not have the same estimand or assumptions.

Likewise, there is no need to use a preliminary equal-variance test to decide whether a Welch t test is permissible; Welch’s method is generally a reasonable default for comparing two independent means. Bartlett’s test is sensitive to nonnormality, whereas Fligner–Killeen is more robust, but no variance test replaces residual plots and subject-matter judgment.

9.3 Independence Usually Cannot Be Tested from the Same Data

Independence is determined primarily by sampling, randomization, follow-up, and the data hierarchy. If the same patient contributes two eyes, multiple lesions, or multiple hospitalizations, a histogram cannot reveal that the wrong unit of analysis was used. Identify the hierarchy from the data dictionary and study workflow, and use a method that represents the correlation structure.

10 Multiple Comparisons and Selective Reporting

10.1 Why Repeated Testing Accumulates False-Positive Risk

When many independent tests are performed within the same family at $\alpha=0.05$, the probability of at least one false positive increases. Multiplicity arises not only from pairwise comparisons but also from multiple outcomes, time points, subgroups, thresholds, model versions, and repeated examination of the data.

raw_p_values <- c(
  "Primary outcome" = 0.012,
  "Secondary outcome A" = 0.041,
  "Secondary outcome B" = 0.007,
  "Secondary outcome C" = 0.18,
  "Exploratory outcome" = 0.049
)

p_adjust_table <- data.frame(
  Test = names(raw_p_values),
  `Raw p value` = unname(raw_p_values),
  Bonferroni = p.adjust(raw_p_values, method = "bonferroni"),
  Holm = p.adjust(raw_p_values, method = "holm"),
  `BH FDR` = p.adjust(raw_p_values, method = "BH"),
  check.names = FALSE
)

knitr::kable(
  p_adjust_table,
  digits = 3,
  caption = "Examples of three common multiplicity adjustments"
)

Examples of three common multiplicity adjustments

	Test	Raw p value	Bonferroni	Holm	BH FDR
Primary outcome	Primary outcome	0.012	0.060	0.048	0.030
Secondary outcome A	Secondary outcome A	0.041	0.205	0.123	0.061
Secondary outcome B	Secondary outcome B	0.007	0.035	0.035	0.030
Secondary outcome C	Secondary outcome C	0.180	0.900	0.180	0.180
Exploratory outcome	Exploratory outcome	0.049	0.245	0.123	0.061

• Holm: Controls the family-wise error rate, is usually more powerful than simple Bonferroni adjustment, and is a practical general choice for confirmatory adjustment.
• Bonferroni: Is simple and transparent but may be conservative when there are many or highly correlated tests.
• Tukey: Is suitable for all pairwise mean comparisons among groups in an ANOVA.
• Dunnett-type comparisons: Are better aligned with the objective when several groups are compared only with one common control.
• Benjamini–Hochberg (BH): Controls the false discovery rate and is commonly used for exploratory or high-dimensional analyses; it does not control the probability of at least one false positive.

The strongest approach is not to search after the fact for the least stringent adjustment, but to define the primary outcome, primary time point, primary contrast, and test family in the protocol. When p values are adjusted, confidence intervals should be compatible with the same multiplicity strategy.

11 Superiority, Equivalence, and Noninferiority

11.1 Nonsignificant Does Not Mean the Same

Failure to reject “zero difference” in a conventional two-sided superiority test means only that the data provide insufficient evidence of a nonzero difference. It does not prove that two treatments are equivalent or that no clinically important difference exists.

Equivalence requires prespecified lower and upper equivalence bounds and evidence that the entire effect interval lies inside those bounds. At significance level $\alpha$, two one-sided tests (TOST) are equivalent to checking whether the corresponding $100(1-2\alpha)\%$ confidence interval lies entirely between the equivalence bounds.

mean_difference <- mean(new_values) - mean(standard_values)
standard_error <- sqrt(var(new_values) / n1 + var(standard_values) / n0)
welch_df <- (var(new_values) / n1 + var(standard_values) / n0)^2 /
  ((var(new_values) / n1)^2 / (n1 - 1) +
     (var(standard_values) / n0)^2 / (n0 - 1))

lower_equivalence_bound <- -6
upper_equivalence_bound <- 6
t_lower <- (mean_difference - lower_equivalence_bound) / standard_error
t_upper <- (mean_difference - upper_equivalence_bound) / standard_error
p_lower <- pt(t_lower, df = welch_df, lower.tail = FALSE)
p_upper <- pt(t_upper, df = welch_df, lower.tail = TRUE)
tost_p <- max(p_lower, p_upper)

ci90 <- mean_difference +
  c(-1, 1) * qt(0.95, df = welch_df) * standard_error

data.frame(
  `Mean difference` = mean_difference,
  `Lower 90% CI` = ci90[1],
  `Upper 90% CI` = ci90[2],
  `Lower equivalence bound` = lower_equivalence_bound,
  `Upper equivalence bound` = upper_equivalence_bound,
  `TOST p value` = tost_p,
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "Illustrative TOST equivalence example")

Illustrative TOST equivalence example

Mean difference	Lower 90% CI	Upper 90% CI	Lower equivalence bound	Upper equivalence bound	TOST p value
4.55	2.97	6.13	-6	6	0.066

Equivalence bounds must come from clinical judgment, prior evidence, and the protocol; they must not be chosen after viewing the data merely to make the interval “fit.” An effect can be statistically different from 0 and also lie within a prespecified broad equivalence interval. “There is a difference” and “the difference is small enough to be considered equivalent” answer different questions.

Noninferiority instead concerns only one prespecified direction—for example, demonstrating that the efficacy of a new treatment is not worse than that of the control by more than a margin $\Delta$. The direction, margin, analysis populations, and roles of intention-to-treat and per-protocol analyses should all be defined in advance. Switching to a one-sided test or choosing a new margin after seeing the results undermines error-rate control.

12 Tests, Regression Models, and Complex Designs

Many basic tests can be expressed as regression models:

Basic test	Corresponding model perspective	Extensions
Independent two-sample t test	Linear model with only a binary group indicator	Add baseline values, covariates, interactions, and nonlinear terms
One-way ANOVA	Linear model with group as a factor predictor	Prespecified contrasts, unbalanced designs, and covariate adjustment
Two proportions / 2×2 table	Binary regression model	Adjust for confounding; estimate OR, RR, or RD; assess interactions
Poisson rate test	Poisson model with person-time as an offset	Multivariable rate ratios and methods for overdispersion or clustering
Log-rank test	Unadjusted comparison of survival curves	Cox, AFT, RMST, and covariate adjustment

test_as_model <- lm(sbp_reduction ~ treatment, data = trial_data)
model_ci <- confint(test_as_model)["treatmentNew treatment", ]

data.frame(
  Method = "Linear model: new treatment vs standard care",
  `Mean difference` = coef(test_as_model)["treatmentNew treatment"],
  `Lower 95% CI` = model_ci[1],
  `Upper 95% CI` = model_ci[2],
  `p value` = summary(test_as_model)$coefficients[
    "treatmentNew treatment",
    "Pr(>|t|)"
  ],
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "Linear-model expression of a two-group mean comparison")

Linear-model expression of a two-group mean comparison

	Method	Mean difference	Lower 95% CI	Upper 95% CI	p value
treatmentNew treatment	Linear model: new treatment vs standard care	4.55	2.66	6.43	0

The conventional standard error from ordinary lm() corresponds to an equal-variance model, so its numerical results need not be identical to those from a Welch t test. Appropriate variance estimators and model structures can relax this assumption.

When a study contains multiple factors, center effects, clustering, nonlinearity, interactions, censoring, or incomplete follow-up, a single basic test is often insufficient. Also remember that a small p value in one subgroup and a large p value in another do not prove that effects differ between subgroups; estimate the interaction directly and report its interval.

13 Sample Size, Power, and Clinical Importance

13.1 Prospective Planning

Sample size should be planned in advance around the primary target estimand, minimum clinically important difference, anticipated variability, significance level, target power, allocation ratio, loss to follow-up, and design effect.

continuous_power <- power.t.test(
  delta = 5,
  sd = 9,
  sig.level = 0.05,
  power = 0.80,
  type = "two.sample",
  alternative = "two.sided"
)

binary_power <- power.prop.test(
  p1 = 0.55,
  p2 = 0.70,
  sig.level = 0.05,
  power = 0.80,
  alternative = "two.sided"
)

data.frame(
  Scenario = c(
    "Difference of 5 between two independent means, SD 9",
    "Two independent proportions, 0.55 vs 0.70"
  ),
  `Complete cases required per group` = ceiling(
    c(continuous_power$n, binary_power$n)
  ),
  check.names = FALSE
) |>
  knitr::kable(caption = "Sample-size examples under simplified assumptions")

Sample-size examples under simplified assumptions

Scenario	Complete cases required per group
Difference of 5 between two independent means, SD 9	52
Two independent proportions, 0.55 vs 0.70	163

These functions apply to simplified designs with two independent groups. Paired, clustered, repeated-measures, and survival designs—as well as unequal allocation, multiple primary outcomes, or complex models—require specialized formulas or simulation. The calculated number of complete observations should also be inflated for a plausible loss-to-follow-up rate.

Using the observed effect after results are available to calculate “post hoc power” is not recommended as an explanation for a nonsignificant result, because it is usually little more than a restatement of the p value. Instead, examine the effect estimate and 95% CI: does the interval exclude clinically important benefit, harm, neither, or both?

14 Interpreting p Values and Confidence Intervals Correctly

14.1 What a p Value Can and Cannot Say

A p value is the probability, under the null hypothesis and the statistical model assumptions, of obtaining the observed result or a result at least as incompatible with the null hypothesis. It is not:

• the probability that the null hypothesis is true;
• the probability that the result was “caused by chance”;
• a measure of effect size or clinical importance;
• proof that the study is free of bias, confounding, or measurement problems;
• the probability that another study will necessarily reproduce the result.

With $\alpha=0.05$ fixed, a Type I error occurs when a true null hypothesis is rejected; a Type II error occurs when the null is not rejected despite an effect as defined by the study. Power is the probability of correctly rejecting the null hypothesis for a specified true effect and set of design conditions.

A confidence interval places effect size and uncertainty on the same scale. The frequentist interpretation of a 95% confidence interval is that, under repeated sampling with the same design and procedure, about 95% of intervals constructed this way would cover the true parameter in the long run. It does not mean that “this particular observed interval has a 95% probability of containing the true value.”

14.2 One-Sided and Two-Sided Tests

Medical research generally uses two-sided tests because unexpected harm or an effect in the opposite direction can be equally important. A one-sided test may be reasonable only when an effect in the opposite direction truly requires no distinction scientifically or for decision-making, the direction was documented in the protocol before the data were examined, and regulatory or field-specific standards permit it. Changing from two-sided to one-sided after seeing the result selectively reduces the p value.

A better reading order: First inspect the denominator and descriptive statistics in each group; next examine the effect estimate in the prespecified direction and its 95% CI; then consider the p value and multiplicity adjustment; finally assess the design, data quality, clinical importance, and generalizability.

15 How to Report Test Results

15.1 Minimum Reporting Checklist

An auditable medical-statistics result should state at least:

• the study design, unit of analysis, and grouping or pairing structure;
• the actual number analyzed, missing count, and key denominator in each group;
• descriptive statistics matched to the outcome scale;
• the contrast direction, effect size, units, and 95% CI;
• the full name of the test or model and whether it was one- or two-sided;
• the exact p value (write p < 0.001 for very small values, not p = 0);
• multiplicity handling, assumption checks, and any sensitivity analyses;
• the scope of applicability and limitations from the design, bias, or model.

15.2 Reusable Reporting Template

In [population and design], [outcome] was [descriptive statistics with denominators] in [Group A] and [Group B], respectively. The [effect measure] in the prespecified direction was [estimate and units] (95% CI [lower, upper]); [full test name, one- or two-sided, and adjustment method] gave $p=$[value]. The analysis was based on [actual sample size] complete/available records and is subject to [specific assumptions, missingness, bias, or generalizability limitations].

15.3 Continuous-Outcome Example from This Tutorial

report_difference <- mean(new_values) - mean(standard_values)
report_ci <- -rev(welch_result$conf.int)
report_p <- format_p(welch_result$p.value)

cat(
  paste0(
    "The new treatment and standard care groups included ", n1, " and ", n0,
    " participants, respectively. Their mean reductions in systolic blood ",
    "pressure were ",
    sprintf("%.1f", mean(new_values)), " (SD ",
    sprintf("%.1f", sd(new_values)), ") and ",
    sprintf("%.1f", mean(standard_values)), " (SD ",
    sprintf("%.1f", sd(standard_values)), ") mmHg, respectively. ",
    "The mean difference for new treatment minus standard care was ",
    sprintf("%.1f", report_difference), " mmHg (95% CI ",
    sprintf("%.1f", report_ci[1]), " to ",
    sprintf("%.1f", report_ci[2]), "; two-sided Welch t test p ",
    ifelse(
      welch_result$p.value < 0.001,
      "< 0.001",
      paste0("= ", report_p)
    ),
    "). This estimate comes from simulated complete data; clinical ",
    "interpretation must also consider the prespecified minimum clinically ",
    "important difference."
  )
)

## The new treatment and standard care groups included 160 and 160 participants, respectively. Their mean reductions in systolic blood pressure were 11.6 (SD 8.9) and 7.1 (SD 8.2) mmHg, respectively. The mean difference for new treatment minus standard care was 4.5 mmHg (95% CI 2.7 to 6.4; two-sided Welch t test p < 0.001). This estimate comes from simulated complete data; clinical interpretation must also consider the prespecified minimum clinically important difference.

16 Quick Guide to Common Errors

Incorrect practice	Why it is a problem	Better approach
Write “the groups are the same” whenever p > 0.05	Failure to reject the null does not prove equality	Report the effect and CI; use prespecified bounds for an equivalence question
Report only that a difference is “statistically significant”	This omits direction, magnitude, units, and precision	Report the effect on the original scale, its CI, and clinical importance
Run Shapiro first and then choose a method mechanically	A two-stage rule ignores the estimand and robustness	Define the question and design first, then use plots and sensitivity analyses
Call Mann–Whitney a universal test of medians	In general, it tests rank distributions	State the additional shape assumption or describe the target accurately
Treat paired data as independent samples	This ignores within-patient dependence	Use a paired test or repeated-measures model
Run three unadjusted tests at three time points	This increases the family-wise false-positive rate	Fit an overall model first, then use prespecified adjusted contrasts
Delete an outlier because it changes significance	This creates outcome-driven data handling	Investigate its source, follow the protocol, and perform sensitivity analyses
Treat an OR as an RR	They can differ substantially when the outcome is common	Report RD, RR, or OR as appropriate to the design and name the scale
Test every baseline variable for significance in an RCT	Postrandomization differences arise by chance by design	Describe baseline variables; adjust as prespecified for prognostic value
Claim interaction because one subgroup is significant and another is not	A difference between two p values is not a p value for the difference	Estimate the interaction directly and report its CI
Switch to a one-sided test after seeing the result	This invalidates the prespecified error rate	Decide the direction and rationale in the protocol
Hide the sample size used in different analyses	Complete-case analysis may change the sample and target population	Report the actual denominator and missingness for every analysis
Keep changing tests, thresholds, or outcomes until something is significant	This produces selective reporting and unstable conclusions	Preregister the primary analysis and disclose the full exploratory process
Treat correlation as agreement or causation	The three concepts answer different questions	Use an agreement method and discuss causation based on the design

17 Exercises and Answers

17.1 Exercise 1: Two Independent Means

A randomized trial compares 8-week change in LDL-C between two groups, with each participant belonging to only one group. The group variances do not appear identical, and the question concerns the difference in mean change. What is the most direct basic test, and what should be reported?

View answer

Use Welch’s independent two-sample t test. Report the sample size, mean, and SD in each group, together with the difference in mean change in the prespecified direction, its 95% CI, and a two-sided p value. A significant Welch p value is not a reason to omit the difference in clinically meaningful units.

17.2 Exercise 2: Left and Right Eyes from the Same Patient

Investigators compare intraocular pressure in the left and right eyes of the same patients using an independent-samples t test. What is the fundamental error?

View answer

The left and right eyes belong to the same patient, so the observations are not independent. If every patient contributes exactly one pair and the question concerns a left–right difference, analyze the within-patient difference. More complex designs involving eye, treatment, repeated visits, or enrollment of both eyes require a model that represents eyes nested within patients.

17.3 Exercise 3: Overall Test for Three Groups

The p value from a three-group ANOVA is 0.003. Can you write that “all three groups differ from one another”?

View answer

No. The omnibus test shows only that at least one group mean is incompatible with the null hypothesis that all means are equal. Estimate prespecified contrasts or use a method matched to the comparison family, such as Tukey, Dunnett, or Holm adjustment, and report simultaneous intervals or adjusted inferences.

17.4 Exercise 4: Sparse Adverse Events

Two groups contain 25 participants each. One has no severe adverse events and the other has three. Why is it insufficient to report only a Fisher exact-test p value?

View answer

The small sample and zero event count make the effect estimate very imprecise. Also report the denominator and risk in each group, a risk difference or suitable relative effect with an interval, and the estimation issue caused by the zero cell. A nonsignificant result does not exclude clinically important harm.

17.5 Exercise 5: Positivity Before and After an Intervention

The same patients receive the same binary test before and after an intervention. Which test uses the correct paired structure?

View answer

The McNemar test. If there are very few discordant pairs, an exact binomial test can be applied to the two directions of discordance. Also report the marginal proportions before and after and their risk difference, not only the p value.

17.6 Exercise 6: Very High Correlation

A new instrument and a standard instrument have Pearson $r=0.97$. Is this sufficient to establish that the instruments are interchangeable?

View answer

No. High correlation can occur in the presence of a fixed or proportional bias. Examine how paired differences vary across the measurement range, report Bland–Altman bias and limits of agreement, and consider the ICC, repeatability, and acceptable error range for the intended use.

17.7 Exercise 7: Crossing Survival Curves

Two Kaplan–Meier curves cross at 6 months, after which the direction of the difference reverses. Why might a single log-rank p value be inadequate?

View answer

Crossing suggests that the effect changes over time, which can limit both the log-rank test’s advantage under proportional-hazards-type differences and the interpretation of a single constant HR. Show the curves, risk table, and absolute survival probabilities at key times; assess time-varying effects; and consider RMST or piecewise effects that match the study question.

17.8 Exercise 8: A Nonsignificant Result

The estimated treatment effect is a risk difference of $-2\%$, with a 95% CI from $-18\%$ to $14\%$ and $p=0.81$. Can you conclude that “there is no treatment effect”?

View answer

No. The interval contains potentially meaningful benefit and harm, which indicates imprecision. You may state that the data do not show a clear difference, but you cannot establish no effect. Discuss sample size, event count, clinically important bounds, and possible study biases.

18 Quick Reference

18.1 Common R Code

Analysis objective	Base R or recommended-package code pattern	Report first
One-sample mean	t.test(x, mu = value)	Difference between the mean and reference value, 95% CI
Two independent means	t.test(y ~ group)	Mean difference, 95% CI
Paired means	t.test(after, before, paired = TRUE)	Mean within-person difference, 95% CI
Two independent rank distributions	wilcox.test(y ~ group, conf.int = TRUE)	Distributions, location shift, and CI
Paired rank distributions	wilcox.test(after, before, paired = TRUE)	Difference distribution and location shift
Multiple group means	oneway.test(y ~ group, var.equal = FALSE)	Prespecified between-group mean differences
Conventional ANOVA	aov(y ~ group)	Omnibus test, specific contrasts, and effect size
Multiple group rank distributions	kruskal.test(y ~ group)	Omnibus rank difference and adjusted comparisons
Complete-block repeated measures	friedman.test(y ~ time \| id)	Time-point differences and adjusted comparisons
One proportion	binom.test(x, n)	Proportion and exact CI
Two independent proportions	prop.test(x, n); chisq.test(tab)	RD, RR, or OR with CI
Sparse 2×2 table	fisher.test(tab)	Denominators, absolute risks, OR, and CI
Paired binary outcome	mcnemar.test(tab)	Marginal risk difference and discordant pairs
Trend in ordered proportions	prop.trend.test(x, n)	Proportions at each level and trend test
Stratified 2×2 tables	mantelhaen.test(array)	Stratum-specific results and common OR with CI
Two incidence rates	poisson.test(events, person_time)	Rates and rate ratio with CI
Pearson correlation	cor.test(x, y, method = "pearson")	$r$ and CI
Rank correlation	cor.test(x, y, method = "spearman")	$\rho$ and an appropriate interval
Survival-curve comparison	survival::survdiff(Surv(time, event) ~ group)	KM probabilities, risk sets, effect, and CI
Multiple-p-value adjustment	p.adjust(p, method = "holm")	Adjustment method, adjusted p values, and simultaneous CIs

18.2 Null Values at a Glance

Effect scale	Null value	Typical interpretation
Mean difference, risk difference, rate difference, correlation coefficient	0	No difference or association on that scale
Risk ratio, rate ratio, odds ratio, hazard ratio	1	The groups are equal on the relative scale
Sensitivity, specificity, predictive values	No universal null value	Compare with the intended use, threshold, and benchmark

Final Analysis Checklist

Before releasing results, confirm that:

• the research question, population, time window, and target estimand were specified in advance;
• the unit of analysis and the independent, paired, repeated, clustered, or censored structure were identified correctly;
• outcome coding, contrast direction, units, reference group, and denominators were verified;
• data quality, missingness, outliers, and group distributions were examined graphically;
• method assumptions match the study design, and the method was not chosen mechanically from a single normality p value;
• every primary conclusion includes an effect size on the original scale and a 95% CI;
• the one- or two-sided direction, exact method, and adjustment for each p value are stated;
• multiplicity from multiple outcomes, time points, subgroups, and model choices was handled or disclosed transparently;
• statistical significance was not misrepresented as clinical importance, and nonsignificance was not misrepresented as equality;
• OR, RR, HR, rate ratio, correlation, and causation were not conflated;
• the actual analysis sample size, missing records, and sensitivity analyses were reported;
• code, random seeds, software versions, and the data-processing trail are sufficient for reproducibility.

sessionInfo()

## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Tahoe 26.5.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
## 
## time zone: America/Edmonton
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.39   R6_2.6.1        fastmap_1.2.0   Matrix_1.7-5   
##  [5] xfun_0.60       lattice_0.22-9  splines_4.6.1   cachem_1.1.0   
##  [9] knitr_1.51      htmltools_0.5.9 rmarkdown_2.31  lifecycle_1.0.5
## [13] cli_3.6.6       grid_4.6.1      sass_0.4.10     jquerylib_0.1.4
## [17] compiler_4.6.1  tools_4.6.1     evaluate_1.0.5  bslib_0.12.0
## [21] survival_3.8-6  yaml_2.3.12     rlang_1.3.0     jsonlite_2.0.0

All data and results in this tutorial are simulated and are provided solely for teaching medical statistics. They do not replace statistical consultation during protocol development, clinical judgment, ethics review, or regulatory requirements. Real studies should select analysis methods according to their design, sampling process, data quality, and target estimand.