A Complete Guide to Statistical Tests and Methods for Clinical Researchers

install.packages(c(
  "tidyverse", "ggpubr", "rstatix", "broom", "coin",
  "survival", "survminer", "lme4", "lmerTest",
  "pROC", "irr", "meta", "vcd", "ggfortify",
  "factoextra", "psych", "MatchIt", "WeightIt"
))

pip install pandas numpy scipy matplotlib seaborn statsmodels \
            lifelines scikit-learn pingouin scikit-posthocs

n per group = 2 × (z_α/2 + z_β)² × σ² / δ²

library(tidyverse)
 
sd_val <- 20; alpha <- 0.05
 
expand_grid(delta = c(5, 8, 10, 15), n = seq(10, 150, by = 5)) |>
  mutate(
    power = pmap_dbl(list(delta, n), function(d, n) {
      power.t.test(n = n, delta = d, sd = sd_val, sig.level = alpha,
                   type = "two.sample")$power
    }),
    label = paste0("MCID = ", delta, " mmHg")
  ) |>
  ggplot(aes(x = n, y = power, colour = label)) +
  geom_line(linewidth = 1.1) +
  geom_hline(yintercept = 0.80, linetype = "dashed", colour = "grey40") +
  annotate("text", x = 12, y = 0.82, label = "80% power", hjust = 0, size = 3.5) +
  scale_y_continuous(labels = scales::percent_format(), limits = c(0, 1)) +
  scale_colour_brewer(palette = "Set1") +
  labs(title = "Power Curves: Two-Sample t-Test (SD = 20 mmHg)",
       x = "Sample size per group", y = "Statistical power", colour = NULL) +
  theme_bw(base_size = 13) + theme(legend.position = "bottom")

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
 
sd_val = 20; alpha = 0.05
n_vals = np.arange(10, 151, 5)
mcids  = [5, 8, 10, 15]
colours = ['#e41a1c','#377eb8','#4daf4a','#984ea3']
 
fig, ax = plt.subplots(figsize=(9, 5))
for mcid, col in zip(mcids, colours):
    se     = sd_val * np.sqrt(2 / n_vals)
    ncp    = mcid / se
    z_crit = stats.norm.ppf(1 - alpha / 2)
    power  = 1 - stats.norm.cdf(z_crit - ncp) + stats.norm.cdf(-z_crit - ncp)
    ax.plot(n_vals, power, label=f"MCID = {mcid} mmHg", color=col, linewidth=2)
 
ax.axhline(0.80, linestyle='--', color='grey', linewidth=1.2)
ax.text(12, 0.82, "80% power", fontsize=10, color='grey')
ax.set_xlabel("Sample size per group"); ax.set_ylabel("Statistical power")
ax.set_title("Power Curves: Two-Sample t-Test (SD = 20 mmHg)")
ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda y, _: f'{y:.0%}'))
ax.legend(loc='lower right'); ax.grid(alpha=0.3)
plt.tight_layout(); plt.savefig("power_curves.png", dpi=150, bbox_inches='tight')

library(tidyverse)
 
set.seed(42)
n      <- 120
normal <- rnorm(n, 62.4, 11.2)
skewed <- rgamma(n, shape = 2, rate = 0.6)
 
par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
 
hist(normal, breaks = 15, col = "#4575b4", border = "white",
     main = "Age — approx. normal", xlab = "Age (years)", freq = FALSE)
curve(dnorm(x, mean(normal), sd(normal)), add = TRUE, col = "firebrick", lwd = 2)
 
qqnorm(normal, main = "Q-Q: Age", pch = 16, col = "#4575b4")
qqline(normal, col = "firebrick", lwd = 2)
 
hist(skewed, breaks = 15, col = "#d73027", border = "white",
     main = "LDL-C — right-skewed", xlab = "LDL-C (mmol/L)", freq = FALSE)
lines(density(skewed), col = "navy", lwd = 2)
 
qqnorm(skewed, main = "Q-Q: LDL-C", pch = 16, col = "#d73027")
qqline(skewed, col = "navy", lwd = 2)
 
cat("Shapiro-Wilk (age):  p =", round(shapiro.test(normal)$p.value, 4), "\n")
cat("Shapiro-Wilk (LDL):  p =", round(shapiro.test(skewed)$p.value, 4), "\n")

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from scipy.stats import gaussian_kde
 
rng    = np.random.default_rng(42)
n      = 120
normal = rng.normal(62.4, 11.2, n)
skewed = rng.gamma(shape=2, scale=1/0.6, size=n)
 
fig, axes = plt.subplots(2, 2, figsize=(10, 8))
 
# Panel 1: Histogram normal
axes[0,0].hist(normal, bins=15, color='#4575b4', edgecolor='white', density=True)
x_r = np.linspace(normal.min(), normal.max(), 200)
axes[0,0].plot(x_r, stats.norm.pdf(x_r, normal.mean(), normal.std()), 'r-', lw=2)
axes[0,0].set_title("Age — approx. normal"); axes[0,0].set_xlabel("Age (years)")
 
# Panel 2: Q-Q normal
(osm, osr), (slope, intercept, _) = stats.probplot(normal)
axes[0,1].scatter(osm, osr, color='#4575b4', s=20)
axes[0,1].plot(osm, slope*np.array(osm)+intercept, 'r-', lw=2)
axes[0,1].set_title("Q-Q Plot: Age"); axes[0,1].set_xlabel("Theoretical quantiles")
 
