04 · 统计推断基础

为什么需要检验统计假设？

统计方法如同工具，每种工具都有其适用范围。当我们使用 \(t\) 检验比较两组均值时，实际上是在构建一个基于正态分布理论的统计量。如果数据严重偏离正态分布（如极度偏态或存在异常值），或者两组方差差异悬殊，那么基于 \(t\) 分布计算的 \(p\) 值可能不再可靠，导致错误的统计推断。因此，在进行正式分析之前，检验数据是否满足统计方法的假设条件，是保证分析结果科学性的重要步骤。

残差（Residual）指观测值与模型预测值之间的差异，即 其中 \(y_i\) 为第 \(i\) 个观测值，\(\hat{y}_i\) 为模型预测值。例如，在比较两组均值时，残差是每位患者的实际观测值减去所在组的样本均数；在回归分析中，残差是实际值减去回归线预测值。残差反映了模型未能解释的随机波动，检验其正态性和方差齐性是验证模型假设的重要步骤。

方差齐性（Homoscedasticity）指不同组或不同预测值水平下的残差具有相同的方差。对于两组比较，即 \(\sigma_1^2 = \sigma_2^2\)；对于回归分析，即残差方差不随预测值变化。方差不齐时，标准误会估计不准，导致检验效能下降或第一类错误率偏离名义水平。实务中可用 Levene 检验或 Bartlett 检验判断方差齐性，若方差不齐可采用 Welch 校正或变量变换。

动手练习：检验正态性与方差齐性

使用 Melanoma 数据集，检验肿瘤厚度（thickness）在不同溃疡组（ulcer）中的分布特征：

library(MASS)
data(Melanoma)

# 1. 按溃疡分组计算描述统计
by(Melanoma$thickness, Melanoma$ulcer, summary)

# 2. 绘制两组直方图，观察分布形态
par(mfrow = c(1, 2))
hist(Melanoma$thickness[Melanoma$ulcer == 0], main = "无溃疡组", xlab = "厚度")
hist(Melanoma$thickness[Melanoma$ulcer == 1], main = "有溃疡组", xlab = "厚度")

# 3. 正态性检验（Shapiro-Wilk 检验）
shapiro.test(Melanoma$thickness[Melanoma$ulcer == 0])
shapiro.test(Melanoma$thickness[Melanoma$ulcer == 1])

# 4. 方差齐性检验（Levene 检验，需安装 car 包）
# install.packages("car")
library(car)
leveneTest(thickness ~ factor(ulcer), data = Melanoma)

思考题：两组厚度分布是否近似正态？方差是否齐？若违背假设，你会选择什么检验方法比较两组均值？

输出结果与解读：

估计量的优良性质

一个好的估计量应具备以下性质：

• 无偏性（Unbiasedness）：估计量的期望等于被估计的参数，即 \(E(\hat{\theta}) = \theta\)。无偏性保证了长期重复抽样下，估计不会系统性地高估或低估。
• 有效性（Efficiency）：在所有无偏估计量中，方差最小的估计量最有效。样本均值是总体均值的最有效线性无偏估计（BLUE）。
• 一致性（Consistency）：随着样本量增大，估计量依概率收敛于真实参数值。大样本下，估计误差趋于零。
• 充分性（Sufficiency）：估计量包含了样本中关于参数的全部信息。

课堂讨论：自由度为什么是 n-1？

在计算样本方差时，我们使用 \(n-1\) 而非 \(n\) 作为分母，这涉及"自由度"的概念。当我们用样本均值 \(\bar{x}\) 估计总体均值后，\(n\) 个偏差 \((x_i - \bar{x})\) 中仅有 \(n-1\) 个是独立的，因为它们的和恒为零（\(\sum(x_i - \bar{x}) = 0\)）。使用 \(n-1\) 作为分母，可以使样本方差成为总体方差的无偏估计。若使用 \(n\) 作分母，会系统性地低估总体方差，尤其在样本量较小时偏差更明显。

集中趋势反映数据的"中心位置"或"典型水平"，是描述一组数据向某一中心值靠拢的程度。均值是最常用的集中趋势指标，适用于近似对称分布；当数据存在极端值或明显偏态时，中位数能更稳健地代表"典型"水平。

离散程度反映数据的分散性或变异性，描述个体观测值偏离中心的幅度。标准差与原始数据同单位，便于解释；方差是标准差的平方，在统计推导中更常用（如方差分析）。离散程度大意味着个体差异大、数据分布分散；离散程度小则表明数据较为集中、个体差异小。

