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

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

library(MASS)
data(Melanoma)

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

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

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

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

动手练习：结果解读

方法选择建议：由于正态性和方差齐性均被违背，建议使用 Welch t 检验 或 Mann-Whitney U 检验

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

任务：手动计算肿瘤厚度的样本均值、方差和标准差

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

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

# 与 R 内置函数对比
mean(x); var(x); sd(x)

动手练习：计算结果

思考：为什么方差要用 n-1 而不是 n 作分母？若存在极端大值，均值和标准差会如何变化？中位数呢？

动手练习：计算置信区间

计算肿瘤厚度总体均值的 95% 置信区间，并比较不同置信水平下区间的宽度

library(MASS)
data(Melanoma)
# 计算描述统计
x <- Melanoma$thickness
cat("样本量：", length(x), "\n")
cat("样本均值：", round(mean(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")
# 比较不同置信水平
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[2] - ci90[1], 3), "\n")
cat("99% CI 宽度：", round(ci99[2] - ci99[1], 3), "\n")

思考：置信水平越高，区间越宽还是越窄？为什么？

动手练习：置信区间结果

关键结论：置信水平越高，区间越宽。因为更高的置信度需要更大的"容错空间"来确保捕获真实参数。

动手练习：两组比较检验

任务：比较不同溃疡组的肿瘤厚度，练习选择合适的检验方法

library(MASS)
data(Melanoma)
# 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）
library(effectsize)
cohens_d(thickness ~ ulcer, data = Melanoma)

动手练习：两组比较结果

方法选择：正态+方差齐→标准 t 检验；正态+方差不齐→Welch t 检验；非正态→Mann-Whitney U

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

任务：分析肿瘤厚度与年龄之间的关系

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)

动手练习：相关与回归结果

解读：年龄与厚度存在弱正相关，但 R² 仅 0.02，说明年龄只能解释厚度变异的 2%，需考虑其他因素。

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

任务：对分类变量进行描述统计分析

library(MASS)
data(Melanoma)
# 1. 性别分布（二分类变量）
table(Melanoma$sex)  # 频数
prop.table(table(Melanoma$sex)) * 100  # 构成比（%）
# 2. 溃疡状态分布
table(Melanoma$ulcer)
prop.table(table(Melanoma$ulcer)) * 100
# 3. 生存状态（多分类）: status: 1=死于黑色素瘤, 2=存活, 3=死于其他原因
table(Melanoma$status)
prop.table(table(Melanoma$status)) * 100
# 4. 创建有序变量：按厚度分组
Melanoma$thickness_grp <- cut(Melanoma$thickness,breaks = c(0, 2, 5, Inf),
                          labels = c("薄", "中", "厚"),right = FALSE)
table(Melanoma$thickness_grp)

动手练习：分类变量结果

注意：分类变量分析前需明确变量类型（二分类/有序/无序），不同类型采用不同统计方法。

动手练习：率的置信区间

任务：计算死亡率和溃疡发生率的置信区间

library(MASS)
data(Melanoma)

# 1. 计算黑色素瘤相关死亡率（status = 1）
n_death <- sum(Melanoma$status == 1)
n_total <- nrow(Melanoma)
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)

动手练习：率的 CI 结果

方法选择：Wilson 方法是率置信区间的推荐默认方法；小样本或 p 接近 0/1 时考虑精确方法。

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

任务：检验溃疡状态与性别之间是否存在关联

library(MASS)
data(Melanoma)

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

# 2. 计算期望频数并检查条件
chisq_test <- chisq.test(tb)
print("期望频数：")
print(round(chisq_test$expected, 1))
print(paste("最小期望频数：", round(min(chisq_test$expected), 2)))

# 3. 卡方检验
print(chisq_test)

# 4. Fisher 精确检验（对比）
fisher_test <- fisher.test(tb)
print(fisher_test)

动手练习：卡方检验结果

条件检查：所有期望频数 ≥ 5，满足卡方检验条件；若期望频数 < 5 的格子超过 20%，应使用 Fisher 精确检验。

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

任务：分析溃疡状态与死亡风险的关系

library(MASS)
data(Melanoma)

# 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(tb)

# 3. 手动计算 OR = (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("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)
summary(fit)

动手练习：OR 与 Logistic 结果

注意：OR ≠ RR，当结局发生率较高时（>10%），OR 会高估 RR。

动手练习：多重检验校正

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

# 模拟 20 个基因的 p 值（假设其中 3 个有真实效应）
set.seed(123)
p_values <- c(runif(17, 0, 1), runif(3, 0, 0.01))
p_sorted <- sort(p_values)

# 1. 不校正：直接用 0.05 阈值
cat("不校正（α = 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")

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

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

动手练习：多重检验校正结果

FDR 控制：BH 方法控制假发现率 FDR = E[假阳性数/总发现数]，允许一定比例的假阳性以换取更高的发现率。

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

任务：演示训练/测试划分及其重要性

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)

# 5. 在训练集和测试集上分别评估
train_pred <- predict(fit, type = "response")
test_pred <- predict(fit, newdata = test_data, type = "response")

动手练习：训练/测试结果

关键原则：永远不能用训练数据评估模型最终性能，必须使用独立的测试数据。

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

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

library(MASS)
data(Melanoma)

# 1. 定义"处理"为溃疡状态，结局为死亡
Melanoma$treat <- Melanoma$ulcer
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. 按倾向得分分组查看处理分配
by(Melanoma$ps, Melanoma$treat, summary)

# 5. 检查共同支撑（重叠）
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")

动手练习：PSM 结果与假设

注意：未调整的比较（粗差异 0.18）可能有偏，PSM 旨在模拟随机化，但仍依赖无未测混杂假设。

动手练习：置换检验实现

任务：比较两组肿瘤厚度的置换检验

library(MASS)
data(Melanoma)

# 1. 提取数据
x <- Melanoma$thickness
g <- Melanoma$ulcer

# 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")

动手练习：置换检验结果

置换检验优势：不依赖正态性假设，通过随机置换组别标签构建零分布，适用于小样本或分布未知的情况。

谢谢

返回演示列表 · 进入 04 讲义