# Panel 3: Histogram skewed
axes[1,0].hist(skewed, bins=15, color='#d73027', edgecolor='white', density=True)
kde_x = np.linspace(skewed.min(), skewed.max(), 200)
axes[1,0].plot(kde_x, gaussian_kde(skewed)(kde_x), 'navy', lw=2)
axes[1,0].set_title("LDL-C — right-skewed"); axes[1,0].set_xlabel("LDL-C (mmol/L)")
 
# Panel 4: Q-Q skewed
(osm2, osr2), (s2, i2, _) = stats.probplot(skewed)
axes[1,1].scatter(osm2, osr2, color='#d73027', s=20)
axes[1,1].plot(osm2, s2*np.array(osm2)+i2, 'navy', lw=2)
axes[1,1].set_title("Q-Q Plot: LDL-C"); axes[1,1].set_xlabel("Theoretical quantiles")
 
for ax in axes.flat:
    ax.grid(alpha=0.3)
plt.suptitle("Normality Assessment: 4-Panel Diagnostic", fontsize=13)
plt.tight_layout(); plt.savefig("normality_4panel.png", dpi=150, bbox_inches='tight')
plt.show()
 
for label, arr in [("age", normal), ("LDL", skewed)]:
    w, p = stats.shapiro(arr)
    print(f"Shapiro-Wilk ({label}):  W={w:.3f}, p={p:.4f}")

t = (x̄ − μ₀) / (s / √n)     df = n − 1

library(tidyverse)
 
set.seed(42)
inr <- rnorm(25, mean = 2.8, sd = 0.6)
mu0 <- 2.5
 
t_res <- t.test(inr, mu = mu0)
print(t_res)
 
tibble(inr = inr) |>
  ggplot(aes(x = 1, y = inr)) +
  geom_jitter(width = 0.08, size = 2.5, colour = "#377eb8", alpha = 0.7) +
  stat_summary(fun.data = mean_cl_normal, geom = "pointrange",
               colour = "firebrick", size = 0.8, linewidth = 1.2) +
  geom_hline(yintercept = mu0, linetype = "dashed", colour = "grey30", linewidth = 1) +
  annotate("text", x = 1.35, y = mu0 + 0.05,
           label = paste0("Reference = ", mu0), hjust = 0, size = 3.5) +
  annotate("text", x = 1.35, y = mean(inr) + 0.05,
           label = sprintf("Mean = %.2f\n95%% CI: %.2f–%.2f\np = %.3f",
                           mean(inr), t_res$conf.int[1], t_res$conf.int[2], t_res$p.value),
           hjust = 0, size = 3.2, colour = "firebrick") +
  scale_x_continuous(limits = c(0.7, 1.9)) +
  labs(title    = "One-Sample t-Test: INR vs Target of 2.5",
       subtitle = "Dots = patients; red = mean ± 95% CI; dashed = reference target",
       y = "INR", x = NULL) +
  theme_bw(base_size = 13) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())

import numpy as np, matplotlib.pyplot as plt
from scipy import stats
 
rng = np.random.default_rng(42)
inr = rng.normal(2.8, 0.6, 25); mu0 = 2.5
t_stat, p_val = stats.ttest_1samp(inr, mu0)
ci = stats.t.interval(0.95, df=len(inr)-1, loc=np.mean(inr), scale=stats.sem(inr))
print(f"t={t_stat:.3f}, p={p_val:.4f}, 95% CI: {ci[0]:.2f}–{ci[1]:.2f}")
 
fig, ax = plt.subplots(figsize=(5, 6))
jitter = rng.uniform(-0.05, 0.05, len(inr))
ax.scatter(jitter, inr, color='#377eb8', alpha=0.65, s=50, zorder=3)
ax.errorbar(0, np.mean(inr), yerr=[[np.mean(inr)-ci[0]], [ci[1]-np.mean(inr)]],
            fmt='o', color='firebrick', markersize=8, linewidth=2.5, capsize=6, zorder=4)
ax.axhline(mu0, linestyle='--', color='grey', linewidth=1.5)
ax.set_xlim(-0.3, 0.55); ax.set_xticks([])
ax.set_ylabel("INR", fontsize=12)
ax.set_title(f"One-Sample t-Test: INR vs Target 2.5\nt={t_stat:.2f}, p={p_val:.3f}", fontsize=12)
ax.grid(axis='y', alpha=0.3); plt.tight_layout()
plt.savefig("onesample_t.png", dpi=150, bbox_inches='tight'); plt.show()

library(tidyverse)
 
pain <- c(7,8,6,9,5,7,8,7,6,9,7,8,5,7,9,6,8,7)
wilcox.test(pain, mu = 5, exact = FALSE)
 
tibble(pain = pain) |>
  ggplot(aes(x = 1, y = pain)) +
  geom_violin(fill = "#a6cee3", alpha = 0.6, width = 0.6) +
  geom_jitter(width = 0.06, size = 2.5, colour = "#1f78b4", alpha = 0.8) +
  stat_summary(fun = median, geom = "crossbar",
               colour = "firebrick", width = 0.3, linewidth = 0.9) +
  geom_hline(yintercept = 5, linetype = "dashed", colour = "grey30") +
  annotate("text", x = 1.32, y = 5.15, label = "Reference median = 5",
           hjust = 0, size = 3.5) +
  scale_y_continuous(breaks = 0:10) +
  labs(title = "One-Sample Wilcoxon: NRS Pain vs Reference",
       subtitle = "Violin = distribution; red bar = sample median; dashed = reference",
       y = "NRS Pain Score (0–10)", x = NULL) +
  theme_bw(base_size = 13) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())

import numpy as np, matplotlib.pyplot as plt
from scipy import stats
from scipy.stats import gaussian_kde
 
pain = np.array([7,8,6,9,5,7,8,7,6,9,7,8,5,7,9,6,8,7])
res  = stats.wilcoxon(pain - 5)
print(f"W={res.statistic:.0f}, p={res.pvalue:.4f}")
 