动手练习：计算均值、方差与标准差

使用 Melanoma 数据集，手动计算肿瘤厚度（thickness）的样本均值、方差和标准差，并与 R 内置函数结果对比：

library(MASS)
data(Melanoma)

x <- Melanoma$thickness
n <- length(x)

# 手动计算
x_bar <- sum(x) / n                    # 样本均值
s2 <- sum((x - x_bar)^2) / (n - 1)     # 样本方差（无偏）
s <- sqrt(s2)                          # 样本标准差

cat("手动计算：\n")
cat("均值 =", round(x_bar, 4), "\n")
cat("方差 =", round(s2, 4), "\n")
cat("标准差 =", round(s, 4), "\n\n")

# 与 R 内置函数对比
cat("R 内置函数：\n")
cat("均值 =", round(mean(x), 4), "\n")
cat("方差 =", round(var(x), 4), "\n")
cat("标准差 =", round(sd(x), 4), "\n")

输出结果：

思考题：

• 为什么方差要用 \(n-1\) 而不是 \(n\) 作分母？
• 如果数据中存在一个极端大值（如 50 mm），均值和标准差会如何变化？中位数呢？

抽样分布（Sampling Distribution）

抽样分布是指统计量（如样本均值）在所有可能样本中的概率分布。设想从同一总体中重复抽取无数个大小为 \(n\) 的样本，计算每个样本的均值，这些均值的分布就是样本均值的抽样分布。中心极限定理告诉我们，无论总体分布如何，当 \(n\) 足够大时，样本均值的抽样分布近似正态分布，其均值为总体均值 \(\mu\)，标准差为总体标准差除以 \(\sqrt{n}\)。

影响置信区间宽度的因素

• 置信水平：置信水平越高（如从 95% 提高到 99%），区间越宽。因为需要更大的把握覆盖真实参数。
• 样本量：样本量越大，标准误越小，区间越窄。要使区间宽度减半，需要将样本量增加到原来的 4 倍。
• 数据变异：数据本身的变异（标准差）越大，区间越宽。这反映了估计的内在难度。

x <- c(118, 122, 125, 119, 121, 124, 120)
t.test(x)$conf.int   # 默认 95% CI for mu

# 计算不同置信水平下的区间
t.test(x, conf.level = 0.90)$conf.int  # 90% CI（更窄）
t.test(x, conf.level = 0.99)$conf.int  # 99% CI（更宽）

动手练习：计算置信区间

使用 Melanoma 数据集，计算肿瘤厚度总体均值的 95% 置信区间，并比较不同置信水平下区间的宽度：

library(MASS)
data(Melanoma)

# 计算肿瘤厚度的描述统计
x <- Melanoma$thickness
cat("样本量：", length(x), "\n")
cat("样本均值：", round(mean(x), 3), "mm\n")
cat("样本标准差：", round(sd(x), 3), "mm\n")
cat("标准误：", round(sd(x) / sqrt(length(x)), 4), "mm\n\n")

# 计算 95% 置信区间
ci95 <- t.test(x)$conf.int
cat("95% CI: [", round(ci95[1], 3), ",", round(ci95[2], 3), "]\n")
cat("区间宽度：", round(ci95[2] - ci95[1], 3), "mm\n\n")

# 计算 90% 和 99% 置信区间比较
ci90 <- t.test(x, conf.level = 0.90)$conf.int
ci99 <- t.test(x, conf.level = 0.99)$conf.int
cat("90% CI: [", round(ci90[1], 3), ",", round(ci90[2], 3), "] 宽度：", round(ci90[2] - ci90[1], 3), "\n")
cat("99% CI: [", round(ci99[1], 3), ",", round(ci99[2], 3), "] 宽度：", round(ci99[2] - ci99[1], 3), "\n")

输出结果：

思考题：

• 如果样本量从 205 增加到 1000，置信区间的宽度会如何变化？
• 如何理解"95% 置信区间"的正确含义？为什么不能说"真实参数有 95% 的概率落在这个区间内"？

假设检验的基本步骤

1. 建立假设：明确原假设 \(H_0\)（通常表示"无效应"或"无差异"）和备择假设 \(H_1\)（表示研究者想要证明的效应或差异）。
2. 选择检验统计量：根据数据类型、样本量和研究设计选择合适的统计量（如 \(t\) 统计量、\(F\) 统计量、\(\chi^2\) 统计量等）。
3. 确定显著性水平：通常取 \(\alpha = 0.05\)，表示愿意承担 5% 的第一类错误风险。
4. 计算 \(p\) 值：在原假设成立的前提下，观察到当前样本或更极端样本的概率。
5. 做出决策：若 \(p < \alpha\)，拒绝 \(H_0\)；否则，不拒绝 \(H_0\)。
6. 解释结果：将统计结论转化为实际问题的语言。

