医学论文统计方法怎么选：t 检验、方差分析、卡方、非参数检验、相关与回归

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

R aov，气缸数放在前面（I 型）     P = 0.056
R aov，变速箱放在前面（I 型）     P = 4.9e-07
car::Anova type = 2 / pingouin   P = 0.056
car::Anova type = 3（contr.sum） P = 0.083
statsmodels anova_lm(typ=3)      P = 0.083   # 与 SPSS 默认 III 型对应

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

# 两独立组：airquality 6 月（n=9）与 8 月（n=26）臭氧浓度
aq <- subset(airquality, Month %in% c(6, 8) & !is.na(Ozone)); aq$Month <- factor(aq$Month)
t.test(Ozone ~ Month, data = aq)                    # R 默认 Welch
t.test(Ozone ~ Month, data = aq, var.equal = TRUE)  # Student t，对应 SPSS 第一行
wilcox.test(Ozone ~ Month, data = aq, conf.int = TRUE)
effectsize::hedges_g(Ozone ~ Month, data = aq, pooled_sd = FALSE)

# 配对：sleep 数据，同一受试者两种药物
w <- reshape(sleep, idvar = "ID", timevar = "group", direction = "wide")
t.test(w$extra.2, w$extra.1, paired = TRUE)
wilcox.test(w$extra.2, w$extra.1, paired = TRUE, exact = FALSE)

# 多组：单因素方差分析、Welch ANOVA、事后检验、Kruskal-Wallis + Dunn
fit <- aov(weight ~ group, data = PlantGrowth); summary(fit)
oneway.test(weight ~ group, data = PlantGrowth)     # Welch ANOVA
TukeyHSD(fit)
games_howell_test(PlantGrowth, weight ~ group)      # 方差不齐时的事后检验
effectsize::omega_squared(fit)
aq_all <- subset(airquality, !is.na(Ozone)); aq_all$Month <- factor(aq_all$Month)
kruskal.test(Ozone ~ Month, data = aq_all)
dunn_test(aq_all, Ozone ~ Month, p.adjust.method = "holm")

# 重复测量：Indometh，6 名受试者 5 个时间点
ind <- subset(as.data.frame(Indometh), time %in% c(0.5, 1, 2, 4, 8))
ind$Subject <- factor(as.character(ind$Subject)); ind$time <- factor(ind$time)
rm <- aov_ez(id = "Subject", dv = "conc", data = ind, within = "time")
summary(rm)                       # Mauchly 球形检验与 GG/HF 校正
friedman.test(conc ~ time | Subject, data = ind)

# 分类变量：卡方、Fisher、McNemar
tab <- table(infert$spontaneous > 0, infert$case)
chisq.test(tab)$expected          # 先看期望频数
chisq.test(tab, correct = FALSE); fisher.test(tab)
mcnemar.test(matrix(c(794, 86, 150, 570), 2))

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

# 多重比较校正
p <- c(0.001, 0.008, 0.012, 0.021, 0.035, 0.048, 0.09, 0.2, 0.41, 0.74)
round(cbind(p, bonf = p.adjust(p, "bonferroni"), holm = p.adjust(p, "holm"), BH = p.adjust(p, "BH")), 4)

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

# Python 3.12：pip install scipy statsmodels pingouin pandas lifelines
# 数据先在 R 中导出：write.csv(airquality, "airquality.csv", row.names = FALSE) 等
import numpy as np, pandas as pd, pingouin as pg
from scipy import stats
import statsmodels.formula.api as smf
from statsmodels.stats.multitest import multipletests
from statsmodels.stats.contingency_tables import mcnemar

aq = pd.read_csv("airquality.csv").dropna(subset=["Ozone"])
x6, x8 = aq.Ozone[aq.Month == 6], aq.Ozone[aq.Month == 8]
print(stats.ttest_ind(x6, x8, equal_var=False))   # scipy 默认 equal_var=True，必须显式写 False 才是 Welch
print(pg.ttest(x6, x8, correction=True))           # pingouin 默认 'auto'：只在样本量不等时用 Welch
print(stats.mannwhitneyu(x6, x8))

sleep = pd.read_csv("sleep.csv").pivot(index="ID", columns="group", values="extra")
print(stats.ttest_rel(sleep[2], sleep[1]))
print(stats.wilcoxon(sleep[2], sleep[1], method="approx", correction=True))  # 与 R 默认一致

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

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

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

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

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

library(survival)
labs <- c("bili", "chol", "albumin", "copper", "alk.phos", "ast", "trig", "platelet", "protime")
pb <- na.omit(pbc[, labs])                      # 276 例完整病例
pca <- prcomp(pb, scale. = TRUE)                # 必须标准化：各指标单位不同
summary(pca)                                    # 方差解释比例与累计比例
eig <- pca$sdev^2                               # 特征值
loadings <- sweep(pca$rotation, 2, pca$sdev, "*")   # 载荷 = 特征向量 × sqrt(特征值)，即 SPSS 成分矩阵
round(loadings[, 1:2], 3)
# 平行分析：保留特征值大于随机数据 95 分位数的主成分
set.seed(2026)
sim <- replicate(500, eigen(cor(matrix(rnorm(nrow(pb) * ncol(pb)), nrow(pb))))$values)
rbind(observed = round(eig, 3), random95 = round(apply(sim, 1, quantile, 0.95), 3))

import numpy as np, pandas as pd
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
pb = pd.read_csv("pbc_labs.csv")                # R 中 write.csv(na.omit(pbc[, labs]), ...) 导出
pca = PCA().fit(StandardScaler().fit_transform(pb))
print(pca.explained_variance_ratio_.round(4))   # 与 R summary(pca) 的 Proportion of Variance 相同
eig = np.linalg.eigvalsh(np.corrcoef(pb.T.to_numpy()))[::-1]   # 相关矩阵特征值，与 R、SPSS 相同
loadings = pca.components_.T * np.sqrt(eig)     # 主成分符号可能与 R 相反，解释时看相对方向
print(pd.DataFrame(loadings[:, :2], index=pb.columns, columns=["PC1", "PC2"]).round(3))