rng = np.random.default_rng(42)
fig, ax = plt.subplots(figsize=(5, 6))
y_range = np.linspace(3, 11, 200)
kde     = gaussian_kde(pain, bw_method=0.4)
width   = kde(y_range) * 2.5
ax.fill_betweenx(y_range, -width, width, color='#a6cee3', alpha=0.55)
ax.scatter(rng.uniform(-0.06, 0.06, len(pain)), pain, color='#1f78b4', s=50, alpha=0.8, zorder=3)
ax.hlines(np.median(pain), -0.28, 0.28, color='firebrick', lw=2.5, zorder=4, label=f"Median={np.median(pain):.0f}")
ax.axhline(5, linestyle='--', color='grey', lw=1.5, label="Reference=5")
ax.set_xlim(-0.5, 0.7); ax.set_xticks([])
ax.set_ylabel("NRS Pain Score (0–10)", fontsize=12)
ax.set_title(f"Wilcoxon Signed-Rank vs Reference\nW={res.statistic:.0f}, p={res.pvalue:.3f}", fontsize=12)
ax.legend(loc='upper right'); ax.grid(axis='y', alpha=0.3)
plt.tight_layout(); plt.savefig("wilcoxon_onesample.png", dpi=150, bbox_inches='tight'); plt.show()

# 15 readmissions in 250 procedures; reference = 4%
prop.test(x = 15, n = 250, p = 0.04, correct = FALSE)

from statsmodels.stats.proportion import proportions_ztest, proportion_confint
z_stat, p_val = proportions_ztest(count=15, nobs=250, value=0.04)
ci = proportion_confint(15, 250, alpha=0.05, method='normal')
print(f"z={z_stat:.3f}, p={p_val:.4f}, 95% CI: {ci[0]:.3f}–{ci[1]:.3f}")

library(tidyverse)
obs      <- c(A = 96, B = 16, AB = 6, O = 82)
expected <- c(A = 0.42, B = 0.10, AB = 0.04, O = 0.44)
chisq.test(obs, p = expected)
 
tibble(Group = names(obs),
       Observed = as.numeric(obs),
       Expected = expected * sum(obs)) |>
  pivot_longer(-Group, names_to = "type", values_to = "count") |>
  ggplot(aes(x = Group, y = count, fill = type)) +
  geom_col(position = "dodge", width = 0.6) +
  scale_fill_manual(values = c("Observed" = "#4575b4", "Expected" = "#d73027")) +
  labs(title    = "Chi-Square GoF: ABO Blood Types",
       subtitle = "Cardiac patients (n=200) vs UK population; χ²(3)=3.42, p=0.331",
       x = "Blood type", y = "Count", fill = NULL) +
  theme_bw(base_size = 13) + theme(legend.position = "top")

import numpy as np, matplotlib.pyplot as plt
from scipy.stats import chisquare
 
obs      = np.array([96, 16, 6, 82])
expected = np.array([0.42, 0.10, 0.04, 0.44]) * obs.sum()
groups   = ['A', 'B', 'AB', 'O']
chi2, p  = chisquare(obs, f_exp=expected)
print(f"chi2={chi2:.3f}, p={p:.4f}")
 
x = np.arange(len(groups))
fig, ax = plt.subplots(figsize=(7, 5))
ax.bar(x - 0.2, obs,      width=0.35, label='Observed', color='#4575b4', alpha=0.9)
ax.bar(x + 0.2, expected, width=0.35, label='Expected',  color='#d73027', alpha=0.9)
ax.set_xticks(x); ax.set_xticklabels(groups, fontsize=12)
ax.set_xlabel("Blood type"); ax.set_ylabel("Count")
ax.set_title(f"Chi-Square GoF: ABO Blood Types\nχ²(3)={chi2:.2f}, p={p:.3f}")
ax.legend(); ax.grid(axis='y', alpha=0.3)
plt.tight_layout(); plt.savefig("chisq_gof.png", dpi=150, bbox_inches='tight'); plt.show()

library(tidyverse); library(rstatix)
set.seed(42)
df <- tibble(group = rep(c("ACE inhibitor","Placebo"), each=45),
             sbp   = c(rnorm(45,12.4,8.2), rnorm(45,5.8,7.6)))
 
levene_test(df, sbp ~ group)
t_test(df, sbp ~ group, var.equal = TRUE, detailed = TRUE)
 
ggplot(df, aes(x=group, y=sbp, fill=group)) +
  geom_boxplot(alpha=0.6, outlier.shape=NA, width=0.5) +
  geom_jitter(aes(colour=group), width=0.12, size=1.8, alpha=0.6) +
  stat_compare_means(method="t.test", label="p.signif",
                     comparisons=list(c("ACE inhibitor","Placebo")), size=5) +
  scale_fill_manual(values=c("#4575b4","#d73027")) +
  scale_colour_manual(values=c("#4575b4","#d73027")) +
  labs(title="Independent t-Test: SBP Reduction by Treatment",
       subtitle="Mean ± IQR; individual patients as dots",
       x=NULL, y="SBP reduction (mmHg)") +
  theme_bw(base_size=13) + theme(legend.position="none")

import numpy as np, pandas as pd, matplotlib.pyplot as plt, seaborn as sns
from scipy import stats
 
rng  = np.random.default_rng(42)
ace  = rng.normal(12.4, 8.2, 45); plac = rng.normal(5.8, 7.6, 45)
lev_stat, lev_p = stats.levene(ace, plac)
t_stat, p_val   = stats.ttest_ind(ace, plac, equal_var=True)
print(f"Levene: F={lev_stat:.3f}, p={lev_p:.4f}")
print(f"t={t_stat:.3f}, p={p_val:.4f}")
 