选择检验方法的决策树

• 两组比较：
  • 独立样本 + 正态 + 方差齐：标准 \(t\) 检验（var.equal = TRUE）
  • 独立样本 + 正态 + 方差不齐：Welch \(t\) 检验（默认推荐）
  • 独立样本 + 非正态：Mann-Whitney \(U\) 检验
  • 配对样本：配对 \(t\) 检验（正态）或 Wilcoxon 符号秩检验（非正态）
• 多组比较（k ≥ 3）：
  • 正态 + 方差齐：单因素 ANOVA
  • 非正态或方差不齐：Kruskal-Wallis 检验

library(MASS)
data(Melanoma)

# 以溃疡分组比较肿瘤厚度（两组独立样本）
t.test(thickness ~ ulcer, data = Melanoma)           # Welch
wilcox.test(thickness ~ ulcer, data = Melanoma)      # 非参数

# 按年份分组比较厚度（多组比较示例）
Melanoma$year_grp <- cut(Melanoma$year, breaks = c(1960, 1970, 1975, 1980))
summary(aov(thickness ~ year_grp, data = Melanoma))

动手练习：两组比较检验

使用 Melanoma 数据集，比较不同溃疡组的肿瘤厚度，练习选择合适的检验方法：

library(MASS)
data(Melanoma)

# 1. 先检验正态性和方差齐性（回顾 1.1 节内容）
# 正态性检验
shapiro.test(Melanoma$thickness[Melanoma$ulcer == 0])
shapiro.test(Melanoma$thickness[Melanoma$ulcer == 1])

# 方差齐性检验
library(car)
leveneTest(thickness ~ factor(ulcer), data = Melanoma)

# 2. 由于正态性和方差齐性均被违背，使用 Welch t 检验
welch_result <- t.test(thickness ~ ulcer, data = Melanoma)
print(welch_result)

# 3. 同时使用非参数方法 Mann-Whitney U 检验
wilcox_result <- wilcox.test(thickness ~ ulcer, data = Melanoma)
print(wilcox_result)

# 4. 计算效应量（Cohen's d）
# install.packages("effectsize")
library(effectsize)
cohens_d <- cohens_d(thickness ~ ulcer, data = Melanoma)
print(cohens_d)

输出结果：

思考题：

• Welch t 检验和标准 t 检验有什么区别？为什么推荐优先使用 Welch t 检验？
• p 值很小是否意味着效应很大？为什么需要报告效应量？

Spearman 秩相关 vs Pearson 相关

当数据不满足正态性或存在明显异常值时，可使用 Spearman 秩相关系数。Spearman 相关系数基于变量的秩次计算，不依赖正态假设，对异常值稳健。它度量的是单调关系（不一定是线性的）。在 R 中：cor.test(x, y, method = "spearman")。

回归诊断：检查模型假设

线性回归的有效性依赖于以下假设：

• 线性性：\(Y\) 与 \(X\) 的关系是线性的。检查方法：残差 vs 拟合值图应呈随机分布，无系统模式。
• 独立性：残差相互独立。检查方法：Durbin-Watson 检验（时间序列数据）。
• 正态性：残差服从正态分布。检查方法：残差 Q-Q 图应近似直线。
• 方差齐性：残差方差不随预测值变化。检查方法：残差 vs 拟合值图中残差散布应均匀。

若假设不满足，可考虑变量变换（如对数变换）、加权最小二乘或使用稳健回归方法。

library(MASS)
data(Melanoma)

# 肿瘤厚度与年龄的相关性
cor.test(Melanoma$thickness, Melanoma$age, method = "pearson")

# 线性回归：厚度预测年龄
fit <- lm(age ~ thickness, data = Melanoma)
summary(fit)
confint(fit, "thickness", level = 0.95)

动手练习：相关与回归分析

使用 Melanoma 数据集，分析肿瘤厚度与年龄之间的关系：

library(MASS)
data(Melanoma)

# 1. 绘制散点图观察关系
plot(Melanoma$age, Melanoma$thickness,
     main = "肿瘤厚度与年龄的关系",
     xlab = "年龄（岁）", ylab = "厚度（mm）",
     pch = 19, col = rgb(0, 0.5, 0.5, 0.5))

