在Rstudio中安装本节课所有需要的包，运行以下指令：

install.packages(c(
  "dplyr", "tidyr", "data.table",
  "ggplot2", "patchwork", "RColorBrewer", "factoextra",
  "cluster", "mclust", "dbscan", "Rtsne", "umap",
  "arules", "arulesViz", "arulesSequences",
  "caret", "randomForest", "xgboost", "iml", "DALEX", "pdp", "lime",
  "forecast", "prophet", "fpp2",
  "tm", "SnowballC", "wordcloud", "topicmodels", "tidytext", "jiebaR",
  "recommenderlab",
  "isotree", "dbscan",
  "purrr", "broom", "knitr", "mlbench"
))

章节	主题	核心内容
第九章	无监督学习	聚类方法、降维、异常检测
第十章	关联规则与模式挖掘	Apriori、FP-Growth、频繁项集
第十一章	特征选择与降维	过滤法、包装法、嵌入法
第十二章	模型解释与可解释性	SHAP、LIME、PDP、公平性
第十三章	时间序列机器学习	特征工程、Prophet、交叉验证
第十四章	文本挖掘基础	TF-IDF、主题模型、文本分类
第十五章	推荐系统基础	协同过滤、矩阵分解、评估指标

第九章无监督学习

9.1 聚类方法概述

什么是聚类？

聚类（Clustering） 是一种无监督学习方法，目标是将相似的样本归为同一组（簇），使组内样本相似度高、组间样本相似度低。

与分类的区别：

对比维度	分类	聚类
学习类型	监督学习	无监督学习
标签	有标签	无标签
目标	预测已知类别	发现数据结构
评估	准确率等	轮廓系数等

聚类方法分类

方法类型	代表算法	核心思想	适用场景
划分方法	K-Means	最小化簇内距离	球形簇、规模适中
层次方法	凝聚式聚类	自底向上合并	层次结构、小数据
密度方法	DBSCAN	密度相连区域	任意形状、有噪声
模型方法	GMM	概率分布拟合	概率聚类、软分配
图方法	谱聚类	图分割	非凸形状、复杂结构

意义与使用场景

核心意义：聚类帮助发现数据的内在结构，是数据探索和模式发现的基础工具。

举例：患者分型研究

为什么这种数据适合用聚类？

在医学研究中，往往存在多种疾病亚型，但临床诊断标准可能无法完全区分。例如，糖尿病患者可能有不同的代谢特征组合，传统分类方法可能遗漏某些亚型。聚类分析可以基于患者的多项指标（血糖、胰岛素、BMI、血脂等）自动发现潜在的患者亚群。

研究目的是什么？

发现疾病亚型：识别具有相似特征的患者群体
精准医疗：为不同亚型制定个性化治疗方案
风险分层：识别高风险或特殊预后的患者群体

可能得到什么样的结果？

聚类可能发现3-4个患者亚群：代谢正常型、肥胖相关型、自身免疫型等。每个亚群可能有不同的治疗响应和预后。

需要注意什么？

聚类结果需要结合临床知识进行解释和验证。不同的聚类方法和参数可能得到不同的结果，需要进行稳定性分析。

9.2 K-Means聚类

基本原理

K-Means 是最常用的聚类算法，通过迭代优化将数据划分为K个簇：

算法步骤： 1. 随机初始化K个聚类中心 2. 将每个样本分配到最近的聚类中心 3. 重新计算每个簇的中心（均值） 4. 重复步骤2-3直到收敛

目标函数：

\[J = \sum_{k=1}^{K} \sum_{x_i \in C_k} \|x_i - \mu_k\|^2\]

其中\(\mu_k\)是簇\(C_k\)的中心。

# K-Means聚类示例
set.seed(42)

# 生成模拟数据：三个患者群体
n_per_group <- 50

# 创建三种不同特征的患者群体
patient_data <- data.frame(
  # 特征1：炎症指标
  inflammation = c(
    rnorm(n_per_group, mean = 5, sd = 1),    # 低炎症组
    rnorm(n_per_group, mean = 10, sd = 1.5), # 中炎症组
    rnorm(n_per_group, mean = 15, sd = 2)    # 高炎症组
  ),
  # 特征2：代谢指标
  metabolism = c(
    rnorm(n_per_group, mean = 3, sd = 0.8),
    rnorm(n_per_group, mean = 6, sd = 1),
    rnorm(n_per_group, mean = 9, sd = 1.2)
  ),
  # 真实标签（仅用于验证，聚类时不用）
  true_group = rep(c("低风险", "中风险", "高风险"), each = n_per_group)
)

# 执行K-Means聚类（假设不知道真实分组）
kmeans_result <- kmeans(patient_data[, c("inflammation", "metabolism")], 
                        centers = 3, nstart = 25)

# 将聚类结果添加到数据框
patient_data$cluster <- as.factor(kmeans_result$cluster)

# 查看聚类中心
print("聚类中心：")

## [1] "聚类中心："

print(kmeans_result$centers)

##   inflammation metabolism
## 1    14.936596   9.109263
## 2     4.975035   3.027561
## 3    10.374221   6.118794

# 查看簇大小
print("各簇样本数：")

## [1] "各簇样本数："

print(kmeans_result$size)

## [1] 46 51 53

可视化聚类结果