df = pd.DataFrame({'group':['ACE inhibitor']*45+['Placebo']*45,
                   'sbp'  :np.concatenate([ace,plac])})
 
fig, ax = plt.subplots(figsize=(6,6))
palette = {'ACE inhibitor':'#4575b4','Placebo':'#d73027'}
sns.boxplot(data=df, x='group', y='sbp', palette=palette, width=0.45, fliersize=0, ax=ax)
sns.stripplot(data=df, x='group', y='sbp', palette=palette, jitter=0.1, size=4, alpha=0.5, ax=ax)
y_max = df['sbp'].max() + 3
ax.plot([0,0,1,1],[y_max,y_max+1,y_max+1,y_max], lw=1.2, color='black')
ax.text(0.5, y_max+1.5, f"p={p_val:.4f}", ha='center', fontsize=11)
ax.set_ylabel("SBP Reduction (mmHg)", fontsize=12); ax.set_xlabel("")
ax.set_title(f"Independent t-Test: SBP Reduction\nMean diff={np.mean(ace)-np.mean(plac):.1f} mmHg, t={t_stat:.2f}, p={p_val:.4f}", fontsize=11)
ax.grid(axis='y', alpha=0.3); plt.tight_layout()
plt.savefig("indep_ttest.png", dpi=150, bbox_inches='tight'); plt.show()

set.seed(42)
df_crp <- tibble(group = rep(c("Bacterial","Viral"), c(30,28)),
                 crp   = c(abs(rnorm(30,118.4,94.2)), abs(rnorm(28,22.6,18.7))))
t.test(crp ~ group, data=df_crp, var.equal=FALSE)
 
ggplot(df_crp, aes(x=group, y=crp+1, fill=group)) +
  geom_boxplot(alpha=0.65, width=0.5) +
  geom_jitter(aes(colour=group), width=0.12, size=1.8, alpha=0.6) +
  scale_y_log10(labels=scales::label_comma()) +
  scale_fill_manual(values=c("#e66101","#5e3c99")) +
  scale_colour_manual(values=c("#e66101","#5e3c99")) +
  labs(title="Welch's t-Test: CRP by Infection Type",
       subtitle="Log scale; Levene's p=0.001 → unequal variances → Welch's appropriate",
       x=NULL, y="CRP+1 (mg/L, log scale)") +
  theme_bw(base_size=13) + theme(legend.position="none")

import numpy as np, pandas as pd, matplotlib.pyplot as plt, seaborn as sns
from scipy import stats
 
rng       = np.random.default_rng(42)
bacterial = np.abs(rng.normal(118.4, 94.2, 30))
viral     = np.abs(rng.normal(22.6,  18.7, 28))
lev_stat, lev_p = stats.levene(bacterial, viral)
t_stat, p_val   = stats.ttest_ind(bacterial, viral, equal_var=False)
print(f"Levene p={lev_p:.4f}; Welch's t={t_stat:.3f}, p={p_val:.4f}")
 
df = pd.DataFrame({'Group':['Bacterial']*30+['Viral']*28,
                   'CRP'  :np.concatenate([bacterial, viral])})
fig, ax = plt.subplots(figsize=(6,6))
palette = {'Bacterial':'#e66101','Viral':'#5e3c99'}
sns.boxplot(data=df, x='Group', y='CRP', palette=palette, width=0.45, ax=ax)
sns.stripplot(data=df, x='Group', y='CRP', palette=palette, jitter=0.1, size=4, alpha=0.5, ax=ax)
ax.set_yscale('log'); ax.set_ylabel("CRP (mg/L, log scale)", fontsize=12)
ax.set_title(f"Welch's t-Test: CRP by Infection\nLevene p={lev_p:.3f}, t={t_stat:.2f}, p<0.001")
ax.grid(axis='y', alpha=0.3); plt.tight_layout()
plt.savefig("welch_ttest.png", dpi=150, bbox_inches='tight'); plt.show()

library(tidyverse)
std_care <- c(42,55,58,60,48,70,62,50,75,38,65,72,55,43,68,58,60,70,45,57,80,30)
int_prog <- c(62,72,80,65,75,82,70,68,75,73,65,78,85,62,70,80,72,68,75,82,66,78,62,74)
df_qol   <- tibble(group=rep(c("Standard care","Integrated programme"),c(22,24)),
                   qol  =c(std_care,int_prog))
wilcox.test(qol ~ group, data=df_qol)
 
ggplot(df_qol, aes(x=group, y=qol, fill=group)) +
  geom_violin(alpha=0.5, trim=FALSE) +
  geom_boxplot(width=0.15, outlier.shape=NA, alpha=0.8) +
  geom_jitter(aes(colour=group), width=0.06, size=2, alpha=0.6) +
  scale_fill_manual(values=c("#1b7837","#762a83")) +
  scale_colour_manual(values=c("#1b7837","#762a83")) +
  labs(title="Mann-Whitney U: Quality of Life by Care Group",
       subtitle="U=161.5, p=0.014; violin+box+jitter",
       x=NULL, y="EORTC QoL score (0–100)") +
  theme_bw(base_size=13) + theme(legend.position="none")

import numpy as np, pandas as pd, matplotlib.pyplot as plt, seaborn as sns
from scipy import stats
 
std_care = [42,55,58,60,48,70,62,50,75,38,65,72,55,43,68,58,60,70,45,57,80,30]
int_prog = [62,72,80,65,75,82,70,68,75,73,65,78,85,62,70,80,72,68,75,82,66,78,62,74]
u_stat, p_val = stats.mannwhitneyu(std_care, int_prog, alternative='two-sided')
print(f"U={u_stat:.1f}, p={p_val:.4f}")
 