# 2. 计算 Pearson 相关系数
pearson_cor <- cor.test(Melanoma$thickness, Melanoma$age, method = "pearson")
print(pearson_cor)

# 3. 计算 Spearman 秩相关（对异常值稳健）
spearman_cor <- cor.test(Melanoma$thickness, Melanoma$age, method = "spearman")
print(spearman_cor)

# 4. 简单线性回归
fit <- lm(thickness ~ age, data = Melanoma)
summary(fit)

# 5. 回归诊断
par(mfrow = c(2, 2))
plot(fit)

输出结果：

思考题：

• 相关系数显著是否意味着临床意义显著？如何解读 R² = 0.04？
• Pearson 相关系数和 Spearman 秩相关系数有什么区别？何时应该使用 Spearman？
• "相关不等于因果"——请举例说明两个变量可能相关但不存在因果关系的情况。

分类变量的类型

• 二分类变量（Binary）：只有两个类别，如性别（男/女）、生存状态（存活/死亡）、疾病诊断（阳性/阴性）。
• 名义变量（Nominal）：多个类别且无顺序关系，如血型（A/B/AB/O）、疾病类型、医院科室。
• 有序变量（Ordinal）：多个类别且有顺序关系，如疾病分期、疼痛程度（轻/中/重）、疗效评价（显效/有效/无效）。

不同类型的分类变量需要采用不同的统计方法。名义变量不能计算"平均"或"中位数"，有序变量则可以利用其顺序信息。

有序变量的特殊处理

有序变量除了报告频数和构成比外，还可报告：

• 中位数秩：将有序类别赋值为秩次后计算中位数。
• 累积百分比：从最低类别开始累积的百分比，如"疗效达有效及以上者占 75%"。
• 比例优势模型：利用有序信息进行更有效的统计分析。

注意：有序变量不应简单按名义变量处理，否则会损失信息。

动手练习：分类变量描述统计

使用 Melanoma 数据集，对分类变量进行描述统计分析：

library(MASS)
data(Melanoma)

# 1. 查看数据结构
str(Melanoma)

# 2. 性别分布（二分类变量）
table(Melanoma$sex)  # 频数
prop.table(table(Melanoma$sex)) * 100  # 构成比（%）

# 3. 溃疡状态分布（二分类变量）
table(Melanoma$ulcer)
prop.table(table(Melanoma$ulcer)) * 100

# 4. 生存状态（多分类，需理解编码）
# status: 1=死于黑色素瘤, 2=存活, 3=死于其他原因
table(Melanoma$status)
prop.table(table(Melanoma$status)) * 100

# 5. 创建有序变量示例：按厚度分组
Melanoma$thickness_grp <- cut(Melanoma$thickness,
                              breaks = c(0, 2, 5, Inf),
                              labels = c("薄", "中", "厚"),
                              right = FALSE)
table(Melanoma$thickness_grp)
prop.table(table(Melanoma$thickness_grp)) * 100

# 6. 交叉表：溃疡状态 × 厚度分组
table(Melanoma$ulcer, Melanoma$thickness_grp)
prop.table(table(Melanoma$ulcer, Melanoma$thickness_grp), 1) * 100  # 行百分比

输出结果：

思考题：

• 频数和构成比有什么区别？报告时应该同时报告哪些信息？
• 有序变量和名义变量在描述统计时有什么不同？为什么有序变量可以报告累积百分比？

临床案例：率的估计与报告

某临床试验中，50 例患者接受新疗法，其中 38 例有效。估计有效率的 95% 置信区间：

• 点估计：\(\hat{p} = 38/50 = 76\%\)
• Wilson 95% CI：约 62% ~ 86%

报告方式："新疗法的有效率为 76%（95% CI: 62%–86%）"。置信区间提供了估计不确定性的量化，是临床研究报告的标准要求。

# 不同方法计算率的置信区间
prop.test(38, 50, correct = FALSE)  # Wilson 区间
binom.test(38, 50)                   # Clopper-Pearson 精确区间

# 多个率的同时估计
prop.test(c(38, 25), c(50, 50))      # 比较两组率

动手练习：率的置信区间计算

使用 Melanoma 数据集，计算死亡率和溃疡发生率的置信区间：

library(MASS)
data(Melanoma)