# 可视化聚类结果
p1 <- ggplot2::ggplot(patient_data, ggplot2::aes(x = inflammation, y = metabolism)) +
  ggplot2::geom_point(ggplot2::aes(color = cluster), size = 3, alpha = 0.7) +
  ggplot2::geom_point(data = as.data.frame(kmeans_result$centers), 
                      ggplot2::aes(x = inflammation, y = metabolism),
                      shape = 4, size = 5, stroke = 2, color = "black") +
  ggplot2::labs(
    title = "K-Means聚类结果",
    subtitle = "X标记为聚类中心",
    x = "炎症指标", y = "代谢指标", color = "聚类"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

# 对比真实分组
p2 <- ggplot2::ggplot(patient_data, ggplot2::aes(x = inflammation, y = metabolism)) +
  ggplot2::geom_point(ggplot2::aes(color = true_group), size = 3, alpha = 0.7) +
  ggplot2::labs(
    title = "真实分组",
    x = "炎症指标", y = "代谢指标", color = "真实组别"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

p1 + p2 + patchwork::plot_annotation(
  title = "K-Means聚类 vs 真实分组对比"
)

肘部法则选择K值

肘部法则（Elbow Method）：绘制不同K值对应的簇内平方和（WSS），选择拐点作为最佳K值。

# 肘部法则：计算不同K值的WSS
wss_values <- sapply(1:10, function(k) {
  kmeans(patient_data[, c("inflammation", "metabolism")], 
         centers = k, nstart = 25)$tot.withinss
})

# 创建数据框用于可视化
elbow_df <- data.frame(
  k = 1:10,
  wss = wss_values
)

# 绘制肘部图
ggplot2::ggplot(elbow_df, ggplot2::aes(x = k, y = wss)) +
  ggplot2::geom_line(color = "steelblue", linewidth = 1) +
  ggplot2::geom_point(color = "steelblue", size = 3) +
  ggplot2::geom_vline(xintercept = 3, linetype = "dashed", color = "red") +
  ggplot2::annotate("text", x = 3.3, y = max(wss_values) * 0.9, 
                    label = "肘部点 K=3", color = "red") +
  ggplot2::labs(
    title = "肘部法则选择K值",
    subtitle = "WSS下降速度明显放缓的点为最佳K值",
    x = "聚类数 K", y = "簇内平方和（WSS）"
  ) +
  ggplot2::theme_minimal()

轮廓系数评估聚类质量

轮廓系数（Silhouette Coefficient） 衡量样本与同簇样本的相似度 vs 与其他簇样本的相似度：

\[s(i) = \frac{b(i) - a(i)}{\max(a(i), b(i))}\]

其中： - \(a(i)\)：样本\(i\)与同簇其他样本的平均距离 - \(b(i)\)：样本\(i\)与最近其他簇样本的平均距离

取值范围：\(-1 \leq s(i) \leq 1\)，越接近1表示聚类效果越好。

# 计算轮廓系数
silhouette_result <- cluster::silhouette(
  kmeans_result$cluster,
  dist(patient_data[, c("inflammation", "metabolism")])
)

# 使用factoextra包可视化轮廓系数
factoextra::fviz_silhouette(silhouette_result) +
  ggplot2::labs(
    title = "轮廓系数图",
    subtitle = "正值表示聚类合理，负值表示可能分错簇"
  ) +
  ggplot2::theme_minimal()

##   cluster size ave.sil.width
## 1       1   46          0.53
## 2       2   51          0.72
## 3       3   53          0.56

# 计算平均轮廓系数
avg_silhouette <- mean(silhouette_result[, 3])
print(paste("平均轮廓系数:", round(avg_silhouette, 3)))

## [1] "平均轮廓系数: 0.606"

初始化方法

K-Means对初始中心敏感，常用改进方法：

方法	说明	R实现
随机初始化	随机选择K个样本作为初始中心	`kmeans()`默认
K-Means++	选择距离已选中心较远的点	`kmeans(nstart > 1)`
多次运行	运行多次选择最优结果	`nstart`参数

# 使用nstart参数多次运行，选择最优结果
kmeans_multi <- kmeans(patient_data[, c("inflammation", "metabolism")], 
                       centers = 3, nstart = 50)

print(paste("最优结果的WSS:", round(kmeans_multi$tot.withinss, 2)))

## [1] "最优结果的WSS: 434.39"

意义与使用场景

核心意义：K-Means是最常用的聚类方法，简单高效，适合大规模数据的初步探索。

举例：医院患者分型

为什么这种数据适合用K-Means？

医院积累了大量患者的检验指标数据，如血常规、生化指标等。这些数据通常是连续型数值，且患者群体可能存在自然的分组（如不同疾病类型、不同严重程度）。K-Means可以快速发现这些潜在的患者亚群。

研究目的是什么？

患者分型：根据检验指标自动识别患者亚群
资源配置：不同亚群可能需要不同的医疗资源
辅助诊断：为新患者快速匹配相似病例

可能得到什么样的结果？

K-Means可能将患者分为3-5个亚群，每个亚群有特征性的检验指标组合。医生可以据此制定针对性的诊疗方案。

需要注意什么？

K-Means假设簇是球形的，对于非球形或大小差异大的簇效果不好。需要预先指定K值，建议结合肘部法则和轮廓系数确定。对异常值敏感，建议先进行异常值处理。

9.3 层次聚类

基本原理

层次聚类（Hierarchical Clustering） 通过构建聚类的层次结构（树状图）来组织数据，不需要预先指定聚类数。

两种主要方法：

方法	方向	过程
凝聚式（Agglomerative）	自底向上	每个样本开始为一个簇，逐步合并
分裂式（Divisive）	自顶向下	所有样本开始为一个簇，逐步分裂

凝聚式聚类

常用链接方法：

链接方法	簇间距离定义	特点
单链接	最近点距离	发现链状结构
全链接	最远点距离	发现紧凑簇
平均链接	平均距离	平衡的方法
Ward方法	簇内方差增量	发现球形簇

# 层次聚类示例
set.seed(42)

# 生成模拟数据
n_samples <- 30
hclust_data <- data.frame(
  feature1 = c(rnorm(n_samples/3, 2, 0.5), 
               rnorm(n_samples/3, 5, 0.5), 
               rnorm(n_samples/3, 8, 0.5)),
  feature2 = c(rnorm(n_samples/3, 2, 0.5), 
               rnorm(n_samples/3, 5, 0.5), 
               rnorm(n_samples/3, 2, 0.5))
)

# 计算距离矩阵
dist_matrix <- dist(hclust_data)

# 执行层次聚类（使用Ward方法）
hclust_result <- hclust(dist_matrix, method = "ward.D2")

# 绘制树状图
plot(hclust_result, main = "层次聚类树状图", 
     xlab = "样本", ylab = "距离", sub = "")

# 添加矩形框标记聚类
rect.hclust(hclust_result, k = 3, border = "red")

树状图解读

# 使用factoextra包绘制更美观的树状图
factoextra::fviz_dend(hclust_result, k = 3, 
                      cex = 0.8, 
                      k_colors = c("#2E9FDF", "#00AFBB", "#E7B800"),
                      rect = TRUE, 
                      rect_border = c("#2E9FDF", "#00AFBB", "#E7B800"),
                      rect_fill = TRUE,
                      lower_rect = -0.5) +
  ggplot2::labs(
    title = "层次聚类树状图",
    subtitle = "不同颜色表示不同聚类"
  )

切割树状图获取聚类

# 切割树状图，获取3个聚类
clusters <- cutree(hclust_result, k = 3)

# 将聚类结果添加到数据框
hclust_data$cluster <- as.factor(clusters)

# 可视化聚类结果
ggplot2::ggplot(hclust_data, ggplot2::aes(x = feature1, y = feature2)) +
  ggplot2::geom_point(ggplot2::aes(color = cluster), size = 4, alpha = 0.7) +
  ggplot2::labs(
    title = "层次聚类结果",
    subtitle = "通过切割树状图获得3个聚类",
    x = "特征1", y = "特征2", color = "聚类"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

意义与使用场景

核心意义：层次聚类提供了数据的层次结构视图，便于理解数据的组织方式。

举例：疾病分类研究

为什么这种数据适合用层次聚类？

疾病分类往往具有层次结构，如疾病大类下分亚类，亚类下再分亚型。层次聚类可以揭示这种自然的层次关系，帮助理解疾病的分类体系。

研究目的是什么？

发现疾病的层次分类结构
在不同粒度上理解疾病关系
为临床诊断提供参考框架

可能得到什么样的结果？

层次聚类可能揭示疾病从大类到亚类的自然分组，如自身免疫病→系统性→系统性红斑狼疮。

需要注意什么？

层次聚类一旦合并就不能撤销，可能导致早期错误传播。计算复杂度较高，不适合大规模数据。

9.4 DBSCAN密度聚类

基本原理

DBSCAN（Density-Based Spatial Clustering of Applications with Noise） 是一种基于密度的聚类算法，能够发现任意形状的簇，并识别噪声点。

核心概念：

概念	定义
核心点	ε邻域内至少有MinPts个点的点
边界点	不是核心点，但在某核心点的ε邻域内
噪声点	既不是核心点也不是边界点

参数： - ε（eps）：邻域半径 - MinPts：形成核心点的最小邻居数

# DBSCAN聚类示例
set.seed(42)

# 生成非球形数据：两个同心圆
n_per_circle <- 100

# 外圆
theta_outer <- runif(n_per_circle, 0, 2 * pi)
outer_circle <- data.frame(
  x = 5 * cos(theta_outer) + rnorm(n_per_circle, 0, 0.2),
  y = 5 * sin(theta_outer) + rnorm(n_per_circle, 0, 2),
  true_cluster = "外圆"
)

# 内圆
theta_inner <- runif(n_per_circle, 0, 2 * pi)
inner_circle <- data.frame(
  x = 2 * cos(theta_inner) + rnorm(n_per_circle, 0, 0.2),
  y = 2 * sin(theta_inner) + rnorm(n_per_circle, 0, 0.2),
  true_cluster = "内圆"
)

# 合并数据
dbscan_data <- rbind(outer_circle, inner_circle)

# 执行DBSCAN聚类
dbscan_result <- dbscan::dbscan(dbscan_data[, c("x", "y")], 
                                 eps = 0.5, minPts = 5)

# 查看聚类结果
print("DBSCAN聚类结果：")

## [1] "DBSCAN聚类结果："

print(table(dbscan_result$cluster))

## 
##  0  1  2  3  4 
## 88  8 20 60 24

# 0表示噪声点
print(paste("噪声点数量:", sum(dbscan_result$cluster == 0)))

## [1] "噪声点数量: 88"

可视化DBSCAN结果

# 将聚类结果添加到数据框
dbscan_data$cluster <- as.factor(dbscan_result$cluster)

# 可视化
ggplot2::ggplot(dbscan_data, ggplot2::aes(x = x, y = y)) +
  ggplot2::geom_point(ggplot2::aes(color = cluster), size = 2, alpha = 0.7) +
  ggplot2::scale_color_manual(
    values = c("0" = "gray50", "1" = "steelblue", "2" = "salmon"),
    labels = c("0" = "噪声", "1" = "簇1", "2" = "簇2")
  ) +
  ggplot2::labs(
    title = "DBSCAN聚类结果",
    subtitle = "能够发现任意形状的簇，灰色点为噪声",
    x = "X", y = "Y", color = "聚类"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

K-Means vs DBSCAN对比

# K-Means对同一数据的聚类结果
kmeans_for_compare <- kmeans(dbscan_data[, c("x", "y")], centers = 2, nstart = 25)
dbscan_data$kmeans_cluster <- as.factor(kmeans_for_compare$cluster)

# 对比可视化
p_dbscan <- ggplot2::ggplot(dbscan_data, ggplot2::aes(x = x, y = y)) +
  ggplot2::geom_point(ggplot2::aes(color = cluster), size = 2, alpha = 0.7) +
  ggplot2::labs(title = "DBSCAN", x = "", y = "", color = "聚类") +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "none")

p_kmeans <- ggplot2::ggplot(dbscan_data, ggplot2::aes(x = x, y = y)) +
  ggplot2::geom_point(ggplot2::aes(color = kmeans_cluster), size = 2, alpha = 0.7) +
  ggplot2::labs(title = "K-Means", x = "", y = "", color = "聚类") +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "none")

p_dbscan + p_kmeans + patchwork::plot_annotation(
  title = "DBSCAN vs K-Means：非球形数据聚类对比",
  subtitle = "DBSCAN能正确识别同心圆结构，K-Means无法处理"
)

参数选择

# 使用k-距离图选择eps参数
# 计算每个点到其第k近邻的距离
k <- 5  # minPts - 1
knn_dist <- dbscan::kNNdist(dbscan_data[, c("x", "y")], k = k)

# kNNdist当k为单个值时返回向量，直接使用
# 绘制k-距离图（排序后）
knn_dist_sorted <- sort(knn_dist, decreasing = FALSE)

plot(knn_dist_sorted, type = "l", 
     main = "K-距离图（用于选择eps）",
     xlab = "样本（按距离排序）", ylab = paste0(k, "-NN距离"))
abline(h = 0.5, col = "red", lty = 2)
text(length(knn_dist_sorted) * 0.8, 0.6, "eps = 0.5", col = "red")

意义与使用场景

核心意义：DBSCAN能够发现任意形状的簇，并自动识别噪声点，适合复杂数据结构。

举例：医学影像异常区域检测

为什么这种数据适合用DBSCAN？

医学影像中的病变区域往往形状不规则，且正常组织和异常组织的密度分布不同。DBSCAN可以基于像素特征的密度分布，发现任意形状的异常区域，同时识别噪声点（如伪影）。

研究目的是什么？

发现影像中的异常区域
自动识别噪声和伪影
辅助医生定位病变

可能得到什么样的结果？

DBSCAN可能识别出多个密度较高的区域（疑似病变），同时标记出孤立的噪声点。医生可以重点检查这些区域。

需要注意什么？

参数eps和minPts的选择对结果影响很大。对于密度不均匀的数据，可能需要使用改进算法如HDBSCAN。

9.5 高斯混合模型（GMM）与EM算法

基本原理

高斯混合模型（Gaussian Mixture Model, GMM） 假设数据由多个高斯分布混合生成：

\[p(x) = \sum_{k=1}^{K} \pi_k \mathcal{N}(x | \mu_k, \Sigma_k)\]

其中： - \(\pi_k\)：第k个高斯分量的混合系数（权重） - \(\mu_k\)：第k个高斯分量的均值 - \(\Sigma_k\)：第k个高斯分量的协方差矩阵

EM算法： 1. E步：计算每个样本属于各高斯分量的概率（后验概率） 2. M步：根据后验概率更新参数（均值、协方差、权重） 3. 重复直到收敛

# GMM聚类示例
set.seed(42)

# 生成模拟数据：两个重叠的高斯分布
n_per_cluster <- 100

gmm_data <- data.frame(
  x = c(rnorm(n_per_cluster, mean = 0, sd = 1),
        rnorm(n_per_cluster, mean = 3, sd = 1.5)),
  y = c(rnorm(n_per_cluster, mean = 0, sd = 1),
        rnorm(n_per_cluster, mean = 2, sd = 1))
)

# 使用mclust包进行GMM聚类
gmm_result <- mclust::Mclust(gmm_data, G = 2)  # 指定2个分量

# 查看模型摘要
print(gmm_result)

## 'Mclust' model object: (EVI,2) 
## 
## Available components: 
##  [1] "call"           "data"           "modelName"      "n"             
##  [5] "d"              "G"              "BIC"            "loglik"        
##  [9] "df"             "bic"            "icl"            "hypvol"        
## [13] "parameters"     "z"              "classification" "uncertainty"

可视化GMM结果

# 使用factoextra包可视化
factoextra::fviz_cluster(gmm_result, data = gmm_data,
                         ellipse.type = "confidence",
                         ellipse.level = 0.95,
                         geom = "point",
                         palette = "jco") +
  ggplot2::labs(
    title = "高斯混合模型聚类",
    subtitle = "椭圆表示95%置信区域"
  ) +
  ggplot2::theme_minimal()

软聚类 vs 硬聚类

# GMM提供软聚类（概率分配）
# 查看每个样本属于各簇的概率
prob_df <- as.data.frame(gmm_result$z)
colnames(prob_df) <- c("P(簇1)", "P(簇2)")
prob_df$sample <- 1:nrow(prob_df)
prob_df$assigned_cluster <- gmm_result$classification

# 查看前10个样本的概率分配
head(prob_df, 10)

##       P(簇1)       P(簇2) sample assigned_cluster
## 1  0.9999325 6.753952e-05      1                1
## 2  0.9944366 5.563433e-03      2                1
## 3  0.8401161 1.598839e-01      3                1
## 4  0.3338216 6.661784e-01      4                2
## 5  0.9998823 1.177108e-04      5                1
## 6  0.9998853 1.146762e-04      6                1
## 7  0.9918100 8.190009e-03      7                1
## 8  0.9998381 1.619277e-04      8                1
## 9  0.9728473 2.715268e-02      9                1
## 10 0.9973525 2.647517e-03     10                1

# 可视化概率分布
prob_df_melt <- tidyr::pivot_longer(prob_df, 
                                     cols = c("P(簇1)", "P(簇2)"),
                                     names_to = "Cluster", 
                                     values_to = "Probability")

ggplot2::ggplot(prob_df_melt, ggplot2::aes(x = sample, y = Probability, fill = Cluster)) +
  ggplot2::geom_bar(stat = "identity", position = "stack", width = 1) +
  ggplot2::labs(
    title = "GMM软聚类：样本属于各簇的概率",
    subtitle = "颜色深浅表示概率大小",
    x = "样本", y = "概率"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

意义与使用场景

核心意义：GMM提供概率性的聚类结果，适合需要不确定性估计的场景。

举例：疾病风险分层

为什么这种数据适合用GMM？

疾病风险往往是连续的，患者可能同时具有多种风险特征。GMM的软聚类可以给出患者属于各风险组的概率，更符合医学实际。

研究目的是什么？

风险分层：将患者分配到不同风险组
不确定性估计：提供概率而非硬标签
识别过渡状态：发现介于各风险组之间的患者

可能得到什么样的结果？

GMM可能发现患者以80%概率属于高风险组，20%概率属于中风险组，而非简单的”高风险”标签。

需要注意什么？

GMM假设数据服从高斯分布，对于非高斯数据效果可能不好。需要指定分量数量，可以使用BIC等准则自动选择。

9.6 谱聚类

基本原理

谱聚类（Spectral Clustering） 基于图论，将聚类问题转化为图分割问题：

算法步骤： 1. 构建相似度图（如k近邻图） 2. 计算图的拉普拉斯矩阵 3. 计算拉普拉斯矩阵的前k个特征向量 4. 对特征向量进行K-Means聚类

适用场景：非凸形状的簇、复杂结构数据

# 谱聚类示例
set.seed(42)

# 生成双月形数据
moon_data <- mlbench::mlbench.circle(n = 200, d = 2)
moon_df <- data.frame(
  x = moon_data$x[, 1],
  y = moon_data$x[, 2],
  true_class = moon_data$classes
)

# 使用kernlab包的谱聚类
library(kernlab)
spectral_result <- kernlab::specc(as.matrix(moon_df[, c("x", "y")]), 
                                   centers = 2)

# 添加聚类结果
moon_df$spectral_cluster <- as.factor(spectral_result@.Data)

# 可视化
ggplot2::ggplot(moon_df, ggplot2::aes(x = x, y = y)) +
  ggplot2::geom_point(ggplot2::aes(color = spectral_cluster), size = 2, alpha = 0.7) +
  ggplot2::labs(
    title = "谱聚类结果",
    subtitle = "能够正确识别非凸形状的簇",
    x = "X", y = "Y", color = "聚类"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

意义与使用场景

核心意义：谱聚类能够处理复杂非线性结构，是处理非凸簇的有效方法。

举例：基因表达谱聚类

为什么这种数据适合用谱聚类？

基因表达数据往往具有复杂的非线性结构，基因之间的调控关系可能形成非凸的簇。谱聚类能够捕捉这些复杂结构。

研究目的是什么？

发现共表达基因模块
识别功能相关的基因群
理解基因调控网络结构

可能得到什么样的结果？

谱聚类可能发现具有相似表达模式的基因群，这些基因可能参与相同的生物学过程。

需要注意什么？

谱聚类计算复杂度较高，对大规模数据需要使用近似方法。相似度图的构建方式对结果影响较大。

9.7 聚类评估指标

内部评估指标

不需要真实标签，仅基于数据本身评估聚类质量。

指标	公式/原理	取值范围	最佳值
轮廓系数	\(\frac{b-a}{\max(a,b)}\)	[-1, 1]	接近1
Davies-Bouldin指数	簇内散度/簇间距离	[0, ∞)	接近0
Calinski-Harabasz指数	簇间方差/簇内方差	[0, ∞)	越大越好

# 使用factoextra包计算多种评估指标
# 准备数据
set.seed(42)
eval_data <- data.frame(
  x = c(rnorm(50, 0, 1), rnorm(50, 4, 1), rnorm(50, 8, 1)),
  y = c(rnorm(50, 0, 1), rnorm(50, 4, 1), rnorm(50, 0, 1))
)

# 执行K-Means
km_eval <- kmeans(eval_data, centers = 3, nstart = 25)

# 轮廓系数
silhouette_eval <- cluster::silhouette(km_eval$cluster, dist(eval_data))
avg_sil <- mean(silhouette_eval[, 3])

# Davies-Bouldin指数
db_index <- clusterSim::index.DB(eval_data, km_eval$cluster)$DB

# Calinski-Harabasz指数
ch_index <- fpc::calinhara(eval_data, km_eval$cluster)

# 汇总结果
eval_results <- data.frame(
  指标 = c("轮廓系数", "Davies-Bouldin指数", "Calinski-Harabasz指数"),
  值 = c(avg_sil, db_index, ch_index),
  最佳值 = c("接近1", "接近0", "越大越好")
)

print(eval_results)

##                    指标          值   最佳值
## 1              轮廓系数   0.6778823    接近1
## 2    Davies-Bouldin指数   0.4907522    接近0
## 3 Calinski-Harabasz指数 532.0862203 越大越好

外部评估指标

需要真实标签，评估聚类结果与真实分类的一致性。

指标	说明	取值范围
兰德指数（Rand Index）	样本对分类一致性	[0, 1]
调整兰德指数（ARI）	校正随机情况的RI	[-1, 1]
互信息（MI）	聚类与标签的互信息	[0, ∞)
归一化互信息（NMI）	归一化的MI	[0, 1]

# 假设有真实标签
set.seed(42)
true_labels <- rep(c(1, 2, 3), each = 50)

# 计算调整兰德指数
ari <- mclust::adjustedRandIndex(km_eval$cluster, true_labels)
print(paste("调整兰德指数:", round(ari, 3)))

## [1] "调整兰德指数: 1"

# 计算归一化互信息
nmi <- aricode::NMI(km_eval$cluster, true_labels)
print(paste("归一化互信息:", round(nmi, 3)))

## [1] "归一化互信息: 1"

意义与使用场景

核心意义：聚类评估指标帮助选择最佳聚类方法和参数，确保聚类结果有意义。

举例：确定最佳聚类数

为什么需要评估指标？

在实际应用中，真实的聚类数往往未知。通过比较不同K值下的评估指标，可以选择最优的聚类数。

研究目的是什么？

确定最佳聚类数
比较不同聚类方法
验证聚类结果的稳定性

可能得到什么样的结果？

轮廓系数在K=3时达到最大，表明数据最适合分为3个簇。

需要注意什么？

不同指标可能给出不同的最佳聚类数，需要综合考虑。内部指标不能替代领域知识的验证。

9.8 主成分分析（PCA）降维

基本原理

主成分分析（PCA） 是一种线性降维方法，通过正交变换将高维数据投影到低维空间，保留尽可能多的方差。

目标：找到数据方差最大的方向（主成分）

\[\max_w \text{Var}(Xw) = \max_w w^T \Sigma w\]

其中\(\Sigma\)是协方差矩阵。

# PCA示例：基因表达数据降维
set.seed(42)

# 生成模拟基因表达数据（100个样本，50个基因）
n_samples <- 100
n_genes <- 50

# 创建具有分组结构的基因表达数据
gene_expr <- matrix(rnorm(n_samples * n_genes), nrow = n_samples, ncol = n_genes)

# 添加分组效应
gene_expr[1:50, 1:10] <- gene_expr[1:50, 1:10] + 2  # 第一组样本的前10个基因高表达
gene_expr[51:100, 11:20] <- gene_expr[51:100, 11:20] + 2  # 第二组样本的中间基因高表达

# 转换为数据框
gene_df <- as.data.frame(gene_expr)
colnames(gene_df) <- paste0("Gene_", 1:n_genes)

# 执行PCA
pca_result <- prcomp(gene_df, scale. = TRUE, center = TRUE)

# 查看主成分的方差贡献
var_explained <- pca_result$sdev^2 / sum(pca_result$sdev^2)
var_df <- data.frame(
  PC = paste0("PC", 1:10),
  Variance = var_explained[1:10] * 100,
  Cumulative = cumsum(var_explained[1:10]) * 100
)

print("前10个主成分的方差贡献：")

## [1] "前10个主成分的方差贡献："

print(head(var_df, 10))

##      PC  Variance Cumulative
## 1   PC1 21.864839   21.86484
## 2   PC2  4.922180   26.78702
## 3   PC3  4.356587   31.14361
## 4   PC4  4.213363   35.35697
## 5   PC5  3.772029   39.12900
## 6   PC6  3.487178   42.61618
## 7   PC7  3.416547   46.03272
## 8   PC8  3.143081   49.17580
## 9   PC9  3.090488   52.26629
## 10 PC10  2.841843   55.10813

可视化PCA结果

# 碎石图：方差贡献
p1 <- ggplot2::ggplot(var_df, ggplot2::aes(x = PC, y = Variance)) +
  ggplot2::geom_bar(stat = "identity", fill = "steelblue") +
  ggplot2::geom_line(ggplot2::aes(group = 1), color = "red") +
  ggplot2::labs(
    title = "碎石图：各主成分方差贡献",
    x = "主成分", y = "方差贡献（%）"
  ) +
  ggplot2::theme_minimal()

# 累积方差贡献
p2 <- ggplot2::ggplot(var_df, ggplot2::aes(x = PC, y = Cumulative)) +
  ggplot2::geom_line(color = "steelblue", linewidth = 1) +
  ggplot2::geom_point(color = "steelblue", size = 2) +
  ggplot2::geom_hline(yintercept = 80, linetype = "dashed", color = "red") +
  ggplot2::labs(
    title = "累积方差贡献",
    x = "主成分", y = "累积方差贡献（%）"
  ) +
  ggplot2::theme_minimal()

p1 + p2 + patchwork::plot_annotation(
  title = "PCA方差分析"
)

PCA双标图

# 添加分组信息
gene_df$group <- factor(rep(c("Group_A", "Group_B"), each = 50))

# 创建PCA得分数据框
pca_scores <- as.data.frame(pca_result$x[, 1:2])
pca_scores$group <- gene_df$group

# 绘制PCA得分图
ggplot2::ggplot(pca_scores, ggplot2::aes(x = PC1, y = PC2)) +
  ggplot2::geom_point(ggplot2::aes(color = group), size = 3, alpha = 0.7) +
  ggplot2::labs(
    title = "PCA得分图：样本在前两个主成分上的分布",
    subtitle = paste0("PC1解释", round(var_explained[1] * 100, 1), 
                      "%方差，PC2解释", round(var_explained[2] * 100, 1), "%方差"),
    x = paste0("PC1 (", round(var_explained[1] * 100, 1), "%)"),
    y = paste0("PC2 (", round(var_explained[2] * 100, 1), "%)"),
    color = "分组"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

意义与使用场景

核心意义：PCA是最常用的降维方法，能够减少数据维度、去除噪声、可视化高维数据。

举例：基因表达谱可视化

为什么这种数据适合用PCA？

基因表达数据通常有数千个基因（高维），但样本数相对较少。PCA可以将高维数据降维到2-3维，便于可视化和发现样本分组模式。

研究目的是什么？

数据可视化：将高维数据投影到2D/3D空间
降噪：去除方差较小的成分（通常是噪声）
特征提取：用主成分作为新特征

可能得到什么样的结果？

PCA可能揭示样本在PC1-PC2平面上的自然分组，不同疾病类型的患者可能形成明显的簇。

需要注意什么？

PCA是线性方法，对于非线性结构可能效果不好。主成分的解释性可能较差，需要结合载荷矩阵分析。

9.9 流形学习：t-SNE与UMAP

t-SNE原理

t-SNE（t-Distributed Stochastic Neighbor Embedding） 是一种非线性降维方法，特别适合高维数据的可视化。

核心思想： 1. 在高维空间计算样本间的相似度（高斯分布） 2. 在低维空间计算样本间的相似度（t分布） 3. 最小化两个相似度分布的KL散度

# t-SNE示例
set.seed(42)

# 使用Rtsne包进行t-SNE降维
library(Rtsne)

# 对基因表达数据进行t-SNE
tsne_result <- Rtsne::Rtsne(gene_df[, 1:n_genes], 
                             dims = 2, 
                             perplexity = 30,
                             verbose = FALSE)

# 创建结果数据框
tsne_df <- as.data.frame(tsne_result$Y)
colnames(tsne_df) <- c("tSNE1", "tSNE2")
tsne_df$group <- gene_df$group

# 可视化t-SNE结果
ggplot2::ggplot(tsne_df, ggplot2::aes(x = tSNE1, y = tSNE2)) +
  ggplot2::geom_point(ggplot2::aes(color = group), size = 3, alpha = 0.7) +
  ggplot2::labs(
    title = "t-SNE降维可视化",
    subtitle = "非线性降维，能够保持局部结构",
    x = "t-SNE维度1", y = "t-SNE维度2", color = "分组"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

UMAP原理

UMAP（Uniform Manifold Approximation and Projection） 是另一种流行的非线性降维方法，比t-SNE更快，且能更好地保持全局结构。

# UMAP示例
library(umap)

# 执行UMAP降维
umap_result <- umap::umap(gene_df[, 1:n_genes])

# 创建结果数据框
umap_df <- as.data.frame(umap_result$layout)
colnames(umap_df) <- c("UMAP1", "UMAP2")
umap_df$group <- gene_df$group

# 可视化UMAP结果
ggplot2::ggplot(umap_df, ggplot2::aes(x = UMAP1, y = UMAP2)) +
  ggplot2::geom_point(ggplot2::aes(color = group), size = 3, alpha = 0.7) +
  ggplot2::labs(
    title = "UMAP降维可视化",
    subtitle = "比t-SNE更快，保持全局结构更好",
    x = "UMAP维度1", y = "UMAP维度2", color = "分组"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

PCA vs t-SNE vs UMAP对比

# 三种方法的对比可视化
p_pca <- ggplot2::ggplot(pca_scores, ggplot2::aes(x = PC1, y = PC2)) +
  ggplot2::geom_point(ggplot2::aes(color = group), size = 2, alpha = 0.7) +
  ggplot2::labs(title = "PCA", x = "", y = "") +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "none")

p_tsne <- ggplot2::ggplot(tsne_df, ggplot2::aes(x = tSNE1, y = tSNE2)) +
  ggplot2::geom_point(ggplot2::aes(color = group), size = 2, alpha = 0.7) +
  ggplot2::labs(title = "t-SNE", x = "", y = "") +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "none")

p_umap <- ggplot2::ggplot(umap_df, ggplot2::aes(x = UMAP1, y = UMAP2)) +
  ggplot2::geom_point(ggplot2::aes(color = group), size = 2, alpha = 0.7) +
  ggplot2::labs(title = "UMAP", x = "", y = "") +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "none")

p_pca + p_tsne + p_umap + patchwork::plot_annotation(
  title = "降维方法对比：PCA vs t-SNE vs UMAP"
)

意义与使用场景

核心意义：非线性降维方法能够揭示高维数据的复杂结构，是数据探索和可视化的有力工具。

举例：单细胞RNA测序数据分析

为什么这种数据适合用t-SNE/UMAP？

单细胞RNA测序数据通常有数千个基因，但每个细胞只表达其中一部分。数据具有高度非线性的结构，传统线性方法难以捕捉。t-SNE和UMAP能够揭示细胞亚群和分化轨迹。

研究目的是什么？

细胞类型识别：发现不同的细胞亚群
分化轨迹推断：理解细胞发育过程
异质性分析：发现稀有细胞类型

可能得到什么样的结果？

t-SNE/UMAP图可能显示多个明显的细胞簇，对应不同的细胞类型或状态。

需要注意什么？

t-SNE和UMAP的结果受参数影响较大，不同运行可能得到不同结果。距离和簇大小在图中可能不反映真实情况，需要结合其他方法验证。

9.10 自编码器降维概念

基本原理

自编码器（Autoencoder） 是一种神经网络架构，用于学习数据的有效压缩表示：

结构： - 编码器：将输入压缩到低维潜在空间 - 解码器：从潜在表示重建输入

目标：最小化重建误差 \(\|x - \hat{x}\|^2\)

输入层 → 编码层 → 潜在层（低维） → 解码层 → 输出层

与传统降维方法的比较

方法	类型	优点	缺点
PCA	线性	快速、可解释	只能捕捉线性结构
t-SNE	非线性	可视化效果好	计算慢、不可逆
UMAP	非线性	快速、保持全局结构	参数敏感
自编码器	非线性	可学习复杂特征	需要大量数据、黑箱

医学应用场景

自编码器在医学领域的应用包括：

医学图像压缩：压缩MRI、CT图像，减少存储空间
异常检测：重建误差大的样本可能是异常
特征学习：学习医学图像的深层特征用于分类

注意：自编码器涉及深度学习，本教程仅介绍概念，具体实现可参考深度学习框架如Keras、PyTorch。

9.11 异常检测方法

异常检测概述

异常检测（Anomaly Detection） 识别与大多数数据显著不同的样本。在医学领域，异常可能对应疾病、误诊、数据错误等。

孤立森林（Isolation Forest）

原理：异常点更容易被孤立（需要更少的分割即可分离）

# 孤立森林示例
set.seed(42)

# 生成正常数据和异常数据
normal_data <- data.frame(
  feature1 = rnorm(200, mean = 0, sd = 1),
  feature2 = rnorm(200, mean = 0, sd = 1),
  type = "正常"
)

# 添加异常点
anomaly_data <- data.frame(
  feature1 = c(rnorm(10, mean = 5, sd = 0.5), rnorm(5, mean = -4, sd = 0.5)),
  feature2 = c(rnorm(10, mean = 5, sd = 0.5), rnorm(5, mean = -4, sd = 0.5)),
  type = "异常"
)

# 合并数据
all_data <- rbind(normal_data, anomaly_data)

# 使用isotree包实现孤立森林
library(isotree)
iso_forest <- isotree::isolation.forest(all_data[, c("feature1", "feature2")])

# 预测异常分数
all_data$anomaly_score <- predict(iso_forest, all_data[, c("feature1", "feature2")])

# 可视化异常分数
ggplot2::ggplot(all_data, ggplot2::aes(x = feature1, y = feature2)) +
  ggplot2::geom_point(ggplot2::aes(color = anomaly_score), size = 2) +
  ggplot2::scale_color_gradient(low = "blue", high = "red") +
  ggplot2::labs(
    title = "孤立森林异常检测",
    subtitle = "红色点表示高异常分数",
    x = "特征1", y = "特征2", color = "异常分数"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

局部异常因子（LOF）

原理：比较样本的局部密度与邻居的局部密度

# LOF异常检测
library(dbscan)

# 计算LOF分数
lof_scores <- dbscan::lof(all_data[, c("feature1", "feature2")], k = 5)

# 添加LOF分数到数据框
all_data$lof_score <- lof_scores

# 可视化LOF分数
ggplot2::ggplot(all_data, ggplot2::aes(x = feature1, y = feature2)) +
  ggplot2::geom_point(ggplot2::aes(color = lof_score), size = 2) +
  ggplot2::scale_color_gradient(low = "blue", high = "red") +
  ggplot2::labs(
    title = "LOF异常检测",
    subtitle = "LOF > 1表示密度低于邻居（可能是异常）",
    x = "特征1", y = "特征2", color = "LOF分数"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

单类SVM

原理：学习正常数据的边界，边界外的点为异常

# 单类SVM示例
library(e1071)

# 训练单类SVM（只使用正常数据）
ocsvm_model <- e1071::svm(
  x = as.matrix(normal_data[, c("feature1", "feature2")]),
  type = "one-classification",
  nu = 0.05,  # 预期异常比例
  kernel = "radial"
)

# 预测
all_data$ocsvm_pred <- predict(ocsvm_model, all_data[, c("feature1", "feature2")])

# 可视化
ggplot2::ggplot(all_data, ggplot2::aes(x = feature1, y = feature2)) +
  ggplot2::geom_point(ggplot2::aes(color = ocsvm_pred), size = 2) +
  ggplot2::scale_color_manual(values = c("FALSE" = "red", "TRUE" = "blue"),
                               labels = c("FALSE" = "异常", "TRUE" = "正常")) +
  ggplot2::labs(
    title = "单类SVM异常检测",
    x = "特征1", y = "特征2", color = "预测结果"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

意义与使用场景

核心意义：异常检测能够自动识别异常样本，在医学领域有重要应用。

举例：医疗数据质量控制

为什么需要异常检测？

医疗数据可能存在录入错误、设备故障、极端病例等异常情况。异常检测可以自动识别这些情况，提高数据质量。

研究目的是什么？

数据清洗：识别和修正数据录入错误
设备监控：检测医疗设备异常
疑难病例发现：识别需要特别关注的病例

可能得到什么样的结果？

异常检测可能发现某些患者的检验指标异常偏离正常范围，需要进一步核实。

需要注意什么？

异常不一定是错误，可能是真实的极端情况。需要结合领域知识判断异常的性质。

第十章关联规则与模式挖掘

10.1 关联规则基本概念

什么是关联规则？

关联规则（Association Rules） 发现数据项之间的有趣关联关系，最著名的例子是”啤酒与尿布”的购物篮分析。

基本概念：

概念	定义	说明
项集	项的集合	{A, B, C}
支持度	项集出现的频率	P(A∩B)
置信度	条件概率	P(B
提升度	规则的有效性	P(B

三个核心指标

支持度（Support）：

\[\text{Support}(A \rightarrow B) = \frac{|A \cap B|}{|D|}\]

表示A和B同时出现的频率。

置信度（Confidence）：

\[\text{Confidence}(A \rightarrow B) = \frac{\text{Support}(A \cup B)}{\text{Support}(A)}\]

表示在A出现的条件下B出现的概率。

提升度（Lift）：

\[\text{Lift}(A \rightarrow B) = \frac{\text{Confidence}(A \rightarrow B)}{\text{Support}(B)}\]

Lift > 1：A和B正相关
Lift = 1：A和B独立
Lift < 1：A和B负相关

# 关联规则示例：医学共病分析
# 创建模拟的疾病共病数据

set.seed(42)

# 假设有1000个患者的诊断记录
n_patients <- 1000

# 创建疾病诊断矩阵（患者 × 疾病）
diseases <- c("高血压", "糖尿病", "冠心病", "高血脂", "肥胖", "肾病")

# 生成模拟数据：某些疾病倾向于同时出现
diagnosis_data <- matrix(0, nrow = n_patients, ncol = length(diseases))
colnames(diagnosis_data) <- diseases

# 模拟疾病关联
for (i in 1:n_patients) {
  # 基础概率
  diagnosis_data[i, "高血压"] <- rbinom(1, 1, 0.3)
  diagnosis_data[i, "糖尿病"] <- rbinom(1, 1, 0.2)
  
  # 高血压患者更容易有冠心病
  if (diagnosis_data[i, "高血压"] == 1) {
    diagnosis_data[i, "冠心病"] <- rbinom(1, 1, 0.4)
  } else {
    diagnosis_data[i, "冠心病"] <- rbinom(1, 1, 0.1)
  }
  
  # 糖尿病患者更容易有肾病
  if (diagnosis_data[i, "糖尿病"] == 1) {
    diagnosis_data[i, "肾病"] <- rbinom(1, 1, 0.3)
  } else {
    diagnosis_data[i, "肾病"] <- rbinom(1, 1, 0.05)
  }
  
  # 其他疾病
  diagnosis_data[i, "高血脂"] <- rbinom(1, 1, 0.25)
  diagnosis_data[i, "肥胖"] <- rbinom(1, 1, 0.2)
}

# 转换为数据框
diagnosis_df <- as.data.frame(diagnosis_data)

# 查看数据结构
head(diagnosis_df)

##   高血压 糖尿病 冠心病 高血脂 肥胖 肾病
## 1      1      1      0      0    0    1
## 2      1      0      1      0    0    0
## 3      1      0      0      1    0    0
## 4      0      0      1      1    1    0
## 5      0      0      0      0    1    0
## 6      1      1      0      0    1    0

# 计算各疾病的患病率（支持度）
disease_support <- colMeans(diagnosis_df)
print("各疾病患病率：")

## [1] "各疾病患病率："

print(round(disease_support, 3))

## 高血压 糖尿病 冠心病 高血脂   肥胖   肾病 
##  0.301  0.227  0.183  0.256  0.214  0.122

意义与使用场景

核心意义：关联规则帮助发现数据中的隐含模式，是数据挖掘的基础方法。

举例：疾病共病模式发现

为什么这种数据适合用关联规则？

疾病之间存在复杂的共病关系，某些疾病可能倾向于同时出现。关联规则可以发现这些隐含的共病模式，帮助医生理解疾病之间的关系。

研究目的是什么？

发现共病模式：识别经常同时出现的疾病组合
风险预警：发现某种疾病后，预警可能的并发症
病理机制探索：为疾病关联机制研究提供线索

可能得到什么样的结果？

可能发现”高血压→冠心病”的强关联规则，提示高血压患者需要重点关注心血管健康。

需要注意什么？

关联关系不等于因果关系。发现的规则需要结合医学知识进行解释和验证。

10.2 Apriori算法

算法原理

Apriori算法 是最经典的关联规则挖掘算法，基于频繁项集的先验性质：

核心思想：如果某个项集是频繁的，那么它的所有子集也是频繁的。

算法步骤： 1. 扫描数据库，找出频繁1-项集（L1） 2. 由L1生成候选2-项集，找出频繁2-项集（L2） 3. 重复直到无法生成更大的频繁项集 4. 由频繁项集生成关联规则

# 使用arules包进行关联规则挖掘
library(arules)

# 将数据转换为transactions格式
# 首先将数据转换为适合arules的格式

# 创建交易数据：每个患者的疾病列表
transactions_list <- apply(diagnosis_df, 1, function(x) {
  diseases[x == 1]
})

# 转换为transactions对象
trans <- as(transactions_list, "transactions")

# 查看事务摘要
summary(trans)

## transactions as itemMatrix in sparse format with
##  1000 rows (elements/itemsets/transactions) and
##  6 columns (items) and a density of 0.2171667 
## 
## most frequent items:
##  高血压  高血脂  糖尿病    肥胖  冠心病 (Other) 
##     301     256     227     214     183     122 
## 
## element (itemset/transaction) length distribution:
## sizes
##   0   1   2   3   4   5 
## 258 366 232 107  33   4 
## 
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   1.000   1.303   2.000   5.000 
## 
## includes extended item information - examples:
##   labels
## 1   肥胖
## 2 高血压
## 3 高血脂

# 执行Apriori算法
rules <- arules::apriori(trans, 
                          parameter = list(
                            support = 0.05,    # 最小支持度
                            confidence = 0.5,  # 最小置信度
                            minlen = 2,        # 规则最小长度
                            maxlen = 4         # 规则最大长度
                          ))

## Apriori
## 
## Parameter specification:
##  confidence minval smax arem  aval originalSupport maxtime support minlen
##         0.5    0.1    1 none FALSE            TRUE       5    0.05      2
##  maxlen target  ext
##       4  rules TRUE
## 
## Algorithmic control:
##  filter tree heap memopt load sort verbose
##     0.1 TRUE TRUE  FALSE TRUE    2    TRUE
## 
## Absolute minimum support count: 50 
## 
## set item appearances ...[0 item(s)] done [0.00s].
## set transactions ...[6 item(s), 1000 transaction(s)] done [0.00s].
## sorting and recoding items ... [6 item(s)] done [0.00s].
## creating transaction tree ... done [0.00s].
## checking subsets of size 1 2 3 done [0.00s].
## writing ... [2 rule(s)] done [0.00s].
## creating S4 object  ... done [0.00s].

# 查看规则摘要
summary(rules)

## set of 2 rules
## 
## rule length distribution (lhs + rhs):sizes
## 2 
## 2 
## 
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       2       2       2       2       2       2 
## 
## summary of quality measures:
##     support         confidence        coverage           lift      
##  Min.   :0.0810   Min.   :0.5847   Min.   :0.1220   Min.   :1.943  
##  1st Qu.:0.0875   1st Qu.:0.6045   1st Qu.:0.1373   1st Qu.:2.188  
##  Median :0.0940   Median :0.6243   Median :0.1525   Median :2.434  
##  Mean   :0.0940   Mean   :0.6243   Mean   :0.1525   Mean   :2.434  
##  3rd Qu.:0.1005   3rd Qu.:0.6441   3rd Qu.:0.1678   3rd Qu.:2.679  
##  Max.   :0.1070   Max.   :0.6639   Max.   :0.1830   Max.   :2.925  
##      count      
##  Min.   : 81.0  
##  1st Qu.: 87.5  
##  Median : 94.0  
##  Mean   : 94.0  
##  3rd Qu.:100.5  
##  Max.   :107.0  
## 
## mining info:
##   data ntransactions support confidence
##  trans          1000    0.05        0.5
##                                                                                                       call
##  arules::apriori(data = trans, parameter = list(support = 0.05, confidence = 0.5, minlen = 2, maxlen = 4))

查看和筛选规则

# 查看前10条规则
arules::inspect(head(rules, 10))

##     lhs         rhs      support confidence coverage lift     count
## [1] {肾病}   => {糖尿病} 0.081   0.6639344  0.122    2.924821  81  
## [2] {冠心病} => {高血压} 0.107   0.5846995  0.183    1.942523 107

# 按提升度排序
rules_by_lift <- sort(rules, by = "lift", decreasing = TRUE)

# 查看提升度最高的规则
print("提升度最高的10条规则：")

## [1] "提升度最高的10条规则："

arules::inspect(head(rules_by_lift, 10))

##     lhs         rhs      support confidence coverage lift     count
## [1] {肾病}   => {糖尿病} 0.081   0.6639344  0.122    2.924821  81  
## [2] {冠心病} => {高血压} 0.107   0.5846995  0.183    1.942523 107

# 筛选特定规则的子集
# 例如：找出所有涉及"糖尿病"的规则
diabetes_rules <- subset(rules, subset = items %in% "糖尿病")

print("涉及糖尿病的规则：")

## [1] "涉及糖尿病的规则："

arules::inspect(head(diabetes_rules, 5))

##     lhs       rhs      support confidence coverage lift     count
## [1] {肾病} => {糖尿病} 0.081   0.6639344  0.122    2.924821 81

可视化关联规则

library(arulesViz)

# 方法1：散点图
plot(rules, measure = c("support", "confidence"), shading = "lift",
     main = "关联规则散点图\n颜色深浅表示提升度")

# 方法2：网络图（显示规则结构）
plot(head(rules_by_lift, 20), method = "graph",
     main = "关联规则网络图（前20条）")

## Available control parameters (with default values):
## layout    =  stress
## circular  =  FALSE
## ggraphdots    =  NULL
## edges     =  <environment>
## nodes     =  <environment>
## nodetext  =  <environment>
## colors    =  c("#EE0000FF", "#EEEEEEFF")
## engine    =  ggplot2
## max   =  100
## verbose   =  FALSE

# 方法3：平行坐标图
plot(head(rules_by_lift, 15), method = "paracoord",
     main = "关联规则平行坐标图")

意义与使用场景

核心意义：Apriori算法能够高效发现数据中的频繁模式和关联规则。

举例：药物相互作用发现

为什么这种数据适合用Apriori？

药物不良反应数据中，某些药物组合可能导致特定的不良反应。Apriori可以发现这些药物-不良反应的关联模式。

研究目的是什么？

药物安全监测：发现潜在的药物相互作用
不良反应预警：识别高风险药物组合
临床决策支持：为医生用药提供参考

可能得到什么样的结果？

可能发现”药物A + 药物B → 不良反应C”的规则，提示这两种药物联用需要谨慎。

需要注意什么？

Apriori需要多次扫描数据库，对大规模数据效率较低。最小支持度和置信度的设置对结果影响很大。

10.3 FP-Growth算法

算法原理

FP-Growth（Frequent Pattern Growth） 是Apriori的改进算法，不需要生成候选项集：

核心思想：使用FP树（Frequent Pattern Tree）压缩存储数据，直接从FP树挖掘频繁项集。

优点： 1. 只需扫描数据库两次 2. 不生成候选项集 3. 效率比Apriori高

缺点： 1. FP树可能占用大量内存 2. 实现较复杂

# FP-Growth在R中通过rCBA包实现
# 这里使用arules包的eclat函数作为替代（类似效果）

# 使用eclat算法挖掘频繁项集
frequent_itemsets <- arules::eclat(trans, 
                                    parameter = list(
                                      support = 0.05,
                                      minlen = 2,
                                      maxlen = 4
                                    ))

## Eclat
## 
## parameter specification:
##  tidLists support minlen maxlen            target  ext
##     FALSE    0.05      2      4 frequent itemsets TRUE
## 
## algorithmic control:
##  sparse sort verbose
##       7   -2    TRUE
## 
## Absolute minimum support count: 50 
## 
## create itemset ... 
## set transactions ...[6 item(s), 1000 transaction(s)] done [0.00s].
## sorting and recoding items ... [6 item(s)] done [0.00s].
## creating bit matrix ... [6 row(s), 1000 column(s)] done [0.00s].
## writing  ... [7 set(s)] done [0.00s].
## Creating S4 object  ... done [0.00s].

# 查看频繁项集
print("频繁项集：")

## [1] "频繁项集："

arules::inspect(head(sort(frequent_itemsets, by = "support"), 10))

##     items            support count
## [1] {高血压, 冠心病} 0.107   107  
## [2] {肾病, 糖尿病}   0.081    81  
## [3] {高血压, 糖尿病} 0.077    77  
## [4] {高血压, 高血脂} 0.076    76  
## [5] {肥胖, 高血脂}   0.056    56  
## [6] {肥胖, 高血压}   0.053    53  
## [7] {高血脂, 糖尿病} 0.050    50

# 从频繁项集生成规则
rules_from_fp <- arules::ruleInduction(frequent_itemsets, trans,
                                        confidence = 0.5)

# 查看生成的规则
print("从频繁项集生成的规则：")

## [1] "从频繁项集生成的规则："

arules::inspect(head(sort(rules_from_fp, by = "lift"), 10))

##     lhs         rhs      support confidence lift     itemset
## [1] {肾病}   => {糖尿病} 0.081   0.6639344  2.924821 1      
## [2] {冠心病} => {高血压} 0.107   0.5846995  1.942523 4

Apriori vs FP-Growth对比

对比维度	Apriori	FP-Growth
扫描次数	多次	2次
候选项集	需要生成	不需要
内存使用	较低	较高（FP树）
适用场景	稀疏数据	密集数据
实现复杂度	简单	复杂

10.4 频繁项集挖掘

频繁项集的类型

类型	定义	示例
频繁项集	支持度≥阈值的项集	{高血压, 糖尿病}
闭频繁项集	不存在超集有相同支持度	{高血压, 糖尿病}支持度=0.1，无超集支持度=0.1
最大频繁项集	不存在频繁超集	{高血压, 糖尿病, 冠心病}

# 挖掘不同类型的频繁项集

# 普通频繁项集
freq_items <- arules::eclat(trans, parameter = list(support = 0.05))

## Eclat
## 
## parameter specification:
##  tidLists support minlen maxlen            target  ext
##     FALSE    0.05      1     10 frequent itemsets TRUE
## 
## algorithmic control:
##  sparse sort verbose
##       7   -2    TRUE
## 
## Absolute minimum support count: 50 
## 
## create itemset ... 
## set transactions ...[6 item(s), 1000 transaction(s)] done [0.00s].
## sorting and recoding items ... [6 item(s)] done [0.00s].
## creating bit matrix ... [6 row(s), 1000 column(s)] done [0.00s].
## writing  ... [13 set(s)] done [0.00s].
## Creating S4 object  ... done [0.00s].

print(paste("普通频繁项集数量:", length(freq_items)))

## [1] "普通频繁项集数量: 13"

# 闭频繁项集
closed_items <- arules::eclat(trans, parameter = list(support = 0.05, 
                                                       target = "closed frequent itemsets"))

## Eclat
## 
## parameter specification:
##  tidLists support minlen maxlen                   target  ext
##     FALSE    0.05      1     10 closed frequent itemsets TRUE
## 
## algorithmic control:
##  sparse sort verbose
##       7   -2    TRUE
## 
## Absolute minimum support count: 50 
## 
## create itemset ... 
## set transactions ...[6 item(s), 1000 transaction(s)] done [0.00s].
## sorting and recoding items ... [6 item(s)] done [0.00s].
## creating bit matrix ... [6 row(s), 1000 column(s)] done [0.00s].
## writing  ... [13 set(s)] done [0.00s].
## Creating S4 object  ... done [0.00s].

print(paste("闭频繁项集数量:", length(closed_items)))

## [1] "闭频繁项集数量: 13"

# 最大频繁项集
maximal_items <- arules::eclat(trans, parameter = list(support = 0.05, 
                                                        target = "maximally frequent itemsets"))

## Eclat
## 
## parameter specification:
##  tidLists support minlen maxlen                      target  ext
##     FALSE    0.05      1     10 maximally frequent itemsets TRUE
## 
## algorithmic control:
##  sparse sort verbose
##       7   -2    TRUE
## 
## Absolute minimum support count: 50 
## 
## create itemset ... 
## set transactions ...[6 item(s), 1000 transaction(s)] done [0.00s].
## sorting and recoding items ... [6 item(s)] done [0.00s].
## creating bit matrix ... [6 row(s), 1000 column(s)] done [0.00s].
## writing  ... [7 set(s)] done [0.00s].
## Creating S4 object  ... done [0.00s].

print(paste("最大频繁项集数量:", length(maximal_items)))

## [1] "最大频繁项集数量: 7"

# 可视化对比
itemset_counts <- data.frame(
  类型 = c("普通频繁项集", "闭频繁项集", "最大频繁项集"),
  数量 = c(length(freq_items), length(closed_items), length(maximal_items))
)

ggplot2::ggplot(itemset_counts, ggplot2::aes(x = 类型, y = 数量)) +
  ggplot2::geom_bar(stat = "identity", fill = "steelblue") +
  ggplot2::geom_text(ggplot2::aes(label = 数量), vjust = -0.5) +
  ggplot2::labs(
    title = "不同类型频繁项集数量对比",
    subtitle = "闭频繁项集和最大频繁项集是压缩表示",
    x = "项集类型", y = "数量"
  ) +
  ggplot2::theme_minimal()

10.5 序列模式挖掘基础

什么是序列模式？

序列模式（Sequential Pattern） 发现事件发生的时间顺序规律：

示例： - 购物序列：{牛奶} → {面包} → {鸡蛋} - 疾病发展：{高血压} → {冠心病} → {心衰}

与关联规则的区别： - 关联规则：发现同时出现的关系 - 序列模式：发现先后出现的关系

# 序列模式挖掘示例
library(arulesSequences)

# 创建模拟的疾病发展序列数据
# 假设跟踪患者的诊断历史

set.seed(42)
n_patients_seq <- 100

# 创建序列数据框
seq_data <- data.frame(
  sequenceID = rep(1:n_patients_seq, each = 3),  # 患者ID
  eventID = rep(1:3, n_patients_seq),             # 时间点
  items = ""                                       # 诊断
)

# 生成模拟序列
for (i in 1:n_patients_seq) {
  idx <- ((i-1)*3 + 1):(i*3)
  
  # 时间点1：初始诊断
  initial <- sample(c("高血压", "糖尿病", "高血脂"), 1)
  seq_data$items[idx[1]] <- initial
  
  # 时间点2：进展（与初始诊断相关）
  if (initial == "高血压") {
    seq_data$items[idx[2]] <- sample(c("冠心病", "肾病", "无进展"), 1, 
                                      prob = c(0.4, 0.2, 0.4))
  } else if (initial == "糖尿病") {
    seq_data$items[idx[2]] <- sample(c("肾病", "视网膜病变", "无进展"), 1,
                                      prob = c(0.3, 0.2, 0.5))
  } else {
    seq_data$items[idx[2]] <- sample(c("冠心病", "无进展"), 1,
                                      prob = c(0.3, 0.7))
  }
  
  # 时间点3：进一步进展
  seq_data$items[idx[3]] <- sample(c("心衰", "脑卒中", "终末期肾病", "稳定"), 1,
                                    prob = c(0.2, 0.2, 0.1, 0.5))
}

# 转换为transactions格式
seq_trans <- as(seq_data, "transactions")

# 查看序列数据结构
print("序列数据示例：")

## [1] "序列数据示例："

head(seq_data, 9)

##   sequenceID eventID  items
## 1          1       1 高血压
## 2          1       2   肾病
## 3          1       3   稳定
## 4          2       1 高血压
## 5          2       2 无进展
## 6          2       3   心衰
## 7          3       1 糖尿病
## 8          3       2 无进展
## 9          3       3   心衰

序列模式挖掘

# 使用cSPADE算法挖掘序列模式
# 注意：需要正确格式的数据

# 简单示例：计算序列支持度
seq_support <- seq_data %>%
  dplyr::group_by(sequenceID) %>%
  dplyr::summarise(
    seq = paste(items, collapse = " → "),
    .groups = "drop"
  )

# 统计常见序列
seq_freq <- seq_support %>%
  dplyr::count(seq, sort = TRUE)

print("最常见的疾病发展序列：")

## [1] "最常见的疾病发展序列："

head(seq_freq, 10)

## # A tibble: 10 × 2
##    seq                          n
##    <chr>                    <int>
##  1 高血脂 → 无进展 → 稳定      13
##  2 糖尿病 → 无进展 → 稳定       8
##  3 高血压 → 冠心病 → 稳定       7
##  4 高血脂 → 无进展 → 脑卒中     6
##  5 糖尿病 → 肾病 → 稳定         5
##  6 高血压 → 肾病 → 稳定         5
##  7 糖尿病 → 无进展 → 心衰       4
##  8 糖尿病 → 肾病 → 脑卒中       4
##  9 高血压 → 无进展 → 稳定       4
## 10 高血脂 → 冠心病 → 稳定       4

# 可视化
ggplot2::ggplot(head(seq_freq, 10), ggplot2::aes(x = reorder(seq, n), y = n)) +
  ggplot2::geom_bar(stat = "identity", fill = "steelblue") +
  ggplot2::coord_flip() +
  ggplot2::labs(
    title = "最常见的疾病发展序列",
    x = "序列", y = "频次"
  ) +
  ggplot2::theme_minimal()

意义与使用场景

核心意义：序列模式挖掘发现事件的时间顺序规律，对疾病进展预测有重要价值。

举例：疾病进展轨迹分析

为什么这种数据适合用序列模式挖掘？

疾病发展往往遵循一定的时间顺序，如高血压→冠心病→心衰。序列模式挖掘可以发现这些进展轨迹。

研究目的是什么？

疾病进展建模：理解疾病的发展过程
早期干预：识别高风险进展路径
预后预测：预测患者可能的疾病发展

可能得到什么样的结果？

可能发现”高血压→冠心病→心衰”是一条常见的疾病进展路径，提示需要早期干预。

需要注意什么？

序列数据需要准确的时间戳记录。数据质量和随访完整性对结果影响很大。

第十一章特征选择与降维

11.1 特征选择概述

为什么需要特征选择？

特征选择（Feature Selection） 从原始特征中选择最相关的子集，目的是：

目的	说明
降低维度	减少计算复杂度
提高性能	去除噪声特征，提高模型性能
防止过拟合	减少特征数量，降低模型复杂度
提高可解释性	特征更少，模型更易理解

特征选择方法分类

方法类型	原理	优点	缺点
过滤法	基于统计指标筛选	快速、独立于模型	忽略特征间交互
包装法	基于模型性能筛选	考虑特征组合	计算量大
嵌入法	模型训练时自动选择	效率高	依赖特定模型

意义与使用场景

核心意义：特征选择是提高模型性能和可解释性的关键步骤。

举例：基因标志物筛选

为什么这种数据适合用特征选择？

基因表达数据通常有数千到数万个基因，但只有少数与疾病相关。特征选择可以从海量基因中筛选出真正有诊断价值的标志物。

研究目的是什么？

标志物发现：识别与疾病相关的关键基因
模型简化：减少特征数量，提高模型可解释性
成本降低：减少需要检测的基因数量

可能得到什么样的结果？

可能从10000个基因中筛选出20个关键基因，这些基因组合可以准确预测疾病状态。

需要注意什么？

特征选择可能导致信息丢失，需要在模型性能和特征数量之间权衡。

11.2 过滤法

方差阈值

原理：删除方差过小的特征（变化不大的特征信息量小）。

# 方差阈值示例
set.seed(42)

# 创建模拟数据：100个样本，20个特征
n <- 100
p <- 20

filter_data <- matrix(rnorm(n * p), nrow = n, ncol = p)
colnames(filter_data) <- paste0("Feature_", 1:p)

# 让某些特征方差很小
filter_data[, 1] <- rnorm(n, mean = 5, sd = 0.01)  # 低方差
filter_data[, 2] <- rnorm(n, mean = 5, sd = 0.02)  # 低方差

filter_df <- as.data.frame(filter_data)

# 计算各特征的方差
variances <- apply(filter_df, 2, var)

# 创建方差数据框
var_df <- data.frame(
  Feature = names(variances),
  Variance = variances
)

# 可视化方差分布
ggplot2::ggplot(var_df, ggplot2::aes(x = reorder(Feature, Variance), y = Variance)) +
  ggplot2::geom_bar(stat = "identity", fill = "steelblue") +
  ggplot2::geom_hline(yintercept = 0.1, color = "red", linetype = "dashed") +
  ggplot2::coord_flip() +
  ggplot2::labs(
    title = "各特征方差分布",
    subtitle = "红色虚线为阈值，低于阈值的特征将被删除",
    x = "特征", y = "方差"
  ) +
  ggplot2::theme_minimal()

# 选择方差大于阈值的特征
threshold <- 0.1
selected_features <- names(variances[variances > threshold])
print(paste("方差阈值筛选后保留的特征数:", length(selected_features)))

## [1] "方差阈值筛选后保留的特征数: 18"

互信息

原理：衡量特征与目标变量之间的信息共享量。

# 互信息示例
# 使用caret包计算特征重要性

set.seed(42)

# 创建分类数据
n <- 200
mi_data <- data.frame(
  feature1 = c(rnorm(n/2, 0, 1), rnorm(n/2, 2, 1)),
  feature2 = c(rnorm(n/2, 0, 1), rnorm(n/2, 2, 1)),
  feature3 = rnorm(n),  # 无信息
  feature4 = rnorm(n),  # 无信息
  target = factor(rep(c("ClassA", "ClassB"), each = n/2))
)

# 使用caret计算特征重要性
# 设置训练控制
ctrl <- caret::trainControl(method = "cv", number = 5)

# 训练简单模型获取重要性
model_imp <- caret::train(target ~ ., data = mi_data, 
                          method = "rf", 
                          trControl = ctrl,
                          importance = TRUE,
                          ntree = 50)

# 获取特征重要性
importance_scores <- caret::varImp(model_imp)

# 可视化
ggplot2::ggplot(importance_scores, ggplot2::aes(x = reorder(rownames(importance_scores), Overall), y = Overall)) +
  ggplot2::geom_bar(stat = "identity", fill = "steelblue") +
  ggplot2::coord_flip() +
  ggplot2::labs(
    title = "特征重要性（基于随机森林）",
    x = "特征", y = "重要性分数"
  ) +
  ggplot2::theme_minimal()

卡方检验

原理：检验类别特征与目标变量之间的独立性。

# 卡方检验示例
# 创建类别数据
set.seed(42)
n <- 200

chi_data <- data.frame(
  symptom1 = sample(c("有", "无"), n, replace = TRUE, prob = c(0.3, 0.7)),
  symptom2 = sample(c("有", "无"), n, replace = TRUE, prob = c(0.4, 0.6)),
  disease = sample(c("患病", "健康"), n, replace = TRUE, prob = c(0.2, 0.8))
)

# 让symptom1与疾病相关
chi_data$symptom1[chi_data$disease == "患病"] <- sample(c("有", "无"), 
                                                         sum(chi_data$disease == "患病"), 
                                                         replace = TRUE, 
                                                         prob = c(0.7, 0.3))

# 计算卡方检验
chi_test1 <- chisq.test(table(chi_data$symptom1, chi_data$disease))
chi_test2 <- chisq.test(table(chi_data$symptom2, chi_data$disease))

# 汇总结果
chi_results <- data.frame(
  Feature = c("symptom1", "symptom2"),
  ChiSquare = c(chi_test1$statistic, chi_test2$statistic),
  PValue = c(chi_test1$p.value, chi_test2$p.value)
)

print("卡方检验结果：")

## [1] "卡方检验结果："

print(chi_results)

##    Feature    ChiSquare       PValue
## 1 symptom1 23.004376576 1.616330e-06
## 2 symptom2  0.006823627 9.341654e-01

11.3 包装法

递归特征消除（RFE）

原理：反复构建模型，每次删除最不重要的特征。

# RFE示例
set.seed(42)

# 创建模拟数据
n <- 150
p <- 20

rfe_data <- matrix(rnorm(n * p), nrow = n, ncol = p)
colnames(rfe_data) <- paste0("Feature_", 1:p)

# 让前5个特征与目标相关
target <- rfe_data[, 1] + rfe_data[, 2] + rfe_data[, 3] + 
          rnorm(n, 0, 0.5)
target <- factor(ifelse(target > median(target), "High", "Low"))

rfe_df <- as.data.frame(rfe_data)
rfe_df$target <- target

# 使用caret的RFE
# 定义控制参数
rfe_ctrl <- caret::rfeControl(
  functions = caret::rfFuncs,  # 使用随机森林
  method = "cv",
  number = 5
)

# 执行RFE
rfe_result <- caret::rfe(
  x = rfe_df[, 1:p],
  y = rfe_df$target,
  sizes = c(1:10, 15, 20),
  rfeControl = rfe_ctrl
)

# 查看结果
print("RFE结果：")

## [1] "RFE结果："

print(rfe_result)

## 
## Recursive feature selection
## 
## Outer resampling method: Cross-Validated (5 fold) 
## 
## Resampling performance over subset size:
## 
##  Variables Accuracy  Kappa AccuracySD KappaSD Selected
##          1   0.6933 0.3867    0.07601 0.15202         
##          2   0.7133 0.4267    0.05055 0.10111         
##          3   0.8067 0.6133    0.04944 0.09888        *
##          4   0.7800 0.5600    0.07303 0.14606         
##          5   0.7733 0.5467    0.09545 0.19090         
##          6   0.7600 0.5200    0.06412 0.12824         
##          7   0.7733 0.5467    0.02789 0.05578         
##          8   0.7800 0.5600    0.05055 0.10111         
##          9   0.7600 0.5200    0.04346 0.08692         
##         10   0.7733 0.5467    0.04944 0.09888         
##         15   0.7867 0.5733    0.06912 0.13824         
##         20   0.7733 0.5467    0.10382 0.20763         
## 
## The top 3 variables (out of 3):
##    Feature_1, Feature_2, Feature_3

# 可视化
plot(rfe_result, main = "递归特征消除结果")

前向/后向特征选择

前向选择：从空集开始，逐步添加特征 后向消除：从全集开始，逐步删除特征

# 使用step函数进行逐步回归
set.seed(42)

# 创建线性回归数据
n <- 100
step_data <- data.frame(
  x1 = rnorm(n),
  x2 = rnorm(n),
  x3 = rnorm(n),  # 噪声
  x4 = rnorm(n),  # 噪声
  y = rnorm(n)
)

# 让y与x1, x2相关
step_data$y <- 2 * step_data$x1 - 1.5 * step_data$x2 + rnorm(n, 0, 0.5)

# 全模型
full_model <- lm(y ~ x1 + x2 + x3 + x4, data = step_data)

# 后向消除
backward_model <- step(full_model, direction = "backward", trace = 0)

# 查看最终模型
print("后向消除后的模型：")

## [1] "后向消除后的模型："

summary(backward_model)

## 
## Call:
## lm(formula = y ~ x1 + x2, data = step_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.93302 -0.36500  0.02268  0.34716  1.65764 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.008394   0.052881  -0.159    0.874    
## x1           1.981572   0.050792  39.013   <2e-16 ***
## x2          -1.616309   0.058498 -27.630   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.526 on 97 degrees of freedom
## Multiple R-squared:  0.9581, Adjusted R-squared:  0.9573 
## F-statistic:  1110 on 2 and 97 DF,  p-value: < 2.2e-16

# 前向选择
null_model <- lm(y ~ 1, data = step_data)
forward_model <- step(null_model, 
                       scope = list(lower = null_model, upper = full_model),
                       direction = "forward", trace = 0)

print("前向选择后的模型：")

## [1] "前向选择后的模型："

summary(forward_model)

## 
## Call:
## lm(formula = y ~ x1 + x2, data = step_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.93302 -0.36500  0.02268  0.34716  1.65764 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.008394   0.052881  -0.159    0.874    
## x1           1.981572   0.050792  39.013   <2e-16 ***
## x2          -1.616309   0.058498 -27.630   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.526 on 97 degrees of freedom
## Multiple R-squared:  0.9581, Adjusted R-squared:  0.9573 
## F-statistic:  1110 on 2 and 97 DF,  p-value: < 2.2e-16

11.4 嵌入法

L1正则化特征选择

原理：Lasso回归自动将不重要特征的系数压缩为0。

# Lasso特征选择示例
library(glmnet)

set.seed(42)

# 创建高维数据
n <- 100
p <- 50

lasso_data <- matrix(rnorm(n * p), nrow = n, ncol = p)

# 只有前10个特征与目标相关
true_coef <- c(rep(1, 10), rep(0, 40))
target <- lasso_data %*% true_coef + rnorm(n, 0, 0.5)

# 执行Lasso回归
lasso_model <- glmnet::glmnet(lasso_data, target, alpha = 1)  # alpha=1为Lasso

# 交叉验证选择lambda
cv_lasso <- glmnet::cv.glmnet(lasso_data, target, alpha = 1)

# 可视化
plot(cv_lasso, main = "Lasso交叉验证")

# 获取最优lambda下的系数
optimal_coef <- coef(cv_lasso, s = "lambda.min")

# 查看非零系数
non_zero_features <- which(optimal_coef != 0) - 1  # -1因为第一行是截距
print(paste("Lasso选择的特征数:", length(non_zero_features[non_zero_features > 0])))

## [1] "Lasso选择的特征数: 39"

print("选择的特征：")

## [1] "选择的特征："

print(non_zero_features[non_zero_features > 0])

##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 14 15 18 19 20 21 23 24 25 26 27 28 29
## [26] 30 31 32 33 36 38 39 41 42 43 44 45 49 50

基于树模型的特征重要性

# 随机森林特征重要性
set.seed(42)

# 创建数据
n <- 200
rf_data <- data.frame(
  feature1 = rnorm(n),
  feature2 = rnorm(n),
  feature3 = rnorm(n),
  feature4 = rnorm(n),
  feature5 = rnorm(n),
  target = factor(rep(c("A", "B"), each = n/2))
)

# 让某些特征与目标相关
rf_data$feature1[rf_data$target == "B"] <- rf_data$feature1[rf_data$target == "B"] + 1
rf_data$feature3[rf_data$target == "B"] <- rf_data$feature3[rf_data$target == "B"] - 1

# 训练随机森林
rf_model <- randomForest::randomForest(target ~ ., data = rf_data, importance = TRUE)

# 获取特征重要性
importance_df <- as.data.frame(randomForest::importance(rf_model))

# 可视化
ggplot2::ggplot(importance_df, ggplot2::aes(x = reorder(rownames(importance_df), MeanDecreaseGini), y = MeanDecreaseGini)) +
  ggplot2::geom_bar(stat = "identity", fill = "steelblue") +
  ggplot2::coord_flip() +
  ggplot2::labs(
    title = "随机森林特征重要性（Mean Decrease Gini）",
    x = "特征", y = "重要性"
  ) +
  ggplot2::theme_minimal()

11.5 特征重要性解释

排列重要性

原理：随机打乱某个特征的值，观察模型性能下降程度。

# 排列重要性示例
# 使用iml包

set.seed(42)

# 训练模型
model <- randomForest::randomForest(target ~ ., data = rf_data)

# 创建预测器对象
predictor <- iml::Predictor$new(model, data = rf_data[, -6], y = rf_data$target)

# 计算排列重要性
perm_importance <- iml::FeatureImp$new(predictor, loss = "ce")

# 可视化
plot(perm_importance) +
  ggplot2::labs(
    title = "排列重要性",
    subtitle = "打乱特征后模型性能下降越多，特征越重要"
  )

SHAP值

原理：基于博弈论，计算每个特征对预测的贡献。

# SHAP值示例
# 使用iml包

# 创建SHAP解释器
shap_explainer <- iml::Shapley$new(predictor, x.interest = rf_data[1, -6])

# 可视化单个样本的SHAP值
plot(shap_explainer) +
  ggplot2::labs(
    title = "SHAP值解释",
    subtitle = "每个特征对预测的贡献"
  )

意义与使用场景

核心意义：特征重要性解释帮助理解模型决策依据，提高模型可解释性。

举例：疾病风险预测模型解释

为什么需要特征重要性解释？

医生需要理解模型为什么预测某患者为高风险，才能信任并采纳模型建议。特征重要性可以告诉医生哪些因素导致了高风险预测。

研究目的是什么？

模型解释：理解模型的决策逻辑
特征验证：确认重要特征符合医学知识
发现新线索：可能发现新的风险因素

可能得到什么样的结果？

模型可能显示”年龄”和”血压”是最重要的预测因素，这与医学知识一致，增加了模型的可信度。

需要注意什么？

特征重要性不等于因果关系。相关性强的特征可能被模型高度依赖，但不一定是真正的病因。

11.6 PCA作为特征提取

特征选择 vs 特征提取

对比维度	特征选择	特征提取
方法	选择原始特征子集	创建新特征
可解释性	保持原始特征含义	新特征难以解释
信息保留	可能丢失信息	压缩信息
典型方法	过滤法、包装法	PCA、LDA

# PCA特征提取示例
set.seed(42)

# 创建高维数据
n <- 100
p <- 20

pca_feat_data <- matrix(rnorm(n * p), nrow = n, ncol = p)
colnames(pca_feat_data) <- paste0("Feature_", 1:p)

# 添加分组结构
pca_feat_data[1:50, 1:5] <- pca_feat_data[1:50, 1:5] + 2

# 执行PCA
pca_feat <- prcomp(pca_feat_data, scale. = TRUE)

# 使用前5个主成分作为新特征
pca_features <- pca_feat$x[, 1:5]

# 查看新特征的方差贡献
var_explained <- pca_feat$sdev^2 / sum(pca_feat$sdev^2)
print("前5个主成分的方差贡献：")

## [1] "前5个主成分的方差贡献："

print(round(cumsum(var_explained)[5], 3))

## [1] 0.466

11.7 LDA作为特征提取

LDA特征提取原理

线性判别分析（LDA） 寻找最能区分类别的投影方向，可用于特征提取。

# LDA特征提取示例
set.seed(42)

# 创建带标签的数据
n <- 100
lda_feat_data <- data.frame(
  x1 = c(rnorm(n/2, 0, 1), rnorm(n/2, 3, 1)),
  x2 = c(rnorm(n/2, 0, 1), rnorm(n/2, 3, 1)),
  x3 = rnorm(n),  # 噪声特征
  x4 = rnorm(n),  # 噪声特征
  class = factor(rep(c("A", "B"), each = n/2))
)

# 执行LDA
lda_feat <- MASS::lda(class ~ x1 + x2 + x3 + x4, data = lda_feat_data)

# 查看判别函数
print("LDA判别函数系数：")

## [1] "LDA判别函数系数："

print(lda_feat$scaling)

##            LD1
## x1 0.621599492
## x2 0.820935594
## x3 0.030917345
## x4 0.005188053

# 投影到判别方向
lda_projection <- as.matrix(lda_feat_data[, 1:4]) %*% lda_feat$scaling

# 可视化
lda_feat_data$LD1 <- lda_projection[, 1]

ggplot2::ggplot(lda_feat_data, ggplot2::aes(x = LD1, fill = class)) +
  ggplot2::geom_histogram(position = "identity", alpha = 0.5, bins = 30) +
  ggplot2::labs(
    title = "LDA特征提取结果",
    subtitle = "投影到判别方向后，两类分离更明显",
    x = "判别方向得分", y = "频数", fill = "类别"
  ) +
  ggplot2::theme_minimal()

意义与使用场景

核心意义：特征提取创建新的综合特征，可以在降低维度的同时保留关键信息。

举例：医学影像特征压缩

为什么这种数据适合用特征提取？

医学影像可能提取出数百个纹理、形状特征，其中很多高度相关。PCA可以将这些特征压缩为少数几个综合特征，同时保留大部分信息。

研究目的是什么？

降维：减少特征数量，降低计算复杂度
去相关：消除特征间的冗余
降噪：去除方差较小的成分

可能得到什么样的结果？

从100个影像特征中提取10个主成分，可能保留了90%的信息量，同时大大简化了模型。

需要注意什么？

提取的特征是原始特征的线性组合，可解释性较差。需要权衡降维效果和可解释性。

第十二章模型解释与可解释性

12.1 可解释性基本概念

为什么需要可解释性？

模型可解释性（Model Interpretability） 指人类能够理解模型决策原因的程度。在医学领域尤为重要：

原因	说明
信任建立	医生需要理解模型才能信任
伦理要求	医疗决策需要可解释
监管合规	医疗AI需要通过审核
错误诊断	发现模型的潜在问题

全局解释 vs 局部解释

类型	定义	目的
全局解释	解释模型整体行为	理解模型如何工作
局部解释	解释单个预测	理解特定预测的原因

可解释性方法分类

方法	适用模型	特点
特征重要性	所有模型	全局解释
PDP/ICE	所有模型	全局/局部解释
SHAP	所有模型	局部解释
LIME	所有模型	局部解释
替代模型	复杂模型	全局解释

意义与使用场景

核心意义：可解释性是医疗AI落地的关键，帮助建立信任和发现错误。

举例：疾病诊断模型解释

为什么需要可解释性？

医生不会仅凭一个”黑箱”模型的预测就做出诊断决策。他们需要知道模型为什么做出这个预测，哪些因素是关键，这有助于医生判断模型预测是否合理。

研究目的是什么？

建立信任：让医生理解并信任模型
发现错误：识别模型的潜在偏见
辅助决策：提供决策依据

可能得到什么样的结果？

模型可能显示”肿瘤大小”和”边界清晰度”是诊断恶性的主要因素，这与医学知识一致。

需要注意什么？

可解释性方法本身也可能有偏见，需要多种方法交叉验证。

12.2 特征重要性（Feature Importance）

基本原理

特征重要性 衡量每个特征对模型预测的贡献程度。

常用方法： 1. 基于模型的内置重要性（如随机森林的Gini重要性） 2. 排列重要性（Permutation Importance）

# 特征重要性示例
set.seed(42)

# 创建模拟疾病预测数据
n <- 300
disease_data <- data.frame(
  age = round(rnorm(n, 50, 15)),
  blood_pressure = rnorm(n, 120, 20),
  cholesterol = rnorm(n, 200, 40),
  bmi = rnorm(n, 25, 5),
  glucose = rnorm(n, 100, 20),
  smoking = factor(sample(c("Yes", "No"), n, replace = TRUE)),
  disease = factor(sample(c("Disease", "Healthy"), n, replace = TRUE))
)

# 让某些特征与疾病相关
disease_data$disease <- factor(ifelse(
  disease_data$blood_pressure > 130 & disease_data$glucose > 110,
  "Disease", "Healthy"
))

# 训练随机森林模型
rf_model <- randomForest::randomForest(disease ~ ., data = disease_data, 
                                        importance = TRUE, ntree = 100)

# 获取特征重要性
importance_matrix <- randomForest::importance(rf_model)

# 创建数据框用于可视化
imp_df <- data.frame(
  Feature = rownames(importance_matrix),
  MeanDecreaseAccuracy = importance_matrix[, "MeanDecreaseAccuracy"],
  MeanDecreaseGini = importance_matrix[, "MeanDecreaseGini"]
)

# 可视化
p1 <- ggplot2::ggplot(imp_df, ggplot2::aes(x = reorder(Feature, MeanDecreaseAccuracy), y = MeanDecreaseAccuracy)) +
  ggplot2::geom_bar(stat = "identity", fill = "steelblue") +
  ggplot2::coord_flip() +
  ggplot2::labs(title = "Mean Decrease Accuracy", x = "", y = "") +
  ggplot2::theme_minimal()

p2 <- ggplot2::ggplot(imp_df, ggplot2::aes(x = reorder(Feature, MeanDecreaseGini), y = MeanDecreaseGini)) +
  ggplot2::geom_bar(stat = "identity", fill = "salmon") +
  ggplot2::coord_flip() +
  ggplot2::labs(title = "Mean Decrease Gini", x = "", y = "") +
  ggplot2::theme_minimal()

p1 + p2 + patchwork::plot_annotation(
  title = "随机森林特征重要性"
)

排列重要性

# 使用DALEX包计算排列重要性
library(DALEX)

# 创建解释器
explainer_rf <- DALEX::explain(
  model = rf_model,
  data = disease_data[, -7],
  y = as.numeric(disease_data$disease == "Disease"),
  label = "Random Forest",
  verbose = FALSE
)

# 计算排列重要性
set.seed(42)
vip_rf <- DALEX::model_parts(explainer_rf, 
                              loss_function = DALEX::loss_root_mean_square)

# 可视化
plot(vip_rf) +
  ggplot2::labs(
    title = "排列重要性（DALEX）",
    subtitle = "打乱特征后模型性能下降越多，特征越重要"
  )

12.3 部分依赖图（PDP）与个体条件期望（ICE）

部分依赖图原理

部分依赖图（Partial Dependence Plot, PDP） 展示特征对预测的平均影响：

\[\text{PD}(x_s) = \frac{1}{n} \sum_{i=1}^{n} f(x_s, x_{c}^{(i)})\]

其中\(x_s\)是感兴趣的特征，\(x_c\)是其他特征。

# PDP示例
library(pdp)

# 计算血压的部分依赖
pdp_bp <- pdp::partial(rf_model, 
                        pred.var = "blood_pressure",
                        train = disease_data[, -7],
                        prob = TRUE)

# 可视化
ggplot2::ggplot(pdp_bp, ggplot2::aes(x = blood_pressure, y = yhat)) +
  ggplot2::geom_line(color = "steelblue", linewidth = 1) +
  ggplot2::labs(
    title = "部分依赖图：血压对疾病风险的影响",
    x = "血压", y = "预测概率"
  ) +
  ggplot2::theme_minimal()

个体条件期望（ICE）

ICE图 展示每个样本的预测曲线，比PDP提供更详细的信息。

# ICE图示例
ice_bp <- pdp::partial(rf_model, 
                        pred.var = "blood_pressure",
                        train = disease_data[, -7],
                        ice = TRUE,
                        prob = TRUE)

# 可视化ICE曲线
ggplot2::ggplot(ice_bp, ggplot2::aes(x = blood_pressure, y = yhat, group = yhat.id)) +
  ggplot2::geom_line(alpha = 0.1, color = "steelblue") +
  ggplot2::labs(
    title = "ICE图：每个样本的血压-风险曲线",
    subtitle = "每条线代表一个样本",
    x = "血压", y = "预测概率"
  ) +
  ggplot2::theme_minimal()

意义与使用场景

核心意义：PDP和ICE展示特征与预测的关系，帮助理解模型行为。

举例：药物剂量-效应关系分析

为什么这种数据适合用PDP/ICE？

药物剂量与疗效的关系可能非线性，PDP可以展示这种关系，帮助确定最佳剂量范围。

研究目的是什么？

理解特征影响：展示特征如何影响预测
发现非线性：识别非线性关系
确定阈值：找到关键决策点

可能得到什么样的结果？

PDP可能显示血压超过130后，疾病风险急剧上升，提示这是关键阈值。

需要注意什么？

PDP假设特征独立，如果特征相关，结果可能有误导。

12.4 SHAP值

SHAP原理

SHAP（SHapley Additive exPlanations） 基于博弈论，计算每个特征对预测的边际贡献：

核心思想：考虑所有可能的特征组合，计算特征加入时的平均贡献。

\[\phi_i = \sum_{S \subseteq F \setminus \{i\}} \frac{|S|!(|F|-|S|-1)!}{|F|!} [f(S \cup \{i\}) - f(S)]\]

# SHAP值示例
# 使用iml包

set.seed(42)

# 创建预测器
predictor <- iml::Predictor$new(
  rf_model, 
  data = disease_data[, -7], 
  y = disease_data$disease,
  class = "Disease"
)

# 计算单个样本的SHAP值
shap_single <- iml::Shapley$new(predictor, x.interest = disease_data[1, -7])

# 可视化
plot(shap_single) +
  ggplot2::labs(
    title = "SHAP值：单个样本的特征贡献",
    subtitle = "正值增加疾病概率，负值降低疾病概率"
  )

多样本SHAP分析

# 计算多个样本的SHAP值
shap_values <- data.frame()

for (i in sample(1:nrow(disease_data), 20)) {
  shap_i <- iml::Shapley$new(predictor, x.interest = disease_data[i, -7])
  shap_df <- shap_i$results
  shap_df$sample_id <- i
  shap_values <- rbind(shap_values, shap_df)
}

# 可视化SHAP值分布
ggplot2::ggplot(shap_values, ggplot2::aes(x = feature, y = phi, fill = phi > 0)) +
  ggplot2::geom_boxplot() +
  ggplot2::scale_fill_manual(values = c("TRUE" = "red", "FALSE" = "blue")) +
  ggplot2::coord_flip() +
  ggplot2::labs(
    title = "SHAP值分布",
    subtitle = "正值增加疾病概率，负值降低疾病概率",
    x = "特征", y = "SHAP值"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "none")

意义与使用场景

核心意义：SHAP提供最全面的特征贡献解释，是当前最流行的可解释性方法。

举例：个体化治疗方案解释

为什么这种数据适合用SHAP？

每个患者的疾病风险因素组合不同，SHAP可以为每个患者提供个性化的风险因素解释。

研究目的是什么？

个性化解释：为每个患者提供定制化解释
特征贡献量化：精确计算每个因素的贡献
干预建议：识别可干预的风险因素

可能得到什么样的结果？

对于某患者，SHAP可能显示”高血压贡献+15%，高血糖贡献+10%，年龄贡献+5%“。

需要注意什么？

SHAP计算量大，对大规模数据可能较慢。结果依赖于模型质量。

12.5 LIME

LIME原理

LIME（Local Interpretable Model-agnostic Explanations） 通过在局部拟合简单模型来解释复杂模型：

核心思想：在待解释样本附近生成扰动数据，用简单模型拟合局部行为。

# LIME示例
library(lime)
library(caret)

# 为LIME准备数据（仅使用数值特征）
lime_data <- disease_data[, c("age", "blood_pressure", "cholesterol", "bmi", "glucose")]
lime_labels <- disease_data$disease

# 使用caret训练模型（lime对caret模型有更好的支持）
caret_model <- caret::train(
  x = lime_data,
  y = lime_labels,
  method = "rf",
  ntree = 100,
  trControl = caret::trainControl(method = "none"),
  tuneGrid = data.frame(mtry = 2)
)

# 创建LIME解释器
lime_explainer <- lime::lime(
  x = lime_data,
  model = caret_model
)

# 解释单个样本
lime_explanation <- lime::explain(
  x = lime_data[1:3, ],
  explainer = lime_explainer,
  n_features = 5,
  n_labels = 1
)

# 可视化
plot_features(lime_explanation)

LIME vs SHAP对比

对比维度	LIME	SHAP
理论基础	局部线性近似	博弈论
计算速度	较快	较慢
一致性	可能不稳定	保证一致性
可解释性	局部解释	全局+局部解释

12.6 替代模型

基本原理

替代模型（Surrogate Model） 用简单可解释的模型近似复杂模型：

步骤： 1. 用复杂模型预测所有样本 2. 用简单模型（如决策树）拟合预测结果 3. 解释简单模型

# 替代模型示例
set.seed(42)

# 获取复杂模型的预测
complex_predictions <- predict(rf_model, disease_data[, -7], type = "prob")[, 2]

# 用决策树拟合预测结果
surrogate_tree <- rpart::rpart(
  complex_predictions ~ age + blood_pressure + cholesterol + bmi + glucose + smoking,
  data = disease_data,
  method = "anova",
  control = rpart::rpart.control(maxdepth = 3)
)

# 可视化决策树
rpart.plot::rpart.plot(surrogate_tree, 
                        main = "替代模型：决策树近似随机森林",
                        type = 2,
                        extra = 101)

12.7 模型公平性评估

公平性概念

模型公平性（Model Fairness） 确保模型对不同群体没有歧视：

公平性指标	定义
人口统计学均等	不同群体的预测率相同
机会均等	不同群体的真阳性率相同
预测均等	不同群体的精确率相同

# 公平性评估示例
# 假设我们关注模型对不同年龄组的公平性

# 创建年龄组
disease_data$age_group <- cut(disease_data$age, 
                               breaks = c(0, 40, 60, 100),
                               labels = c("Young", "Middle", "Senior"))

# 获取预测
disease_data$predicted <- predict(rf_model, disease_data[, -7])

# 计算各年龄组的性能
fairness_metrics <- disease_data %>%
  dplyr::group_by(age_group) %>%
  dplyr::summarise(
    accuracy = mean(predicted == disease),
    disease_rate = mean(disease == "Disease"),
    predicted_rate = mean(predicted == "Disease")
  )

print("各年龄组的模型表现：")

## [1] "各年龄组的模型表现："

print(fairness_metrics)

## # A tibble: 3 × 4
##   age_group accuracy disease_rate predicted_rate
##   <fct>        <dbl>        <dbl>          <dbl>
## 1 Young            1       0.103          0.103 
## 2 Middle           1       0.0452         0.0452
## 3 Senior           1       0.104          0.104

# 可视化
ggplot2::ggplot(fairness_metrics, ggplot2::aes(x = age_group, y = accuracy, fill = age_group)) +
  ggplot2::geom_bar(stat = "identity") +
  ggplot2::labs(
    title = "模型公平性评估：各年龄组准确率",
    x = "年龄组", y = "准确率"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "none")

意义与使用场景

核心意义：公平性评估确保模型对所有患者群体公平，避免医疗歧视。

举例：不同种族的疾病预测公平性

为什么需要公平性评估？

如果模型主要用某一种族的数据训练，可能对其他种族预测不准，导致医疗不公平。

研究目的是什么？

发现偏见：识别模型对某些群体的歧视
改进模型：调整模型提高公平性
合规审查：满足监管要求

可能得到什么样的结果？

可能发现模型对老年患者的准确率较低，需要收集更多老年患者数据重新训练。

需要注意什么？

公平性和准确性可能存在权衡，需要根据具体场景平衡。

第十三章时间序列机器学习

13.1 时间序列监督学习范式

时间序列特点

时间序列数据 具有时间依赖性，样本间存在顺序关系：

特点	说明
时间依赖	当前值与历史值相关
趋势	长期上升或下降趋势
季节性	周期性重复模式
周期性	非固定周期的波动
不规则性	随机波动

滑动窗口构建特征

核心思想：将时间序列转换为监督学习问题，用历史值预测未来值。

# 时间序列监督学习示例
set.seed(42)

# 创建模拟的时间序列数据：患者每日血糖监测
n_days <- 200
dates <- seq(as.Date("2023-01-01"), by = "day", length.out = n_days)

# 生成带有趋势和季节性的血糖数据
trend <- seq(100, 130, length.out = n_days)
seasonal <- 10 * sin(2 * pi * (1:n_days) / 30)  # 月周期
noise <- rnorm(n_days, 0, 5)
glucose <- trend + seasonal + noise

ts_data <- data.frame(
  date = dates,
  glucose = glucose
)

# 可视化原始时间序列
ggplot2::ggplot(ts_data, ggplot2::aes(x = date, y = glucose)) +
  ggplot2::geom_line(color = "steelblue") +
  ggplot2::labs(
    title = "患者每日血糖监测时间序列",
    x = "日期", y = "血糖值"
  ) +
  ggplot2::theme_minimal()

# 构建滑动窗口特征
create_lagged_features <- function(data, lag = 7) {
  n <- nrow(data)
  result <- data.frame(date = data$date[(lag+1):n])
  
  for (i in 1:lag) {
    result[[paste0("lag_", i)]] <- data$glucose[(lag+1-i):(n-i)]
  }
  
  result$target <- data$glucose[(lag+1):n]
  return(result)
}

# 创建滞后特征
lagged_data <- create_lagged_features(ts_data, lag = 7)

# 查看数据结构
head(lagged_data)

##         date    lag_1    lag_2    lag_3    lag_4    lag_5    lag_6    lag_7
## 1 2023-01-08 118.4074 109.7337 111.2846 111.0480 107.9950 101.3946 108.9339
## 2 2023-01-09 110.5272 118.4074 109.7337 111.2846 111.0480 107.9950 101.3946
## 3 2023-01-10 120.8087 110.5272 118.4074 109.7337 111.2846 111.0480 107.9950
## 4 2023-01-11 109.7035 120.8087 110.5272 118.4074 109.7337 111.2846 111.0480
## 5 2023-01-12 115.4633 109.7035 120.8087 110.5272 118.4074 109.7337 111.2846
## 6 2023-01-13 118.9694 115.4633 109.7035 120.8087 110.5272 118.4074 109.7337
##      target
## 1 110.52720
## 2 120.80871
## 3 109.70347
## 4 115.46333
## 5 118.96937
## 6  98.93211

意义与使用场景

核心意义：滑动窗口将时间序列转换为监督学习问题，可以使用传统机器学习方法。

举例：血糖预测

为什么这种数据适合用滑动窗口？

血糖值具有时间依赖性，今天的血糖与过去几天的血糖相关。滑动窗口可以捕捉这种依赖关系。

研究目的是什么？

短期预测：预测明天的血糖值
异常预警：预测血糖异常波动
干预评估：评估治疗效果

可能得到什么样的结果？

模型可能准确预测未来1-3天的血糖值，帮助患者提前调整饮食和用药。

需要注意什么？

滑动窗口大小需要根据数据特点选择。窗口太小可能丢失信息，太大可能引入噪声。

13.2 时间序列特征工程

滞后特征

# 滞后特征示例
# 添加多个滞后特征
ts_data$glucose_lag1 <- dplyr::lag(ts_data$glucose, 1)
ts_data$glucose_lag7 <- dplyr::lag(ts_data$glucose, 7)

# 可视化滞后关系
p1 <- ggplot2::ggplot(ts_data, ggplot2::aes(x = glucose_lag1, y = glucose)) +
  ggplot2::geom_point(alpha = 0.5) +
  ggplot2::geom_smooth(method = "lm", se = FALSE, color = "red") +
  ggplot2::labs(title = "1天滞后相关性", x = "昨天血糖", y = "今天血糖") +
  ggplot2::theme_minimal()

p2 <- ggplot2::ggplot(ts_data, ggplot2::aes(x = glucose_lag7, y = glucose)) +
  ggplot2::geom_point(alpha = 0.5) +
  ggplot2::geom_smooth(method = "lm", se = FALSE, color = "red") +
  ggplot2::labs(title = "7天滞后相关性", x = "7天前血糖", y = "今天血糖") +
  ggplot2::theme_minimal()

p1 + p2 + patchwork::plot_annotation(title = "滞后特征与目标的关系")

滚动统计特征

# 滚动统计特征
# 计算滚动均值和滚动标准差
ts_data$rolling_mean_7 <- zoo::rollmean(ts_data$glucose, k = 7, fill = NA, align = "right")
ts_data$rolling_sd_7 <- zoo::rollapply(ts_data$glucose, width = 7, FUN = sd, fill = NA, align = "right")

# 可视化滚动统计
ggplot2::ggplot(ts_data, ggplot2::aes(x = date)) +
  ggplot2::geom_line(ggplot2::aes(y = glucose), color = "gray50", alpha = 0.5) +
  ggplot2::geom_line(ggplot2::aes(y = rolling_mean_7), color = "steelblue", linewidth = 1) +
  ggplot2::geom_ribbon(ggplot2::aes(ymin = rolling_mean_7 - rolling_sd_7, 
                                     ymax = rolling_mean_7 + rolling_sd_7),
                       alpha = 0.2, fill = "steelblue") +
  ggplot2::labs(
    title = "滚动统计特征：7天均值和标准差",
    subtitle = "阴影区域为均值±标准差",
    x = "日期", y = "血糖值"
  ) +
  ggplot2::theme_minimal()

季节特征

# 提取季节特征
ts_data$day_of_week <- weekdays(ts_data$date)
ts_data$day_of_month <- as.numeric(format(ts_data$date, "%d"))
ts_data$month <- as.numeric(format(ts_data$date, "%m"))
ts_data$is_weekend <- ts_data$day_of_week %in% c("Saturday", "Sunday")

# 可视化周内模式
weekly_pattern <- ts_data %>%
  dplyr::group_by(day_of_week) %>%
  dplyr::summarise(
    mean_glucose = mean(glucose, na.rm = TRUE),
    sd_glucose = sd(glucose, na.rm = TRUE)
  )

# 设置星期顺序
weekly_pattern$day_of_week <- factor(weekly_pattern$day_of_week, 
                                      levels = c("Monday", "Tuesday", "Wednesday", 
                                                "Thursday", "Friday", "Saturday", "Sunday"))

ggplot2::ggplot(weekly_pattern, ggplot2::aes(x = day_of_week, y = mean_glucose)) +
  ggplot2::geom_bar(stat = "identity", fill = "steelblue") +
  ggplot2::geom_errorbar(ggplot2::aes(ymin = mean_glucose - sd_glucose, 
                                       ymax = mean_glucose + sd_glucose),
                         width = 0.2) +
  ggplot2::labs(
    title = "周内血糖模式",
    x = "星期", y = "平均血糖"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1))

13.3 时间序列交叉验证

时间序列分割

重要原则：不能用未来数据预测过去，必须保持时间顺序。

# 时间序列交叉验证示例
library(forecast)

# 创建时间序列对象
ts_glucose <- ts(ts_data$glucose, frequency = 30)  # 假设月周期

# 使用forecast包的时间序列CV
# 创建训练集和测试集
train_size <- floor(n_days * 0.8)
train_data <- ts_data$glucose[1:train_size]
test_data <- ts_data$glucose[(train_size + 1):n_days]

# 滚动预测
h <- 7  # 预测步长
n_windows <- 5  # 滚动窗口数

predictions <- data.frame()

for (i in 1:n_windows) {
  # 训练集结束点
  train_end <- train_size + (i - 1) * h
  
  if (train_end + h > n_days) break
  
  # 训练模型
  train_subset <- ts_data$glucose[1:train_end]
  fit <- forecast::auto.arima(train_subset)
  
  # 预测
  fc <- forecast::forecast(fit, h = h)
  
  # 存储结果
  predictions <- rbind(predictions, data.frame(
    window = i,
    date = ts_data$date[(train_end + 1):(train_end + h)],
    actual = ts_data$glucose[(train_end + 1):(train_end + h)],
    predicted = as.numeric(fc$mean)
  ))
}

# 计算预测误差
predictions$error <- predictions$actual - predictions$predicted
rmse <- sqrt(mean(predictions$error^2))
print(paste("RMSE:", round(rmse, 2)))

## [1] "RMSE: 8.4"

# 可视化预测结果
ggplot2::ggplot(predictions, ggplot2::aes(x = date)) +
  ggplot2::geom_line(ggplot2::aes(y = actual, color = "实际值")) +
  ggplot2::geom_line(ggplot2::aes(y = predicted, color = "预测值"), linetype = "dashed") +
  ggplot2::scale_color_manual(values = c("实际值" = "steelblue", "预测值" = "red")) +
  ggplot2::labs(
    title = "时间序列交叉验证预测结果",
    x = "日期", y = "血糖值", color = ""
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

13.4 基于树的时序模型

随机森林应用于时序

# 使用随机森林进行时间序列预测
set.seed(42)

# 准备特征数据
ts_features <- lagged_data[, -c(1, ncol(lagged_data))]  # 移除日期和目标
ts_target <- lagged_data$target

# 划分训练集和测试集
train_idx <- 1:floor(nrow(lagged_data) * 0.8)
train_x <- ts_features[train_idx, ]
train_y <- ts_target[train_idx]
test_x <- ts_features[-train_idx, ]
test_y <- ts_target[-train_idx]

# 训练随机森林
rf_ts <- randomForest::randomForest(
  x = as.matrix(train_x),
  y = train_y,
  ntree = 100,
  importance = TRUE
)

# 预测
rf_predictions <- predict(rf_ts, as.matrix(test_x))

# 评估
rf_rmse <- sqrt(mean((test_y - rf_predictions)^2))
print(paste("随机森林RMSE:", round(rf_rmse, 2)))

## [1] "随机森林RMSE: 6.43"

# 可视化
results_df <- data.frame(
  date = lagged_data$date[-train_idx],
  actual = test_y,
  predicted = rf_predictions
)

ggplot2::ggplot(results_df, ggplot2::aes(x = date)) +
  ggplot2::geom_line(ggplot2::aes(y = actual, color = "实际值")) +
  ggplot2::geom_line(ggplot2::aes(y = predicted, color = "预测值"), linetype = "dashed") +
  ggplot2::scale_color_manual(values = c("实际值" = "steelblue", "预测值" = "red")) +
  ggplot2::labs(
    title = "随机森林时间序列预测",
    x = "日期", y = "血糖值", color = ""
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

13.5 Prophet模型

Prophet原理

Prophet 是Facebook开发的时间序列预测模型，适合具有趋势和季节性的数据：

模型组成： \[y(t) = g(t) + s(t) + h(t) + \epsilon_t\]

\(g(t)\)：趋势项
\(s(t)\)：季节项
\(h(t)\)：节假日效应
\(\epsilon_t\)：误差项

# Prophet模型示例
library(prophet)

# 准备数据（Prophet需要ds和y两列）
prophet_data <- data.frame(
  ds = ts_data$date,
  y = ts_data$glucose
)

# 训练Prophet模型
prophet_model <- prophet::prophet(
  prophet_data,
  yearly.seasonality = FALSE,
  weekly.seasonality = TRUE,
  daily.seasonality = FALSE
)

# 创建未来数据框
future <- prophet::make_future_dataframe(prophet_model, periods = 30)

# 预测
prophet_forecast <- predict(prophet_model, future)

# 可视化
plot(prophet_model, prophet_forecast) +
  ggplot2::labs(
    title = "Prophet预测结果",
    subtitle = "蓝色区域为预测区间",
    x = "日期", y = "血糖值"
  ) +
  ggplot2::theme_minimal()

Prophet组件分解

# 查看Prophet的组件分解
prophet::prophet_plot_components(prophet_model, prophet_forecast)

意义与使用场景

核心意义：Prophet自动处理趋势和季节性，适合医学监测数据的预测。

举例：患者生命体征监测预测

为什么这种数据适合用Prophet？

生命体征数据通常有日周期（如血压早晚不同）、周周期（工作日vs周末）等季节性模式，Prophet可以自动捕捉这些模式。

研究目的是什么？

异常预警：预测正常范围，检测异常值
资源规划：预测医院就诊量
效果评估：评估干预措施的长期效果

可能得到什么样的结果？

Prophet可能预测未来一周的血糖趋势，帮助医生调整治疗方案。

需要注意什么？

Prophet假设趋势和季节性是可加的，对于乘性关系可能需要变换。

13.6 时序聚类与异常检测

时间序列聚类

# 时间序列聚类示例
set.seed(42)

# 创建多个患者的时间序列
n_patients <- 20
n_days <- 100

# 生成不同模式的时间序列
ts_list <- list()
for (i in 1:n_patients) {
  pattern_type <- sample(1:3, 1)
  
  if (pattern_type == 1) {
    # 上升趋势
    ts_list[[i]] <- seq(80, 120, length.out = n_days) + rnorm(n_days, 0, 5)
  } else if (pattern_type == 2) {
    # 下降趋势
    ts_list[[i]] <- seq(120, 80, length.out = n_days) + rnorm(n_days, 0, 5)
  } else {
    # 稳定趋势
    ts_list[[i]] <- rep(100, n_days) + rnorm(n_days, 0, 5)
  }
}

# 转换为矩阵
ts_matrix <- do.call(rbind, ts_list)

# 计算距离矩阵（动态时间规整）
# 简化：使用欧氏距离
dist_matrix <- dist(ts_matrix)

# 层次聚类
hc_result <- hclust(dist_matrix, method = "ward.D2")

# 可视化
plot(hc_result, main = "时间序列聚类树状图", xlab = "患者", ylab = "距离")

时间序列异常检测

# 时间序列异常检测
# 使用预测残差检测异常

# 计算预测残差
residuals <- test_y - rf_predictions

# 定义异常阈值（3倍标准差）
threshold <- 3 * sd(residuals)

# 识别异常点
anomalies <- which(abs(residuals) > threshold)

print(paste("检测到", length(anomalies), "个异常点"))

## [1] "检测到 0 个异常点"

# 可视化异常
results_df$is_anomaly <- abs(residuals) > threshold

ggplot2::ggplot(results_df, ggplot2::aes(x = date)) +
  ggplot2::geom_line(ggplot2::aes(y = actual), color = "steelblue") +
  ggplot2::geom_point(data = dplyr::filter(results_df, is_anomaly), 
                      ggplot2::aes(y = actual), color = "red", size = 3) +
  ggplot2::labs(
    title = "时间序列异常检测",
    subtitle = "红点为检测到的异常值",
    x = "日期", y = "血糖值"
  ) +
  ggplot2::theme_minimal()

意义与使用场景

核心意义：时间序列聚类发现相似模式，异常检测识别异常事件。

举例：患者病情轨迹分析

为什么这种数据适合用时序聚类？

不同患者的病情发展轨迹可能属于不同类型，如”快速恢复型”、“缓慢恢复型”、“恶化型”。时序聚类可以发现这些模式。

研究目的是什么？

模式发现：识别不同的病情发展模式
预后预测：根据轨迹类型预测预后
个体化治疗：针对不同轨迹类型制定方案

可能得到什么样的结果？

聚类可能发现3种血糖控制模式，每种模式需要不同的干预策略。

需要注意什么？

时间序列聚类对距离度量敏感，需要选择合适的距离方法（如DTW）。

第十四章文本挖掘基础

14.1 文本预处理

文本挖掘概述

文本挖掘（Text Mining） 从非结构化文本中提取有价值的信息：

步骤	说明
文本收集	获取原始文本数据
文本预处理	清洗和标准化文本
特征提取	将文本转换为数值特征
模型构建	应用机器学习方法

分词

分词（Tokenization） 将文本切分为词语或标记：

# 中文分词示例
# 使用简单的字符串分割方法演示分词概念
# 实际应用中可使用jiebaR等专门包

# 创建模拟的医疗文本数据
medical_texts <- c(
  "患者主诉头痛三天，伴有恶心呕吐症状",
  "既往有高血压病史五年，规律服药控制",
  "体格检查显示血压偏高，心率正常",
  "建议进行头部CT检查排除颅内病变"
)

# 简单分词：按标点和空格分割
simple_tokenize <- function(text) {
  # 使用正则表达式分割
  tokens <- unlist(strsplit(text, "[，。、；：！？\\s]+"))
  tokens <- tokens[nchar(tokens) > 0]
  paste(tokens, collapse = " | ")
}

# 进行分词
tokenized_texts <- sapply(medical_texts, simple_tokenize)

# 查看分词结果
print("分词结果：")

## [1] "分词结果："

for (i in seq_along(medical_texts)) {
  cat("原文:", medical_texts[i], "\n")
  cat("分词:", tokenized_texts[i], "\n\n")
}

## 原文: 患者主诉头痛三天，伴有恶心呕吐症状 
## 分词: 患者主诉头痛三天 | 伴有恶心呕吐症状 
## 
## 原文: 既往有高血压病史五年，规律服药控制 
## 分词: 既往有高血压病史五年 | 规律服药控制 
## 
## 原文: 体格检查显示血压偏高，心率正常 
## 分词: 体格检查显示血压偏高 | 心率正常 
## 
## 原文: 建议进行头部CT检查排除颅内病变 
## 分词: 建议进行头部CT检查排除颅内病变

停用词处理

停用词（Stop Words） 是高频但信息量小的词，如”的”、“是”等。

# 停用词处理
# 定义停用词列表
stop_words <- c("的", "是", "有", "在", "和", "了", "进行", "显示", "伴有")

# 移除停用词（简单示例）
remove_stopwords <- function(text, stop_words) {
  tokens <- unlist(strsplit(text, "[，。、；：！？\\s]+"))
  tokens <- tokens[!tokens %in% stop_words]
  tokens <- tokens[nchar(tokens) > 0]
  paste(tokens, collapse = " | ")
}

# 应用停用词移除
cleaned_texts <- sapply(medical_texts, function(x) {
  remove_stopwords(x, stop_words)
})

print("移除停用词后：")

## [1] "移除停用词后："

for (i in seq_along(cleaned_texts)) {
  cat(cleaned_texts[i], "\n")
}

## 患者主诉头痛三天 | 伴有恶心呕吐症状 
## 既往有高血压病史五年 | 规律服药控制 
## 体格检查显示血压偏高 | 心率正常 
## 建议进行头部CT检查排除颅内病变

词干提取与词形还原

词干提取（Stemming）：将词还原为词干 词形还原（Lemmatization）：将词还原为词典形式

# 英文文本处理示例
library(SnowballC)

# 英文医疗文本
english_texts <- c(
  "The patient was diagnosed with diabetes",
  "Patients are receiving treatments for various diseases",
  "The medication effectively reduced the symptoms"
)

# 创建语料库
corpus <- tm::Corpus(tm::VectorSource(english_texts))

# 文本预处理：转小写、移除标点、移除停用词
corpus <- tm::tm_map(corpus, tm::content_transformer(tolower))
corpus <- tm::tm_map(corpus, tm::removePunctuation)
corpus <- tm::tm_map(corpus, tm::removeWords, tm::stopwords("en"))

# 词干提取
corpus <- tm::tm_map(corpus, tm::stemDocument)

# 查看处理结果
print("预处理后的文本：")

## [1] "预处理后的文本："

for (i in 1:length(corpus)) {
  cat(corpus[[i]]$content, "\n")
}

## patient diagnos diabet 
## patient receiv treatment various diseas 
## medic effect reduc symptom

意义与使用场景

核心意义：文本预处理是文本挖掘的基础，直接影响后续分析质量。

举例：电子病历文本处理

为什么这种数据适合用文本预处理？

电子病历包含大量非结构化文本，如主诉、病程记录等。预处理可以将这些文本转换为可分析的结构化数据。

研究目的是什么？

信息提取：从病历中提取关键信息
文本分类：自动分类病历类型
知识发现：发现疾病与症状的关联

可能得到什么样的结果？

预处理后可以从”患者主诉头痛三天”中提取”头痛”作为症状，“三天”作为持续时间。

需要注意什么？

中文分词质量对后续分析影响很大，需要选择合适的分词工具和自定义词典。

14.2 词袋模型与TF-IDF

词袋模型（BoW）

词袋模型（Bag of Words） 将文本表示为词频向量，忽略词序：

# 词袋模型示例
# 使用tm包创建文档-词矩阵

# 创建模拟的医疗文本数据集
documents <- c(
  "高血压 糖尿病 冠心病",
  "高血压 冠心病 治疗",
  "糖尿病 肾病 治疗",
  "冠心病 糖尿病 高血压 治疗"
)

# 创建语料库
corpus <- tm::Corpus(tm::VectorSource(documents))

# 创建文档-词矩阵
dtm <- tm::DocumentTermMatrix(corpus)

# 查看矩阵
print("文档-词矩阵：")

## [1] "文档-词矩阵："

tm::inspect(dtm)

## <<DocumentTermMatrix (documents: 4, terms: 5)>>
## Non-/sparse entries: 13/7
## Sparsity           : 35%
## Maximal term length: 3
## Weighting          : term frequency (tf)
## Sample             :
##     Terms
## Docs 高血压 冠心病 肾病 糖尿病 治疗
##    1      1      1    0      1    0
##    2      1      1    0      0    1
##    3      0      0    1      1    1
##    4      1      1    0      1    1

# 转换为普通矩阵查看
dtm_matrix <- as.matrix(dtm)
print(dtm_matrix)

##     Terms
## Docs 冠心病 糖尿病 高血压 治疗 肾病
##    1      1      1      1    0    0
##    2      1      0      1    1    0
##    3      0      1      0    1    1
##    4      1      1      1    1    0

TF-IDF

TF-IDF（Term Frequency-Inverse Document Frequency） 衡量词的重要性：

\[\text{TF-IDF}(t, d) = \text{TF}(t, d) \times \text{IDF}(t)\]

\[\text{IDF}(t) = \log\frac{N}{df(t)}\]

其中： - \(\text{TF}(t, d)\)：词t在文档d中的频率 - \(N\)：文档总数 - \(df(t)\)：包含词t的文档数

# TF-IDF计算
# 创建带TF-IDF权重的文档-词矩阵
dtm_tfidf <- tm::weightTfIdf(dtm)

# 查看TF-IDF矩阵
tfidf_matrix <- as.matrix(dtm_tfidf)
print("TF-IDF矩阵：")

## [1] "TF-IDF矩阵："

print(round(tfidf_matrix, 3))

##     Terms
## Docs 冠心病 糖尿病 高血压  治疗  肾病
##    1  0.138  0.138  0.138 0.000 0.000
##    2  0.138  0.000  0.138 0.138 0.000
##    3  0.000  0.138  0.000 0.138 0.667
##    4  0.104  0.104  0.104 0.104 0.000

# 找出每个文档中最重要的词
for (i in 1:nrow(tfidf_matrix)) {
  top_words <- names(sort(tfidf_matrix[i, ], decreasing = TRUE))[1:2]
  print(paste("文档", i, "最重要的词:", paste(top_words, collapse = ", ")))
}

## [1] "文档 1 最重要的词: 冠心病, 糖尿病"
## [1] "文档 2 最重要的词: 冠心病, 高血压"
## [1] "文档 3 最重要的词: 肾病, 糖尿病"
## [1] "文档 4 最重要的词: 冠心病, 糖尿病"

意义与使用场景

核心意义：TF-IDF能够识别文档中的关键词，是文本分析的基础方法。

举例：医学文献关键词提取

为什么这种数据适合用TF-IDF？

医学文献中，专业术语（如疾病名称、药物名称）在某篇文献中频繁出现，但在所有文献中出现较少，TF-IDF可以识别这些关键词。

研究目的是什么？

关键词提取：自动提取文献关键词
文献检索：提高检索准确性
文献分类：根据关键词分类文献

可能得到什么样的结果？

TF-IDF可能识别出”糖尿病”、“胰岛素抵抗”等作为某篇文献的关键词。

需要注意什么？

TF-IDF忽略词序和语义，对于需要理解上下文的任务可能不够。

14.3 N-gram模型

N-gram原理

N-gram 是连续的N个词的序列，可以捕捉局部词序信息：

N-gram类型	示例
Unigram (1-gram)	“高血压”, “糖尿病”
Bigram (2-gram)	“高血压患者”, “糖尿病治疗”
Trigram (3-gram)	“高血压患者治疗”

# N-gram示例
library(tidytext)

# 创建医疗文本数据框
text_df <- data.frame(
  doc_id = 1:4,
  text = c(
    "高血压患者需要长期治疗",
    "糖尿病患者应控制饮食",
    "高血压和糖尿病都是慢性病",
    "慢性病需要长期管理和治疗"
  )
)

# 分词
tokens <- tidytext::unnest_tokens(text_df, word, text)

# 查看unigram结果
print("Unigram结果：")

## [1] "Unigram结果："

head(tokens, 10)

##    doc_id   word
## 1       1 高血压
## 2       1   患者
## 3       1   需要
## 4       1   长期
## 5       1   治疗
## 6       2 糖尿病
## 7       2   患者
## 8       2     应
## 9       2   控制
## 10      2   饮食

# 创建bigram
bigrams <- text_df %>%
  tidytext::unnest_tokens(bigram, text, token = "ngrams", n = 2)

# 查看bigram结果
print("Bigram结果：")

## [1] "Bigram结果："

head(bigrams, 10)

##    doc_id      bigram
## 1       1 高血压 患者
## 2       1   患者 需要
## 3       1   需要 长期
## 4       1   长期 治疗
## 5       2 糖尿病 患者
## 6       2     患者 应
## 7       2     应 控制
## 8       2   控制 饮食
## 9       3   高血压 和
## 10      3   和 糖尿病

# 统计bigram频率
bigram_counts <- bigrams %>%
  dplyr::count(bigram, sort = TRUE)

print("最常见的bigram：")

## [1] "最常见的bigram："

head(bigram_counts, 10)

##         bigram n
## 1    需要 长期 2
## 2      和 治疗 1
## 3    和 糖尿病 1
## 4      应 控制 1
## 5      患者 应 1
## 6    患者 需要 1
## 7  慢性病 需要 1
## 8    控制 饮食 1
## 9      管理 和 1
## 10 糖尿病 患者 1

14.4 主题模型：LDA

LDA原理

潜在狄利克雷分配（LDA） 是一种生成式概率模型，用于发现文档集合中的潜在主题：

核心假设： 1. 每个文档是多个主题的混合 2. 每个主题是词语的概率分布

# LDA主题模型示例
library(topicmodels)

# 创建模拟的医疗文档集合
medical_docs <- c(
  "高血压 血压 心血管 治疗 药物 患者",
  "糖尿病 血糖 胰岛素 治疗 患者 饮食",
  "高血压 心血管 疾病 风险 预防 治疗",
  "糖尿病 并发症 肾病 视网膜 病变",
  "心血管 冠心病 心脏 血管 疾病",
  "血糖 胰岛素 糖尿病 控制 治疗"
)

# 创建语料库和文档-词矩阵
corpus <- tm::Corpus(tm::VectorSource(medical_docs))
dtm <- tm::DocumentTermMatrix(corpus)

# 训练LDA模型
lda_model <- topicmodels::LDA(dtm, k = 2, control = list(seed = 42))

# 查看主题-词分布
topics <- tidytext::tidy(lda_model, matrix = "beta")

# 找出每个主题最重要的词
top_terms <- topics %>%
  dplyr::group_by(topic) %>%
  dplyr::slice_max(beta, n = 5) %>%
  dplyr::ungroup() %>%
  dplyr::arrange(topic, -beta)

# 可视化
ggplot2::ggplot(top_terms, ggplot2::aes(x = reorder(term, beta), y = beta, fill = factor(topic))) +
  ggplot2::geom_bar(stat = "identity", show.legend = FALSE) +
  ggplot2::facet_wrap(~ topic, scales = "free") +
  ggplot2::coord_flip() +
  ggplot2::labs(
    title = "LDA主题模型：每个主题最重要的词",
    x = "词", y = "概率"
  ) +
  ggplot2::theme_minimal()

意义与使用场景

核心意义：LDA可以发现文档集合中的潜在主题，是文本分析的重要工具。

举例：医学文献主题发现

为什么这种数据适合用LDA？

医学文献数量庞大，人工分类困难。LDA可以自动发现文献中的主要研究主题，帮助研究者快速了解领域动态。

研究目的是什么？

主题发现：识别文献中的主要研究主题
文献分类：根据主题自动分类文献
趋势分析：分析研究主题的时间变化

可能得到什么样的结果？

LDA可能发现”心血管疾病”、“糖尿病治疗”等主题，每个主题包含相关的关键词。

需要注意什么？

主题数量K需要预先指定，可以使用困惑度等方法选择最优K值。

14.5 文本分类

文本分类概述

文本分类 将文档分配到预定义的类别：

# 文本分类示例
# 使用朴素贝叶斯分类器

set.seed(42)

# 创建模拟的医疗文本分类数据
text_class_data <- data.frame(
  text = c(
    "患者血压偏高，诊断为高血压",
    "血糖控制不佳，糖尿病诊断明确",
    "心脏彩超显示心功能下降",
    "胰岛素治疗糖尿病效果良好",
    "高血压药物治疗后血压稳定",
    "冠心病患者需要长期服药"
  ),
  category = factor(c("高血压", "糖尿病", "心血管", "糖尿病", "高血压", "心血管"))
)

# 创建语料库
corpus <- tm::Corpus(tm::VectorSource(text_class_data$text))

# 预处理
corpus <- tm::tm_map(corpus, tm::content_transformer(tolower))
corpus <- tm::tm_map(corpus, tm::removePunctuation)

# 创建文档-词矩阵
dtm <- tm::DocumentTermMatrix(corpus)

# 划分训练集和测试集
train_idx <- sample(1:nrow(text_class_data), size = 4)
test_idx <- setdiff(1:nrow(text_class_data), train_idx)

train_dtm <- dtm[train_idx, ]
test_dtm <- dtm[test_idx, ]

train_labels <- text_class_data$category[train_idx]
test_labels <- text_class_data$category[test_idx]

# 使用朴素贝叶斯分类器
# 训练模型
classifier <- e1071::naiveBayes(as.matrix(train_dtm), train_labels)

# 预测
predictions <- predict(classifier, as.matrix(test_dtm))

# 评估
print("预测结果：")

## [1] "预测结果："

print(data.frame(
  实际 = test_labels,
  预测 = predictions
))

##     实际   预测
## 1 糖尿病 高血压
## 2 心血管 高血压

# 计算准确率
accuracy <- mean(predictions == test_labels)
print(paste("准确率:", round(accuracy, 2)))

## [1] "准确率: 0"

14.6 词嵌入基础

词嵌入概念

词嵌入（Word Embedding） 将词映射到低维向量空间，捕捉语义关系：

方法	特点
Word2Vec	预测上下文或目标词
GloVe	基于全局共现矩阵
FastText	考虑子词信息

核心思想：语义相似的词在向量空间中距离较近。

# 词嵌入概念演示
# 注意：实际应用需要大量文本数据训练

# 模拟词嵌入向量
set.seed(42)

# 创建模拟的医疗词向量
words <- c("高血压", "血压", "心血管", "糖尿病", "血糖", "胰岛素", 
           "治疗", "药物", "患者", "疾病")

# 模拟2维词向量（实际通常是100-300维）
word_vectors <- data.frame(
  word = words,
  dim1 = c(1.2, 1.1, 1.0, -1.0, -1.1, -1.2, 0.1, 0.2, 0.0, 0.1),
  dim2 = c(0.5, 0.4, 0.6, 0.5, 0.4, 0.6, -0.8, -0.7, -0.9, -0.8)
)

# 可视化词向量
ggplot2::ggplot(word_vectors, ggplot2::aes(x = dim1, y = dim2)) +
  ggplot2::geom_point(size = 3, color = "steelblue") +
  ggplot2::geom_text(ggplot2::aes(label = word), vjust = -1, size = 3) +
  ggplot2::labs(
    title = "词嵌入可视化（模拟）",
    subtitle = "语义相似的词距离较近",
    x = "维度1", y = "维度2"
  ) +
  ggplot2::theme_minimal()

意义与使用场景

核心意义：词嵌入捕捉词语的语义关系，是现代NLP的基础。

举例：医学文本语义分析

为什么这种数据适合用词嵌入？

医学文本中，“高血压”和”血压”语义相关，词嵌入可以捕捉这种关系，用于文本相似度计算、信息检索等。

研究目的是什么？

语义搜索：根据语义相似度检索文档
文本聚类：基于语义相似度聚类文档
特征提取：用词向量作为下游任务的特征

可能得到什么样的结果？

词嵌入可能发现”糖尿病”与”胰岛素”距离较近，反映它们的语义关联。

需要注意什么？

词嵌入需要大量文本训练，对于医学领域建议使用领域语料训练或微调。

第十五章推荐系统基础

15.1 推荐系统概述

什么是推荐系统？

推荐系统（Recommender System） 预测用户对物品的偏好，提供个性化推荐：

类型	说明	示例
协同过滤	基于用户行为相似性	Netflix电影推荐
基于内容	基于物品特征相似性	新闻推荐
混合方法	结合多种方法	亚马逊商品推荐

应用场景	说明
治疗方案推荐	根据相似患者推荐治疗方案
药物推荐	根据病情推荐合适药物
健康内容推荐	推荐相关健康知识
医生推荐	推荐合适的医生

意义与使用场景

核心意义：推荐系统帮助用户从海量信息中发现感兴趣的内容。

举例：治疗方案推荐

为什么这种数据适合用推荐系统？

不同患者对同一治疗方案的反应可能不同。推荐系统可以根据相似患者的治疗历史，为新患者推荐可能有效的治疗方案。

研究目的是什么？

个性化治疗：为患者推荐最适合的治疗方案
辅助决策：为医生提供决策支持
提高疗效：选择最可能有效的治疗

可能得到什么样的结果？

系统可能推荐”对于类似患者，方案A的有效率为80%“。

需要注意什么？

医疗推荐需要谨慎，推荐结果应作为参考而非决定。

15.2 协同过滤

基于用户的协同过滤

原理：找到与目标用户相似的用户，推荐这些用户喜欢的物品。

# 协同过滤示例
library(recommenderlab)

# 创建模拟的用户-物品评分矩阵
# 假设是患者-治疗方案评分
set.seed(42)

# 创建评分矩阵
n_users <- 20
n_items <- 10

# 模拟评分数据（1-5分）
ratings <- matrix(sample(c(1:5, NA), n_users * n_items, replace = TRUE, 
                          prob = c(0.1, 0.1, 0.2, 0.3, 0.2, 0.1)),
                  nrow = n_users, ncol = n_items)
colnames(ratings) <- paste0("治疗", 1:n_items)
rownames(ratings) <- paste0("患者", 1:n_users)

# 转换为realRatingMatrix
rating_matrix <- as(ratings, "realRatingMatrix")

# 查看数据摘要
print("评分矩阵摘要：")

## [1] "评分矩阵摘要："

print(dim(rating_matrix))

## [1] 20 10

print("评分分布：")

## [1] "评分分布："

print(table(getRatings(rating_matrix)))

## 
##  1  2  3  4  5 
## 26 14 40 56 33

# 训练基于用户的协同过滤模型
ubcf_model <- recommenderlab::Recommender(
  rating_matrix, 
  method = "UBCF",
  param = list(method = "cosine", nn = 5)
)

# 为前3个患者推荐治疗
recommendations <- recommenderlab::predict(ubcf_model, rating_matrix[1:3], n = 3)

# 查看推荐结果
print("推荐结果：")

## [1] "推荐结果："

rec_list <- as(recommendations, "list")
for (i in 1:3) {
  cat("患者", i, "的推荐治疗:", rec_list[[i]], "\n")
}

## 患者 1 的推荐治疗:  
## 患者 2 的推荐治疗:  
## 患者 3 的推荐治疗: 治疗4

基于物品的协同过滤

原理：找到与用户喜欢的物品相似的物品进行推荐。

# 基于物品的协同过滤
ibcf_model <- recommenderlab::Recommender(
  rating_matrix, 
  method = "IBCF",
  param = list(k = 5)
)

# 推荐
recommendations_ibcf <- recommenderlab::predict(ibcf_model, rating_matrix[1:3], n = 3)

# 查看推荐结果
print("基于物品的协同过滤推荐结果：")

## [1] "基于物品的协同过滤推荐结果："

rec_list_ibcf <- as(recommendations_ibcf, "list")
for (i in 1:3) {
  cat("患者", i, "的推荐治疗:", rec_list_ibcf[[i]], "\n")
}

## 患者 1 的推荐治疗:  
## 患者 2 的推荐治疗:  
## 患者 3 的推荐治疗: 治疗4

15.3 矩阵分解

SVD矩阵分解

原理：将用户-物品矩阵分解为用户矩阵和物品矩阵的乘积：

\[R \approx U \times V^T\]

# 矩阵分解示例
# 使用SVD进行推荐

# 将评分矩阵转换为数值矩阵（NA用0填充）
ratings_filled <- ratings
ratings_filled[is.na(ratings_filled)] <- 0

# 执行SVD
svd_result <- svd(ratings_filled)

# 选择前k个奇异值
k <- 5
U_k <- svd_result$u[, 1:k]
V_k <- svd_result$v[, 1:k]
D_k <- diag(svd_result$d[1:k])

# 重构评分矩阵
ratings_predicted <- U_k %*% D_k %*% t(V_k)

# 查看预测评分
print("预测评分矩阵（前5行5列）：")

## [1] "预测评分矩阵（前5行5列）："

print(round(ratings_predicted[1:5, 1:5], 2))

##      [,1] [,2] [,3] [,4] [,5]
## [1,] 1.32 1.47 4.95 2.34 4.33
## [2,] 0.88 2.63 4.60 2.59 5.23
## [3,] 2.69 1.46 3.86 2.28 4.22
## [4,] 2.29 1.98 1.39 1.73 2.56
## [5,] 3.63 2.77 4.06 3.21 0.97

15.4 基于内容的推荐

基本原理

基于内容的推荐 根据物品特征和用户偏好进行推荐：

# 基于内容的推荐示例
# 假设治疗方案有特征描述

# 创建治疗方案特征
treatment_features <- data.frame(
  治疗 = paste0("治疗", 1:6),
  类型 = c("药物", "手术", "药物", "物理", "药物", "手术"),
  侵入性 = c(1, 5, 1, 2, 1, 4),
  成本 = c(2, 5, 3, 2, 4, 5),
  恢复时间 = c(1, 5, 2, 2, 3, 4)
)

print("治疗方案特征：")

## [1] "治疗方案特征："

print(treatment_features)

##    治疗 类型 侵入性 成本 恢复时间
## 1 治疗1 药物      1    2        1
## 2 治疗2 手术      5    5        5
## 3 治疗3 药物      1    3        2
## 4 治疗4 物理      2    2        2
## 5 治疗5 药物      1    4        3
## 6 治疗6 手术      4    5        4

# 用户偏好（假设用户偏好低侵入性、低成本）
user_preference <- c(侵入性 = -0.5, 成本 = -0.3, 恢复时间 = -0.2)

# 计算每个治疗的得分
scores <- as.matrix(treatment_features[, 3:5]) %*% user_preference
treatment_features$推荐得分 <- scores

# 排序推荐
treatment_features <- treatment_features[order(-treatment_features$推荐得分), ]

print("推荐排序：")

## [1] "推荐排序："

print(treatment_features[, c("治疗", "类型", "推荐得分")])

##    治疗 类型 推荐得分
## 1 治疗1 药物     -1.3
## 3 治疗3 药物     -1.8
## 4 治疗4 物理     -2.0
## 5 治疗5 药物     -2.3
## 6 治疗6 手术     -4.3
## 2 治疗2 手术     -5.0

15.5 冷启动问题

冷启动类型

类型	说明	解决方案
新用户	没有用户历史数据	基于人口统计学推荐
新物品	没有物品评分数据	基于内容推荐
系统冷启动	系统刚上线	使用热门推荐

# 冷启动解决方案示例

# 新用户：基于人口统计学特征推荐
user_profiles <- data.frame(
  用户 = paste0("用户", 1:5),
  年龄段 = c("青年", "中年", "老年", "青年", "中年"),
  性别 = c("男", "女", "男", "女", "男"),
  偏好治疗 = c("治疗1", "治疗2", "治疗3", "治疗1", "治疗2")
)

# 对于新用户，根据人口统计学特征推荐
new_user <- data.frame(年龄段 = "青年", 性别 = "男")

# 找到相似用户
similar_users <- user_profiles %>%
  dplyr::filter(年龄段 == new_user$年龄段, 性别 == new_user$性别)

print("相似用户的偏好：")

## [1] "相似用户的偏好："

print(similar_users)

##    用户 年龄段 性别 偏好治疗
## 1 用户1   青年   男    治疗1

15.6 评估指标

指标	公式	说明
RMSE	\(\sqrt{\frac{1}{n}\sum(r - \hat{r})^2}\)	预测误差
精确率@K	相关推荐数/K	前K个推荐的准确率
召回率@K	相关推荐数/总相关数	找到相关物品的比例
NDCG	考虑排序位置的指标	排序质量

意义与使用场景

核心意义：评估指标帮助选择最优的推荐算法和参数。

举例：治疗方案推荐系统评估

为什么需要评估指标？

不同的推荐算法可能适合不同的场景，需要通过评估指标选择最适合的算法。

研究目的是什么？

算法选择：比较不同算法的性能
参数调优：找到最优参数
系统优化：持续改进推荐质量

可能得到什么样的结果？

评估可能发现UBCF在精确率上优于IBCF，但IBCF计算效率更高。

需要注意什么？

离线评估指标可能与实际用户体验不一致，需要结合在线A/B测试。

总结

本教程涵盖了机器学习的核心内容：

章节	主题	核心方法
第九章	无监督学习	K-Means、DBSCAN、PCA、t-SNE
第十章	关联规则	Apriori、FP-Growth
第十一章	特征选择	过滤法、包装法、嵌入法
第十二章	模型解释	SHAP、LIME、PDP
第十三章	时间序列	滑动窗口、Prophet
第十四章	文本挖掘	TF-IDF、LDA、词嵌入
第十五章	推荐系统	协同过滤、矩阵分解

R语言与机器学习（下）

小天使格礼

目录

第九章 无监督学习

9.1 聚类方法概述

什么是聚类？

聚类方法分类

意义与使用场景

9.2 K-Means聚类

基本原理

可视化聚类结果

肘部法则选择K值

轮廓系数评估聚类质量

初始化方法

意义与使用场景

9.3 层次聚类

基本原理

凝聚式聚类

树状图解读

切割树状图获取聚类

意义与使用场景

9.4 DBSCAN密度聚类

基本原理

可视化DBSCAN结果

K-Means vs DBSCAN对比

参数选择

意义与使用场景

9.5 高斯混合模型（GMM）与EM算法

基本原理

可视化GMM结果

软聚类 vs 硬聚类

意义与使用场景

9.6 谱聚类

基本原理

意义与使用场景

9.7 聚类评估指标

内部评估指标

外部评估指标

意义与使用场景

9.8 主成分分析（PCA）降维

基本原理

可视化PCA结果

PCA双标图

意义与使用场景

9.9 流形学习：t-SNE与UMAP

t-SNE原理

UMAP原理

PCA vs t-SNE vs UMAP对比

意义与使用场景

9.10 自编码器降维概念

基本原理

与传统降维方法的比较

医学应用场景

9.11 异常检测方法

异常检测概述

孤立森林（Isolation Forest）

局部异常因子（LOF）

单类SVM

意义与使用场景

第十章 关联规则与模式挖掘

10.1 关联规则基本概念

什么是关联规则？

三个核心指标

意义与使用场景

10.2 Apriori算法

算法原理

查看和筛选规则

可视化关联规则

意义与使用场景

10.3 FP-Growth算法

算法原理

Apriori vs FP-Growth对比

10.4 频繁项集挖掘

频繁项集的类型

10.5 序列模式挖掘基础

什么是序列模式？

序列模式挖掘

意义与使用场景

第十一章 特征选择与降维

11.1 特征选择概述

第九章无监督学习

第十章关联规则与模式挖掘

第十一章特征选择与降维

第十二章模型解释与可解释性

第十三章时间序列机器学习

第十四章文本挖掘基础