df = pd.DataFrame({'Group':['Standard care']*22+['Integrated programme']*24,
                   'QoL'  :std_care+int_prog})
fig, ax = plt.subplots(figsize=(7,6))
palette = {'Standard care':'#1b7837','Integrated programme':'#762a83'}
sns.violinplot(data=df,x='Group',y='QoL',palette=palette,inner=None,alpha=0.45,ax=ax)
sns.boxplot(data=df,x='Group',y='QoL',palette=palette,width=0.15,fliersize=0,ax=ax)
sns.stripplot(data=df,x='Group',y='QoL',palette=palette,jitter=0.07,size=4,alpha=0.5,ax=ax)
ax.set_xlabel(""); ax.set_ylabel("EORTC QoL Score (0–100)", fontsize=12)
ax.set_title(f"Mann-Whitney U: QoL by Care Group\nU={u_stat:.0f}, p={p_val:.4f}")
ax.grid(axis='y', alpha=0.3); plt.tight_layout()
plt.savefig("mannwhitney.png", dpi=150, bbox_inches='tight'); plt.show()

library(tidyverse)
set.seed(42); n <- 20
pre  <- rnorm(n, 168.4, 35)
post <- pre - rnorm(n, 43.8, 28.4)
t.test(pre, post, paired=TRUE)
 
tibble(id=1:n, pre=pre, post=post) |>
  pivot_longer(c(pre,post), names_to="time", values_to="sodium") |>
  mutate(time=factor(time, c("pre","post"), c("Before","After"))) |>
  ggplot(aes(x=time, y=sodium)) +
  geom_line(aes(group=id), colour="grey60", alpha=0.5) +
  geom_point(aes(colour=time), size=2.5, alpha=0.8) +
  stat_summary(fun=mean, geom="line", aes(group=1), colour="firebrick", linewidth=2) +
  stat_summary(fun.data=mean_cl_normal, geom="pointrange",
               colour="firebrick", size=0.9, linewidth=1.5) +
  scale_colour_manual(values=c("Before"="#4575b4","After"="#d73027")) +
  labs(title="Paired t-Test: Urinary Sodium Before/After Restriction",
       subtitle="Grey=individual trajectories; red=mean ± 95% CI",
       x=NULL, y="24-h urinary sodium (mmol/24h)") +
  theme_bw(base_size=13) + theme(legend.position="none")

import numpy as np, matplotlib.pyplot as plt
from scipy import stats
 
rng  = np.random.default_rng(42); n = 20
pre  = rng.normal(168.4, 35, n)
post = pre - rng.normal(43.8, 28.4, n)
t_stat, p_val = stats.ttest_rel(pre, post)
print(f"t={t_stat:.3f}, p={p_val:.4f}")
 
fig, ax = plt.subplots(figsize=(6,6))
for i in range(n):
    ax.plot([0,1],[pre[i],post[i]], color='grey', alpha=0.4, lw=1)
ax.scatter([0]*n, pre,  color='#4575b4', s=50, alpha=0.75, zorder=3)
ax.scatter([1]*n, post, color='#d73027', s=50, alpha=0.75, zorder=3)
ci_pre  = stats.t.interval(0.95, n-1, np.mean(pre),  stats.sem(pre))
ci_post = stats.t.interval(0.95, n-1, np.mean(post), stats.sem(post))
ax.errorbar([0,1],[np.mean(pre),np.mean(post)],
            yerr=[[np.mean(pre)-ci_pre[0],np.mean(post)-ci_post[0]],
                  [ci_pre[1]-np.mean(pre),ci_post[1]-np.mean(post)]],
            fmt='o-', color='firebrick', markersize=9, lw=2.5, capsize=6, zorder=5)
ax.set_xticks([0,1]); ax.set_xticklabels(['Before','After'], fontsize=12)
ax.set_ylabel("24-h Urinary Sodium (mmol/24h)", fontsize=12)
ax.set_title(f"Paired t-Test: Sodium Restriction\nMean reduction={np.mean(pre-post):.1f} mmol/24h, p<0.001")
ax.grid(axis='y', alpha=0.3); plt.tight_layout()
plt.savefig("paired_ttest.png", dpi=150, bbox_inches='tight'); plt.show()

library(tidyverse)
pre  <- c(7,8,6,9,7,8,6,7,9,8,7,6,8,7,9,6)
post <- c(4,5,3,5,4,4,3,4,5,4,3,3,5,4,4,3)
wilcox.test(pre, post, paired=TRUE)
 
tibble(id=1:16, pre=pre, post=post, improved=ifelse(post<pre,"Improved","Not improved")) |>
  pivot_longer(c(pre,post), names_to="time", values_to="nrs") |>
  mutate(time=factor(time,c("pre","post"),c("Before","After"))) |>
  ggplot(aes(x=time, y=nrs, group=id)) +
  geom_line(aes(colour=rep(c("Improved","Not improved"),each=2)[match(id,1:16)]),
            alpha=0.7, linewidth=1) +
  geom_point(size=2.5) +
  scale_colour_manual(values=c("Improved"="#1b7837","Not improved"="#d73027"), name=NULL) +
  labs(title="Wilcoxon Signed-Rank: Pain Before/After Physiotherapy",
       subtitle="Green=improved; red=not improved",
       x=NULL, y="NRS Pain Score (0–10)") +
  theme_bw(base_size=13) + theme(legend.position="bottom")

import numpy as np, matplotlib.pyplot as plt
from scipy import stats
from matplotlib.lines import Line2D
 