# 1. 计算黑色素瘤相关死亡率（status = 1）
n_death <- sum(Melanoma$status == 1)
n_total <- nrow(Melanoma)
cat("死亡数：", n_death, "，总数：", n_total, "\n")
cat("死亡率：", round(n_death / n_total * 100, 1), "%\n\n")

# 2. 使用 Wilson 方法计算 95% CI
wilson_ci <- prop.test(n_death, n_total, correct = FALSE)
cat("Wilson 95% CI: [", round(wilson_ci$conf.int[1] * 100, 1), "%, ",
    round(wilson_ci$conf.int[2] * 100, 1), "%]\n\n")

# 3. 使用 Clopper-Pearson 精确方法
cp_ci <- binom.test(n_death, n_total)
cat("Clopper-Pearson 95% CI: [", round(cp_ci$conf.int[1] * 100, 1), "%, ",
    round(cp_ci$conf.int[2] * 100, 1), "%]\n\n")

# 4. 比较不同置信水平
wilson_90 <- prop.test(n_death, n_total, correct = FALSE, conf.level = 0.90)
wilson_99 <- prop.test(n_death, n_total, correct = FALSE, conf.level = 0.99)
cat("90% CI: [", round(wilson_90$conf.int[1] * 100, 1), "%, ",
    round(wilson_90$conf.int[2] * 100, 1), "%]\n")
cat("95% CI: [", round(wilson_ci$conf.int[1] * 100, 1), "%, ",
    round(wilson_ci$conf.int[2] * 100, 1), "%]\n")
cat("99% CI: [", round(wilson_99$conf.int[1] * 100, 1), "%, ",
    round(wilson_99$conf.int[2] * 100, 1), "%]\n")

输出结果：

思考题：

• 为什么 Wilson 区间比 Wald 区间更适合小样本？
• 如果死亡率只有 5%（n=10），Wilson 和 Clopper-Pearson 区间会有什么差异？

\(\chi^2\) 检验 vs Fisher 精确检验

特征	\(\chi^2\) 检验	Fisher 精确检验
理论基础	大样本渐近分布	精确概率计算
适用样本量	大样本	任意样本量
计算复杂度	简单	大表时计算量大
推荐使用	期望频数足够时	期望频数小或精确推断

tb <- matrix(c(20, 15, 10, 25), nrow = 2, byrow = TRUE)
chisq.test(tb)
fisher.test(tb)

动手练习：卡方检验与 Fisher 精确检验

使用 Melanoma 数据集，检验溃疡状态与性别之间是否存在关联：

library(MASS)
data(Melanoma)

# 1. 创建列联表：溃疡状态 × 性别
tb <- table(Melanoma$ulcer, Melanoma$sex)
rownames(tb) <- c("无溃疡", "有溃疡")
colnames(tb) <- c("男性", "女性")
print("列联表：")
print(tb)

# 2. 计算期望频数
chisq_test <- chisq.test(tb)
print("期望频数：")
print(round(chisq_test$expected, 1))

# 3. 检查期望频数条件
print(paste("最小期望频数：", round(min(chisq_test$expected), 2)))
print(paste("期望频数 < 5 的格子数：", sum(chisq_test$expected < 5)))

# 4. 卡方检验
print("卡方检验结果：")
print(chisq_test)

# 5. Fisher 精确检验（对比）
fisher_test <- fisher.test(tb)
print("Fisher 精确检验结果：")
print(fisher_test)

输出结果：

思考题：

• 什么情况下应该使用 Fisher 精确检验而不是卡方检验？
• 卡方检验的期望频数条件是什么？如果不满足应该怎么办？

OR 与 RR 的选择原则

• 队列研究 / RCT：优先报告 RR，更直观易解释。
• 病例对照研究：只能估计 OR。
• Logistic 回归：系数对应 OR，适合各种研究设计。
• 罕见病：OR ≈ RR，可直接解释为风险比。

# Logistic 回归示例
library(MASS)
data(Melanoma)

# 创建二分类结局变量（以生存状态为例）
Melanoma$death <- ifelse(Melanoma$status == 1, 1, 0)

# 单变量 Logistic 回归
fit <- glm(death ~ ulcer, family = binomial, data = Melanoma)
summary(fit)

# 计算 OR 及 95% CI
exp(cbind(OR = coef(fit), confint(fit)))

动手练习：OR、RR 与 Logistic 回归

使用 Melanoma 数据集，分析溃疡状态与死亡风险的关系：

library(MASS)
data(Melanoma)

# 1. 创建二分类结局变量（死于黑色素瘤 = 1）
Melanoma$death <- ifelse(Melanoma$status == 1, 1, 0)