pre  = np.array([7,8,6,9,7,8,6,7,9,8,7,6,8,7,9,6])
post = np.array([4,5,3,5,4,4,3,4,5,4,3,3,5,4,4,3])
res  = stats.wilcoxon(pre, post)
print(f"W={res.statistic:.1f}, p={res.pvalue:.4f}")
 
improved = post < pre
fig, ax = plt.subplots(figsize=(6,6))
for i in range(len(pre)):
    col = '#1b7837' if improved[i] else '#d73027'
    ax.plot([0,1],[pre[i],post[i]], color=col, alpha=0.65, lw=1.5)
    ax.scatter([0,1],[pre[i],post[i]], color=col, s=40, zorder=3)
ax.set_xticks([0,1]); ax.set_xticklabels(['Before','After'], fontsize=12)
ax.set_yticks(range(0,11)); ax.set_ylabel("NRS Pain Score (0–10)", fontsize=12)
ax.set_title(f"Wilcoxon Signed-Rank: Knee OA Physio\nW={res.statistic:.0f}, p={res.pvalue:.4f}")
legend_elements=[Line2D([0],[0],color='#1b7837',label=f'Improved ({improved.sum()})'),
                 Line2D([0],[0],color='#d73027',label=f'Not improved ({(~improved).sum()})')]
ax.legend(handles=legend_elements, loc='upper right')
ax.grid(axis='y', alpha=0.3); plt.tight_layout()
plt.savefig("wilcoxon_signed_rank.png", dpi=150, bbox_inches='tight'); plt.show()

library(tidyverse)
ct <- matrix(c(24,36,16,104), nrow=2,
             dimnames=list(c("Smoker","Non-smoker"),c("Pneumonia","No pneumonia")))
chisq.test(ct)
 
as_tibble(ct, rownames="Smoking") |>
  pivot_longer(-Smoking, names_to="Outcome", values_to="n") |>
  group_by(Smoking) |> mutate(pct=n/sum(n)) |>
  ggplot(aes(x=Smoking, y=pct, fill=Outcome)) +
  geom_col(position="fill", width=0.55) +
  geom_text(aes(label=scales::percent(pct,accuracy=0.1)),
            position=position_fill(vjust=0.5), size=4, colour="white") +
  scale_y_continuous(labels=scales::percent_format()) +
  scale_fill_manual(values=c("Pneumonia"="#d73027","No pneumonia"="#4575b4")) +
  labs(title="Chi-Square: Smoking and Post-Op Pneumonia",
       subtitle="χ²(1)=16.58, p<0.001; OR=4.33",
       x=NULL, y="Proportion", fill="Outcome") +
  theme_bw(base_size=13)

import numpy as np, matplotlib.pyplot as plt
from scipy.stats import chi2_contingency
 
ct = np.array([[24,36],[16,104]])
chi2, p, dof, expected = chi2_contingency(ct, correction=False)
print(f"χ²({dof})={chi2:.3f}, p={p:.6f}")
 
groups = ['Smoker','Non-smoker']
pneum  = ct[:,0] / ct.sum(axis=1)
no_pn  = ct[:,1] / ct.sum(axis=1)
fig, ax = plt.subplots(figsize=(6,6))
ax.bar(groups, pneum, label='Pneumonia', color='#d73027', alpha=0.9)
ax.bar(groups, no_pn, bottom=pneum, label='No pneumonia', color='#4575b4', alpha=0.9)
for i,(p0,p1) in enumerate(zip(pneum,no_pn)):
    ax.text(i, p0/2,     f"{p0:.1%}", ha='center', va='center', color='white', fontsize=11)
    ax.text(i, p0+p1/2,  f"{p1:.1%}", ha='center', va='center', color='white', fontsize=11)
ax.set_ylabel("Proportion"); ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda y,_:f'{y:.0%}'))
ax.set_title(f"Chi-Square: Smoking and Post-Op Pneumonia\nχ²(1)={chi2:.2f}, p<0.001")
ax.legend(loc='upper right'); ax.grid(axis='y', alpha=0.3)
plt.tight_layout(); plt.savefig("chisq_independence.png", dpi=150, bbox_inches='tight'); plt.show()

ct_small <- matrix(c(4,1,1,6), nrow=2,
                   dimnames=list(c("Immunosuppressed","Not"),c("Infection","No infection")))
fisher.test(ct_small)
# Expected cells: min=(5×5)/12=2.08 < 5 → Fisher's appropriate

from scipy.stats import fisher_exact
import numpy as np
ct_small = np.array([[4,1],[1,6]])
odds_ratio, p_val = fisher_exact(ct_small, alternative='two-sided')
print(f"OR={odds_ratio:.2f}, p={p_val:.4f}")

ct_mc <- matrix(c(38,12,22,8), nrow=2,
                dimnames=list(c("Pre: Compliant","Pre: Non-compliant"),
                              c("Post: Compliant","Post: Non-compliant")))
mcnemar.test(ct_mc)

from statsmodels.stats.contingency_tables import mcnemar
import numpy as np
ct_mc  = np.array([[38,12],[22,8]])
result = mcnemar(ct_mc, exact=False)
print(f"chi2={result.statistic:.3f}, p={result.pvalue:.4f}")

library(tidyverse); library(rstatix)
set.seed(42)
df_anova <- tibble(
  group = rep(c("Placebo","Low","Medium","High"), each=50),
  vomit = c(rnorm(50,6.8,2.4), rnorm(50,5.1,2.1),
            rnorm(50,3.4,1.8), rnorm(50,2.9,1.7))
) |> mutate(group=factor(group, levels=c("Placebo","Low","Medium","High")))
 
# Levene's and ANOVA
levene_test(df_anova, vomit ~ group)
anova_result <- aov(vomit ~ group, data=df_anova)
summary(anova_result)
 
# Tukey HSD post-hoc
tukey_hsd(df_anova, vomit ~ group)
 
# Plot
ggplot(df_anova, aes(x=group, y=vomit, fill=group)) +
  geom_boxplot(alpha=0.65, outlier.shape=NA, width=0.55) +
  geom_jitter(aes(colour=group), width=0.12, size=1.5, alpha=0.5) +
  stat_compare_means(method="anova", label.y=max(df_anova$vomit)+1) +
  scale_fill_brewer(palette="RdYlGn") +
  scale_colour_brewer(palette="RdYlGn") +
  labs(title="One-Way ANOVA: Vomiting Episodes by Drug Dose",
       subtitle="Post-hoc Tukey HSD; medium vs high not significant (p=0.41)",
       x="Treatment group", y="Vomiting episodes (24h)") +
  theme_bw(base_size=13) + theme(legend.position="none")

import numpy as np, pandas as pd, matplotlib.pyplot as plt, seaborn as sns
from scipy import stats
import scikit_posthocs as sp
 
rng = np.random.default_rng(42)
groups_dict = {'Placebo':rng.normal(6.8,2.4,50), 'Low':rng.normal(5.1,2.1,50),
               'Medium':rng.normal(3.4,1.8,50),  'High':rng.normal(2.9,1.7,50)}
df = pd.DataFrame([(g,v) for g,vals in groups_dict.items() for v in vals],
                  columns=['group','vomit'])
 
# ANOVA
f_stat, p_val = stats.f_oneway(*[v for v in groups_dict.values()])
print(f"F={f_stat:.2f}, p={p_val:.6f}")
 
# Tukey post-hoc via pingouin
import pingouin as pg
posthoc = pg.pairwise_tukey(data=df, dv='vomit', between='group')
print(posthoc[['A','B','diff','p-tukey']])
 
# Plot
fig, ax = plt.subplots(figsize=(8,6))
order  = ['Placebo','Low','Medium','High']
palette = {'Placebo':'#d73027','Low':'#fc8d59','Medium':'#91cf60','High':'#1a9850'}
sns.boxplot(data=df, x='group', y='vomit', order=order, palette=palette,
            width=0.5, fliersize=0, ax=ax)
sns.stripplot(data=df, x='group', y='vomit', order=order, palette=palette,
              jitter=0.1, size=3, alpha=0.4, ax=ax)
ax.set_xlabel("Treatment group"); ax.set_ylabel("Vomiting episodes (24h)")
ax.set_title(f"One-Way ANOVA: Vomiting Episodes by Dose\nF(3,196)={f_stat:.1f}, p<0.001")
ax.grid(axis='y', alpha=0.3); plt.tight_layout()
plt.savefig("oneway_anova.png", dpi=150, bbox_inches='tight'); plt.show()

library(tidyverse); library(rstatix)
set.seed(42)
df_kw <- tibble(
  hospital = rep(c("A","B","C"), c(35,38,33)),
  pain     = c(sample(c(3,4,5,6), 35, replace=TRUE, prob=c(0.3,0.4,0.2,0.1)),
               sample(c(4,5,6,7,8), 38, replace=TRUE, prob=c(0.15,0.25,0.3,0.2,0.1)),
               sample(c(3,4,5,6,7), 33, replace=TRUE, prob=c(0.2,0.25,0.3,0.15,0.1)))
)
 
kruskal.test(pain ~ hospital, data=df_kw)
dunn_test(df_kw, pain ~ hospital, p.adjust.method="bonferroni")
 
ggplot(df_kw, aes(x=hospital, y=pain, fill=hospital)) +
  geom_violin(alpha=0.55, trim=FALSE) +
  geom_boxplot(width=0.12, outlier.shape=NA) +
  geom_jitter(aes(colour=hospital), width=0.08, size=1.8, alpha=0.5) +
  scale_fill_brewer(palette="Set2") + scale_colour_brewer(palette="Set2") +
  scale_y_continuous(breaks=0:10) +
  labs(title="Kruskal-Wallis: Discharge Pain by Hospital",
       subtitle="H(2)=8.74, p=0.013; Dunn's post-hoc: A<B (p=0.010)",
       x="Hospital", y="NRS Pain at Discharge (0–10)") +
  theme_bw(base_size=13) + theme(legend.position="none")

import numpy as np, pandas as pd, matplotlib.pyplot as plt, seaborn as sns
from scipy import stats
import scikit_posthocs as sp
 
rng = np.random.default_rng(42)
data_A = rng.choice([3,4,5,6], 35, p=[0.3,0.4,0.2,0.1])
data_B = rng.choice([4,5,6,7,8], 38, p=[0.15,0.25,0.3,0.2,0.1])
data_C = rng.choice([3,4,5,6,7], 33, p=[0.2,0.25,0.3,0.15,0.1])
 
h_stat, p_val = stats.kruskal(data_A, data_B, data_C)
print(f"H={h_stat:.3f}, p={p_val:.4f}")
 
df = pd.DataFrame({'Hospital':['A']*35+['B']*38+['C']*33,
                   'Pain'    :np.concatenate([data_A,data_B,data_C])})
dunn_result = sp.posthoc_dunn(df, val_col='Pain', group_col='Hospital', p_adjust='bonferroni')
print(dunn_result)
 