# 2. 创建 2×2 表计算 OR
tb <- table(Melanoma$ulcer, Melanoma$death)
rownames(tb) <- c("无溃疡", "有溃疡")
colnames(tb) <- c("存活", "死亡")
print("列联表：")
print(tb)

# 3. 手动计算 OR
# OR = (a/b) / (c/d) = (a*d) / (b*c)
a <- tb[2, 2]; b <- tb[2, 1]; c <- tb[1, 2]; d <- tb[1, 1]
OR_manual <- (a * d) / (b * c)
cat("手动计算的 OR:", round(OR_manual, 2), "\n\n")

# 4. 使用 fisher.test 计算 OR 和 CI
fisher_result <- fisher.test(tb)
cat("Fisher 检验给出的 OR:", round(fisher_result$estimate, 2), "\n")
cat("95% CI: [", round(fisher_result$conf.int[1], 2), ", ",
    round(fisher_result$conf.int[2], 2), "]\n\n")

# 5. Logistic 回归
fit <- glm(death ~ ulcer, family = binomial, data = Melanoma)
print(summary(fit))

# 6. 计算 OR 及 95% CI
or_ci <- exp(cbind(OR = coef(fit), confint(fit)))
print("Logistic 回归 OR 及 95% CI：")
print(round(or_ci, 3))

输出结果：

思考题：

• OR 和 RR 有什么区别？什么情况下 OR 可以近似等于 RR？
• Logistic 回归的优势是什么？如何调整混杂因素？

方法选择的核心原则

统计方法的选择主要取决于：

1. 因变量类型：这是选择方法的首要依据。因变量是连续、二分类、有序分类、名义分类、计数还是生存时间？
2. 自变量类型：自变量是连续还是分类？有多少个自变量？
3. 研究设计：是独立样本还是配对/重复测量？是横断面、队列还是病例对照？
4. 样本量：大样本可以使用渐近方法，小样本可能需要精确检验。

数据挖掘场景的特殊挑战

• 高维特征：基因芯片可同时检测数万个基因，影像组学可提取数百个特征，传统单变量分析面临多重检验问题。
• 模型选择偏差：在同一数据上反复调参、选特征会导致性能估计的"乐观偏倚"。
• 因果推断困难：观察性数据中的混杂因素可能导致虚假关联，需要更审慎的因果推断方法。
• 可重复性危机：许多数据挖掘发现无法在独立数据上复现，需要规范化的验证流程。

FWER vs FDR：如何选择？

• 控制 FWER（Bonferroni/Holm）：适合验证性研究，要求"即使有一个假阳性也不可接受"的场景，如确证性临床试验。
• 控制 FDR（BH）：适合探索性研究，可以容忍一定比例的假阳性以换取更高的发现率，如基因筛选、候选标志物发现。
• 实务建议：探索阶段用 FDR，验证阶段用 FWER；报告时应明确校正方法和阈值。

p <- c(0.001, 0.04, 0.08, 0.15, 0.5)
p.adjust(p, method = "bonferroni")
p.adjust(p, method = "BH")   # Benjamini–Hochberg FDR

动手练习：多重检验校正

模拟多重检验场景，比较不同校正方法的效果：

# 模拟 20 个基因的 p 值（假设其中 3 个有真实效应）
set.seed(123)

# 生成模拟 p 值：大部分来自无效应假设，少数来自有效应假设
p_values <- c(
  runif(17, 0, 1),      # 17 个无真实效应（均匀分布）
  runif(3, 0, 0.01)     # 3 个有真实效应（偏向小值）
)

# 按从小到大排序
p_sorted <- sort(p_values)

cat("原始 p 值（排序后）：\n")
print(round(p_sorted, 4))

# 1. 不校正：直接用 0.05 阈值
cat("\n不校正（α = 0.05）：", sum(p_sorted < 0.05), "个显著\n")

# 2. Bonferroni 校正
p_bonf <- p.adjust(p_sorted, method = "bonferroni")
cat("Bonferroni 校正后：", sum(p_bonf < 0.05), "个显著\n")
cat("阈值：", round(0.05 / 20, 4), "\n\n")

# 3. Holm 校正（逐步校正）
p_holm <- p.adjust(p_sorted, method = "holm")
cat("Holm 校正后：", sum(p_holm < 0.05), "个显著\n\n")

# 4. Benjamini-Hochberg FDR
p_bh <- p.adjust(p_sorted, method = "BH")
cat("BH FDR 校正后：", sum(p_bh < 0.05), "个显著\n")

# 5. 创建比较表
results <- data.frame(
  原始p值 = round(p_sorted, 4),
  Bonferroni = round(p_bonf, 4),
  Holm = round(p_holm, 4),
  BH_FDR = round(p_bh, 4)
)
print(results)

输出结果：

思考题：

• 如果有 20,000 个基因，Bonferroni 阈值是多少？这会带来什么问题？
• FWER 控制和 FDR 控制有什么区别？各适合什么场景？

为什么会出现乐观偏倚？

当我们用数据选择模型或特征时，实际上是让数据"说话"选择最适合它的模型。这种选择过程会利用数据中的随机噪声，导致选出的模型在当前数据上表现特别好，但这种"好"部分是噪声的功劳。当应用到新数据时，噪声不同，表现就会下降。

library(MASS)
data(Melanoma)

set.seed(1)
n <- nrow(Melanoma)
idx <- sample(seq_len(n), size = floor(0.75 * n))
train <- Melanoma[idx, ]
test  <- Melanoma[-idx, ]
# 仅在 train 上拟合，在 test 上算误差或 AUC 等

动手练习：训练/测试划分

使用 Melanoma 数据集，演示训练/测试划分及其重要性：

library(MASS)
data(Melanoma)

# 1. 设置随机种子保证可重复性
set.seed(42)

# 2. 随机划分训练集和测试集（75% 训练，25% 测试）
n <- nrow(Melanoma)
train_idx <- sample(1:n, size = floor(0.75 * n))
train_data <- Melanoma[train_idx, ]
test_data <- Melanoma[-train_idx, ]

cat("训练集样本量：", nrow(train_data), "\n")
cat("测试集样本量：", nrow(test_data), "\n\n")

# 3. 创建二分类结局
train_data$death <- ifelse(train_data$status == 1, 1, 0)
test_data$death <- ifelse(test_data$status == 1, 1, 0)

# 4. 在训练集上拟合 Logistic 回归
fit <- glm(death ~ thickness + age + ulcer, 
           family = binomial, data = train_data)
print(summary(fit))

# 5. 在训练集和测试集上分别评估
# 训练集预测
train_pred <- predict(fit, type = "response")
train_auc <- pROC::auc(train_data$death, train_pred)

# 测试集预测
test_pred <- predict(fit, newdata = test_data, type = "response")
test_auc <- pROC::auc(test_data$death, test_pred)

cat("\n模型性能比较：\n")
cat("训练集 AUC：", round(train_auc, 3), "\n")
cat("测试集 AUC：", round(test_auc, 3), "\n")
cat("差异：", round(train_auc - test_auc, 3), "\n")

输出结果：

思考题：

• 为什么测试集的性能通常低于训练集？
• 如果在测试集上反复调参，会发生什么问题？

定位：因果推断（含 PSM 等方法）常被概括为医学数据挖掘从"描述现状"走向"指导治疗与实践决策"的最高阶段：要回答的是“若采取某干预，结局可能如何”，而不仅是“变量之间是否相关”。达到这一阶段需要在观察性数据中审慎处理混杂、明确识别假设，因而也是整个流程里最难做扎实的一环。

# 需安装 MatchIt；dat 含 treat(0/1)、结局与协变量
# library(MatchIt)
# m <- matchit(treat ~ age + sex + bmi, data = dat, method = "nearest", ratio = 1)
# md <- match.data(m)
# 在 md 上做结局比较或加权回归

动手练习：倾向得分匹配概念

理解倾向得分的计算和匹配原理：

library(MASS)
data(Melanoma)

# 1. 定义"处理"为溃疡状态，结局为死亡
Melanoma$treat <- Melanoma$ulcer  # 0=无溃疡, 1=有溃疡
Melanoma$death <- ifelse(Melanoma$status == 1, 1, 0)

# 2. 计算倾向得分（用 Logistic 回归估计）
ps_model <- glm(treat ~ age + sex, family = binomial, data = Melanoma)
Melanoma$ps <- predict(ps_model, type = "response")

# 3. 查看倾向得分分布
summary(Melanoma$ps)

# 4. 按倾向得分分组查看处理分配
cat("\n倾向得分分布（按处理组）：\n")
by(Melanoma$ps, Melanoma$treat, summary)

# 5. 检查共同支撑（重叠）
cat("\n共同支撑检查：\n")
cat("处理组 PS 范围：[", round(min(Melanoma$ps[Melanoma$treat==1]), 3), ", ",
    round(max(Melanoma$ps[Melanoma$treat==1]), 3), "]\n")