fig, ax = plt.subplots(figsize=(7,6))
palette = {'A':'#66c2a5','B':'#fc8d62','C':'#8da0cb'}
sns.violinplot(data=df,x='Hospital',y='Pain',palette=palette,inner=None,alpha=0.5,ax=ax)
sns.boxplot(data=df,x='Hospital',y='Pain',palette=palette,width=0.12,fliersize=0,ax=ax)
sns.stripplot(data=df,x='Hospital',y='Pain',palette=palette,jitter=0.08,size=4,alpha=0.4,ax=ax)
ax.set_ylabel("NRS Pain at Discharge (0–10)"); ax.set_xlabel("Hospital")
ax.set_title(f"Kruskal-Wallis: Pain by Hospital\nH(2)={h_stat:.2f}, p={p_val:.3f}")
ax.grid(axis='y', alpha=0.3); plt.tight_layout()
plt.savefig("kruskal_wallis.png", dpi=150, bbox_inches='tight'); plt.show()

library(tidyverse); library(rstatix)
set.seed(42); n <- 30
df_rm <- tibble(
  id        = rep(1:n, 4),
  timepoint = rep(c("Baseline","3 months","6 months","12 months"), each=n),
  creatinine = c(rnorm(n,142,38), rnorm(n,138,36), rnorm(n,131,34), rnorm(n,128,33))
) |> mutate(timepoint=factor(timepoint, levels=c("Baseline","3 months","6 months","12 months")))
 
anova_test(df_rm, dv=creatinine, wid=id, within=timepoint)
 
# Mean ± CI line plot
df_rm |>
  group_by(timepoint) |>
  summarise(mean=mean(creatinine), se=sd(creatinine)/sqrt(n()),
            ci_lo=mean-qt(0.975,n()-1)*se, ci_hi=mean+qt(0.975,n()-1)*se) |>
  ggplot(aes(x=timepoint, y=mean, group=1)) +
  geom_ribbon(aes(ymin=ci_lo, ymax=ci_hi), fill="#74add1", alpha=0.3) +
  geom_line(colour="#4575b4", linewidth=1.5) +
  geom_point(colour="#4575b4", size=4) +
  labs(title="Repeated Measures ANOVA: Serum Creatinine over 12 Months",
       subtitle="Mean ± 95% CI; F(3,87)=8.43, p<0.001, η²=0.23",
       x="Time point", y="Serum creatinine (µmol/L)") +
  theme_bw(base_size=13)

import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy import stats
import pingouin as pg
 
rng = np.random.default_rng(42); n = 30
timepoints = ['Baseline','3 months','6 months','12 months']
means = [142,138,131,128]; sds = [38,36,34,33]
df = pd.concat([pd.DataFrame({'id':range(n),'timepoint':tp,
                              'creatinine':rng.normal(m,s,n)})
               for tp,m,s in zip(timepoints,means,sds)])
df['timepoint'] = pd.Categorical(df['timepoint'], categories=timepoints, ordered=True)
 
# RM-ANOVA via pingouin
rm_result = pg.rm_anova(data=df, dv='creatinine', within='timepoint', subject='id', detailed=True)
print(rm_result[['Source','F','p-unc','np2']])
 
# Mean ± CI
summary = df.groupby('timepoint')['creatinine'].agg(['mean','sem'])
summary['ci95'] = summary['sem'] * stats.t.ppf(0.975, n-1)
 
fig, ax = plt.subplots(figsize=(8,5))
x = range(len(timepoints))
ax.fill_between(x, summary['mean']-summary['ci95'], summary['mean']+summary['ci95'],
                alpha=0.25, color='#74add1')
ax.plot(x, summary['mean'], 'o-', color='#4575b4', linewidth=2, markersize=8)
ax.errorbar(x, summary['mean'], yerr=summary['ci95'],
            fmt='none', color='#4575b4', capsize=5, linewidth=1.5)
ax.set_xticks(x); ax.set_xticklabels(timepoints)
ax.set_ylabel("Serum creatinine (µmol/L)"); ax.set_xlabel("")
ax.set_title("Repeated Measures ANOVA: Creatinine over 12 Months\n"
             "Mean ± 95% CI; F(3,87)=8.43, p<0.001")
ax.grid(alpha=0.3); plt.tight_layout()
plt.savefig("rm_anova.png", dpi=150, bbox_inches='tight'); plt.show()

library(tidyverse)
set.seed(42); n <- 20
df_fr <- tibble(
  id        = rep(1:n, 3),
  timepoint = rep(c("Baseline","Week 4","Week 8"), each=n),
  nrs       = c(sample(6:9,n,replace=TRUE), sample(3:7,n,replace=TRUE),
                sample(1:5,n,replace=TRUE))
) |> mutate(timepoint=factor(timepoint,levels=c("Baseline","Week 4","Week 8")))
 
friedman.test(nrs ~ timepoint | id, data=df_fr)
 
# Spaghetti + median
medians <- df_fr |> group_by(timepoint) |> summarise(med=median(nrs))
 
ggplot(df_fr, aes(x=timepoint, y=nrs, group=id)) +
  geom_line(colour="grey60", alpha=0.5, linewidth=0.8) +
  geom_point(colour="grey60", size=1.5, alpha=0.5) +
  geom_line(data=medians, aes(y=med, group=1),
            colour="firebrick", linewidth=2.5, inherit.aes=FALSE,
            mapping=aes(x=timepoint)) +
  geom_point(data=medians, aes(x=timepoint, y=med),
             colour="firebrick", size=5, inherit.aes=FALSE) +
  labs(title="Friedman Test: Pain Scores over 8 Weeks of Biologic Therapy",
       subtitle="Grey=individual patients; red=median trajectory; χ²(2)=28.4, p<0.001",
       x=NULL, y="NRS Pain Score (0–10)") +
  theme_bw(base_size=13)