cat("对照组 PS 范围：[", round(min(Melanoma$ps[Melanoma$treat==0]), 3), ", ",
    round(max(Melanoma$ps[Melanoma$treat==0]), 3), "]\n")

# 6. 未调整的死亡风险比较
cat("\n未调整的比较：\n")
cat("处理组死亡率：", round(mean(Melanoma$death[Melanoma$treat==1]), 3), "\n")
cat("对照组死亡率：", round(mean(Melanoma$death[Melanoma$treat==0]), 3), "\n")
cat("粗差异：", round(mean(Melanoma$death[Melanoma$treat==1]) - 
               mean(Melanoma$death[Melanoma$treat==0]), 3), "\n")

输出结果：

思考题：

• 什么是倾向得分？为什么它可以用于因果推断？
• PSM 的关键假设是什么？"无未测混杂"假设为什么无法直接检验？

置换检验的优势

• 无需分布假设：不依赖正态分布等假设，适合偏态数据或小样本。
• 灵活性强：可以使用任意检验统计量，不限于均值差，也可以是中位数差、复杂评分等。
• 处理复杂流程：可以把“选特征+算统计量”整体包进置换循环，部分缓解逐步挑选带来的名义水平失真。
• 直观易懂：直接模拟零假设下的数据分布，结果易于解释。

set.seed(2)
x <- c(rnorm(30, 0), rnorm(30, 0.4))   # 两组各 30，略有位置差
g <- rep(1:2, each = 30)
obs <- abs(mean(x[g == 2]) - mean(x[g == 1]))
B <- 4999
perm <- replicate(B, {
  g2 <- sample(g)
  abs(mean(x[g2 == 2]) - mean(x[g2 == 1]))
})
(mean(perm >= obs) + 1) / (B + 1)   # 双侧置换 p（含观测值）

动手练习：置换检验实现

使用 Melanoma 数据集，比较两组肿瘤厚度的置换检验：

library(MASS)
data(Melanoma)

# 1. 提取数据
x <- Melanoma$thickness
g <- Melanoma$ulcer  # 0=无溃疡, 1=有溃疡

# 2. 计算观测统计量（两组均值差）
obs_diff <- mean(x[g == 1]) - mean(x[g == 0])
cat("观测均值差：", round(obs_diff, 3), "mm\n\n")

# 3. 置换检验
set.seed(123)
B <- 9999  # 置换次数
perm_diffs <- numeric(B)

for (i in 1:B) {
  # 随机置换组别标签
  g_perm <- sample(g)
  # 计算置换后的统计量
  perm_diffs[i] <- mean(x[g_perm == 1]) - mean(x[g_perm == 0])
}

# 4. 计算双侧 p 值
p_value <- (sum(abs(perm_diffs) >= abs(obs_diff)) + 1) / (B + 1)
cat("置换检验 p 值：", round(p_value, 4), "\n\n")

# 5. 与传统 t 检验比较
t_test_p <- t.test(x ~ g)$p.value
cat("Welch t 检验 p 值：", round(t_test_p, 4), "\n")

# 6. 绘制零分布
cat("\n零分布分位数：\n")
cat("2.5%:", round(quantile(perm_diffs, 0.025), 3), "\n")
cat("50%:", round(quantile(perm_diffs, 0.5), 3), "\n")
cat("97.5%:", round(quantile(perm_diffs, 0.975), 3), "\n")

输出结果：

思考题：

• 置换检验和传统 t 检验有什么区别？什么情况下应该使用置换检验？
• 为什么置换检验不需要正态性假设？它的核心假设是什么？

问题：两组厚度分布是否近似正态？方差是否齐？若违背假设，你会选择什么检验方法比较两组均值？

问题 1：为什么方差要用 n-1 而不是 n 作分母？

问题 2：如果数据中存在一个极端大值（如 50 mm），均值和标准差会如何变化？中位数呢？

问题：95% 置信区间的正确解读是什么？能否说"真实参数有 95% 的概率落在这个区间内"？

问题：在两组比较时，何时选用 Welch t 检验、配对 t 检验、Mann-Whitney U 检验或置换检验？

问题：在基因表达分析中，若对 20,000 个基因分别做 t 检验，不校正会有什么后果？应选用 Bonferroni 还是 FDR 校正？

问题：使用倾向得分匹配（PSM）估计因果效应时，需要哪些关键假设？若存在未测混杂，会有什么后果？