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

install.packages(c(
  "dplyr", "tidyr", "data.table",
  "ggplot2", "patchwork", "RColorBrewer",
  "caret", "mlr3", "tidymodels",
  "MASS", "e1071", "class", "nnet",
  "rpart", "rpart.plot", "randomForest", "gbm", "xgboost",
  "pROC", "ROCR", "yardstick",
  "recipes", "dimRed", "RANN",
  "mlbench", "ISLR",
  "purrr", "broom", "knitr", "kableExtra"
))

目录

章节 主题 核心内容
第一章 机器学习基础概念 定义、学习类型、问题类型、核心概念
第二章 数据预处理与特征工程 特征缩放、编码、缺失值、特征选择
第三章 模型评估与验证 评估指标、混淆矩阵、ROC曲线、交叉验证
第四章 线性模型 线性回归、逻辑回归、正则化、SVM
第五章 树模型与集成学习 决策树、随机森林、GBDT、XGBoost
第六章 支持向量机与核方法 最大间隔、核函数、核技巧、SVR
第七章 最近邻与基于实例的方法 KNN、距离度量、加权KNN
第八章 贝叶斯与概率图模型 朴素贝叶斯、贝叶斯网络、高斯过程

第一章 机器学习基础概念

1.1 机器学习定义与核心思想

什么是机器学习?

机器学习(Machine Learning) 是人工智能的一个分支,其核心思想是让计算机从数据中”学习”规律,而不是由人类显式编程所有规则。

经典定义

美国计算机科学家Tom Mitchell给出的定义:“如果一个计算机程序在某类任务T上的性能P,随着经验E的积累而不断改善,则称该程序对任务T进行了学习。”

这个定义包含三个核心要素:

  1. 任务T(Task):程序要完成的工作,如分类、回归、聚类等
  2. 性能P(Performance):衡量任务完成质量的指标,如准确率、均方误差等
  3. 经验E(Experience):用于学习的数据,即训练样本

核心思想:从数据到决策

传统编程 vs 机器学习:

对比维度 传统编程 机器学习
输入 数据 + 规则 数据 + 答案
输出 答案 规则(模型)
核心逻辑 人类编写规则 算法自动学习规则
适用场景 规则明确、可穷举 规则复杂、难以显式表达

机器学习的本质

机器学习本质上是一个函数逼近问题:

\[\hat{f} \approx f\]

其中: - \(f\):真实的数据生成函数(未知) - \(\hat{f}\):学习到的模型(已知)

我们的目标是从假设空间中找到一个最优的 \(\hat{f}\),使得它在未见数据上的表现尽可能好。

# 机器学习本质示例:拟合一个未知函数
# 假设真实函数 f(x) = sin(x) + 噪声

# 生成模拟数据
set.seed(123)
n <- 100
x <- seq(0, 2 * pi, length.out = n)
y_true <- sin(x)  # 真实函数
y <- y_true + rnorm(n, mean = 0, sd = 0.2)  # 观测值(带噪声)

# 创建数据框
data <- data.frame(x = x, y = y, y_true = y_true)

# 用多项式回归"学习"这个函数
# 尝试不同阶数的多项式
model_3 <- lm(y ~ poly(x, degree = 3), data = data)
model_9 <- lm(y ~ poly(x, degree = 9), data = data)

# 预测
data$pred_3 <- predict(model_3)
data$pred_9 <- predict(model_9)

# 可视化:展示学习过程
ggplot2::ggplot(data, ggplot2::aes(x = x)) +
  ggplot2::geom_point(ggplot2::aes(y = y), alpha = 0.5, color = "gray50") +
  ggplot2::geom_line(ggplot2::aes(y = y_true), color = "black", linewidth = 1, linetype = "dashed") +
  ggplot2::geom_line(ggplot2::aes(y = pred_3, color = "3阶多项式"), linewidth = 1) +
  ggplot2::geom_line(ggplot2::aes(y = pred_9, color = "9阶多项式"), linewidth = 1) +
  ggplot2::scale_color_manual(values = c("3阶多项式" = "blue", "9阶多项式" = "red")) +
  ggplot2::labs(
    title = "机器学习本质:从数据中学习函数",
    subtitle = "黑色虚线为真实函数,灰点为观测数据",
    x = "x", y = "y", color = "模型"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

意义与使用场景

核心意义:机器学习让计算机能够从数据中自动发现规律,解决那些规则复杂、难以显式编程的问题。

举例:基于患者的临床指标预测糖尿病风险

为什么这种数据适合用机器学习?

糖尿病的发生与年龄、BMI、血糖、血压、家族史等多种因素相关,这些因素之间存在复杂的交互作用和非线性关系,难以用简单的规则描述。机器学习可以从大量历史病例数据中自动学习这些复杂关系,建立预测模型。

研究目的是什么?

  1. 风险评估:根据患者的各项指标,预测其患糖尿病的概率
  2. 早期干预:识别高风险人群,进行针对性的预防措施
  3. 辅助诊断:为临床医生提供决策支持

可能得到什么样的结果?

模型可能发现:年龄大于45岁、BMI超过28、空腹血糖高于6.1 mmol/L的患者,糖尿病风险显著升高。这些规律可能不是简单的线性组合,而是复杂的非线性关系。

需要注意什么?

机器学习模型是”黑箱”,其预测结果需要结合医学知识进行解释。同时,模型在训练数据上的表现好,不代表在新患者上也能准确预测,需要验证模型的泛化能力。


1.2 什么时候要用到机器学习

适用场景判断

机器学习并非万能,以下情况适合使用机器学习:

适用条件 说明 示例
规则复杂 问题涉及大量变量和复杂关系 疾病诊断、股票预测
数据充足 有足够的历史数据用于学习 医疗影像、交易记录
模式存在 数据中确实存在可学习的规律 垃圾邮件识别
容忍误差 允许一定程度的预测错误 推荐系统
重复任务 需要对大量相似任务做决策 信用评分

不适合使用机器学习的场景

不适用条件 说明 替代方案
规则简单明确 可以用简单的if-else解决 传统编程
数据不足 样本量太少无法学习 专家系统
需要完美准确 零容忍错误场景 规则引擎
需要完全解释 必须解释每个决策原因 统计模型
实时性要求极高 毫秒级响应要求 优化算法

意义与使用场景

核心意义:判断何时使用机器学习是项目成功的第一步,避免在不适合的场景浪费资源。

举例:医院是否需要用机器学习预测患者住院时长

为什么需要先判断适用性?

医院信息系统积累了大量患者数据,管理者想用机器学习预测住院时长以优化床位管理。但需要先评估:住院时长受病情、治疗方案、并发症等多种因素影响,规则复杂;历史数据充足(数万条记录);确实存在可学习的规律(同类疾病住院时长有分布特征);允许一定预测误差(用于规划而非精确安排)。这些条件都满足,适合使用机器学习。

研究目的是什么?

  1. 床位规划:根据预测结果提前安排床位
  2. 资源调配:预测高峰期,提前增加医护人员
  3. 患者沟通:给患者提供大致住院时长预期

可能得到什么样的结果?

模型可能发现:年龄大于65岁、有并发症、手术类型为某类的患者,平均住院时长显著延长。预测准确率可能达到70-80%,足以支持规划决策。

需要注意什么?

如果数据质量差(缺失值多、记录不准),或预测精度要求极高(如ICU床位必须精确到天),则不适合使用机器学习,应考虑其他方案。

1.3 监督学习 vs 无监督学习 vs 半监督学习 vs 强化学习

四种主要学习范式

机器学习根据学习方式和数据特点,主要分为四种范式:

学习类型 是否有标签 数据特点 典型任务
监督学习 带标签的训练数据 分类、回归
无监督学习 无标签数据 聚类、降维
半监督学习 部分 少量标签+大量无标签 分类
强化学习 延迟奖励 环境交互反馈 决策控制

监督学习(Supervised Learning)

定义:从带标签的训练数据中学习一个映射函数,将输入映射到输出。

数学表达

给定训练集 \(D = \{(x_1, y_1), (x_2, y_2), ..., (x_n, y_n)\}\),学习函数 \(f: X \rightarrow Y\),使得对新输入 \(x_{new}\),能预测其输出 \(y_{new}\)

核心特点: 1. 训练数据包含输入和对应的正确答案(标签) 2. 学习目标是建立输入到输出的映射关系 3. 可以用测试数据评估模型性能

# 监督学习示例:鸢尾花分类
# 加载鸢尾花数据集
data(iris)

# 查看数据结构
head(iris)
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa
# 可视化:不同种类的鸢尾花在特征空间中的分布
ggplot2::ggplot(iris, ggplot2::aes(x = Sepal.Length, y = Petal.Length, color = Species)) +
  ggplot2::geom_point(size = 3, alpha = 0.7) +
  ggplot2::labs(
    title = "监督学习示例:鸢尾花分类",
    subtitle = "已知类别标签,学习分类边界",
    x = "萼片长度", y = "花瓣长度", color = "种类"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

无监督学习(Unsupervised Learning)

定义:从无标签数据中发现隐藏的结构或模式。

数学表达

给定数据集 \(D = \{x_1, x_2, ..., x_n\}\)(无标签),发现数据的内在结构,如聚类或低维表示。

核心特点: 1. 训练数据没有标签 2. 学习目标是发现数据的内在结构 3. 评估标准不如监督学习明确

# 无监督学习示例:K-means聚类
# 使用鸢尾花数据,但假装不知道真实类别

# 只取特征列(去掉标签)
iris_features <- iris[, 1:4]

# K-means聚类(假设不知道有3类)
set.seed(42)
kmeans_result <- kmeans(iris_features, centers = 3, nstart = 10)

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

# 可视化:聚类结果 vs 真实类别
p1 <- ggplot2::ggplot(iris_clustered, ggplot2::aes(x = Sepal.Length, y = Petal.Length, color = cluster)) +
  ggplot2::geom_point(size = 3, alpha = 0.7) +
  ggplot2::labs(title = "无监督学习:K-means聚类结果", x = "萼片长度", y = "花瓣长度", color = "聚类") +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

p2 <- ggplot2::ggplot(iris, ggplot2::aes(x = Sepal.Length, y = Petal.Length, color = Species)) +
  ggplot2::geom_point(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 = "无监督学习 vs 真实标签对比",
  subtitle = "无监督学习在没有标签的情况下发现数据结构"
)

半监督学习(Semi-supervised Learning)

定义:结合少量有标签数据和大量无标签数据进行学习。

适用场景: 1. 标签获取成本高(需要专家标注) 2. 无标签数据容易获取 3. 医学影像诊断(标注需要专业医生)

核心思想: 1. 利用有标签数据学习初步模型 2. 利用无标签数据增强模型(如伪标签、一致性正则化)

# 半监督学习示意图
# 假设只有10%的数据有标签

set.seed(123)
n_labeled <- 15  # 只有15个样本有标签

# 随机选择有标签的样本
labeled_idx <- sample(1:nrow(iris), n_labeled)
iris_semi <- iris
iris_semi$label_status <- ifelse(1:nrow(iris) %in% labeled_idx, "有标签", "无标签")
iris_semi$known_species <- ifelse(1:nrow(iris) %in% labeled_idx, as.character(iris$Species), "未知")

# 可视化半监督学习场景
ggplot2::ggplot(iris_semi, ggplot2::aes(x = Sepal.Length, y = Petal.Length)) +
  ggplot2::geom_point(ggplot2::aes(color = known_species, shape = label_status), size = 3, alpha = 0.7) +
  ggplot2::scale_shape_manual(values = c("有标签" = 16, "无标签" = 1)) +
  ggplot2::labs(
    title = "半监督学习场景示意",
    subtitle = paste0("只有", n_labeled, "个样本有标签(实心点),其余无标签(空心点)"),
    x = "萼片长度", y = "花瓣长度", color = "已知类别", shape = "标签状态"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

强化学习(Reinforcement Learning)

定义:智能体通过与环境交互,根据奖励/惩罚信号学习最优策略。

核心要素: 1. 智能体(Agent):学习者 2. 环境(Environment):智能体交互的对象 3. 动作(Action):智能体的行为 4. 状态(State):环境的当前情况 5. 奖励(Reward):环境对动作的反馈

数学框架:马尔可夫决策过程(MDP)

\[\text{目标:最大化累积奖励} \quad G_t = \sum_{k=0}^{\infty} \gamma^k R_{t+k+1}\]

其中 \(\gamma\) 是折扣因子,\(0 \leq \gamma < 1\)

# 强化学习概念示意图
# 创建一个简单的网格世界示例

grid_world <- data.frame(
  x = rep(1:4, each = 4),
  y = rep(1:4, 4),
  cell_type = c(
    "起点", "普通", "普通", "普通",
    "普通", "障碍", "普通", "普通",
    "普通", "普通", "普通", "陷阱",
    "普通", "普通", "普通", "终点"
  ),
  reward = c(
    0, -1, -1, -1,
    -1, NA, -1, -1,
    -1, -1, -1, -100,
    -1, -1, -1, 100
  )
)

# 添加标签
grid_world$label <- ifelse(is.na(grid_world$reward), "障碍",
                           ifelse(grid_world$reward == 100, "终点(+100)",
                                  ifelse(grid_world$reward == -100, "陷阱(-100)",
                                         ifelse(grid_world$reward == 0, "起点", ""))))

ggplot2::ggplot(grid_world, ggplot2::aes(x = x, y = y)) +
  ggplot2::geom_tile(ggplot2::aes(fill = cell_type), color = "white", linewidth = 1) +
  ggplot2::geom_text(ggplot2::aes(label = label), size = 3) +
  ggplot2::scale_fill_manual(values = c(
    "起点" = "green", "终点" = "gold", "障碍" = "gray50",
    "陷阱" = "red", "普通" = "lightblue"
  )) +
  ggplot2::labs(
    title = "强化学习示例:网格世界",
    subtitle = "智能体从起点出发,学习避开陷阱、到达终点的最优路径",
    x = "", y = "", fill = "格子类型"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(
    axis.text = ggplot2::element_blank(),
    axis.ticks = ggplot2::element_blank(),
    legend.position = "bottom"
  ) +
  ggplot2::coord_equal()

四种学习范式对比总结

学习类型 数据特点 学习目标 典型算法 医学应用
监督学习 输入-标签对 预测标签 决策树、SVM、神经网络 疾病诊断、预后预测
无监督学习 无标签数据 发现结构 K-means、PCA、自编码器 患者分型、基因聚类
半监督学习 少量标签+大量无标签 利用无标签数据提升性能 伪标签、一致性正则化 影像诊断(标注成本高)
强化学习 环境交互反馈 学习最优策略 Q-learning、Policy Gradient 治疗方案优化

1.4 分类问题、回归问题、聚类问题

三大基本问题类型

机器学习任务根据输出类型的不同,主要分为三大类:

问题类型 输出类型 目标 典型算法
分类问题 离散类别 预测样本所属类别 逻辑回归、决策树、SVM
回归问题 连续数值 预测连续值 线性回归、随机森林回归
聚类问题 无(发现结构) 将相似样本分组 K-means、层次聚类

分类问题(Classification)

定义:预测样本属于哪个离散类别。

数学表达

给定输入 \(x\),预测其类别 \(y \in \{C_1, C_2, ..., C_K\}\)

分类类型: 1. 二分类:只有两个类别(如:患病/健康) 2. 多分类:多个类别(如:肿瘤类型I/II/III期)

# 分类问题示例:乳腺癌诊断
# 使用威斯康星乳腺癌数据集
data(BreastCancer, package = "mlbench")

# 数据预处理
bc_data <- BreastCancer %>%
  dplyr::select(-Id) %>%
  dplyr::mutate_all(~ifelse(. == "?", NA, .)) %>%
  tidyr::drop_na() %>%
  dplyr::mutate_at(vars(-Class), as.numeric)

# 查看数据
head(bc_data[, 1:5])
##   Cl.thickness Cell.size Cell.shape Marg.adhesion Epith.c.size
## 1            5         1          1             1            2
## 2            5         4          4             5            7
## 3            3         1          1             1            2
## 4            6         8          8             1            3
## 5            4         1          1             3            2
## 6            8        10         10             8            7
# 可视化:两类样本在特征空间中的分布
ggplot2::ggplot(bc_data, ggplot2::aes(x = Cl.thickness, y = Cell.size, color = Class)) +
  ggplot2::geom_jitter(size = 2, alpha = 0.5, width = 0.2, height = 0.2) +
  ggplot2::labs(
    title = "分类问题示例:乳腺癌诊断",
    subtitle = "根据细胞特征预测良性/恶性",
    x = "细胞厚度", y = "细胞大小", color = "诊断结果"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

回归问题(Regression)

定义:预测连续的数值输出。

数学表达

给定输入 \(x\),预测连续输出 \(y \in \mathbb{R}\)

与分类的区别: - 分类输出是离散的类别标签 - 回归输出是连续的数值

# 回归问题示例:医疗费用预测
# 使用模拟数据演示

set.seed(456)
n_patients <- 200

# 生成模拟数据
medical_cost <- data.frame(
  age = round(runif(n_patients, 18, 80)),
  bmi = round(rnorm(n_patients, 25, 4), 1),
  smoker = sample(c("Yes", "No"), n_patients, replace = TRUE, prob = c(0.3, 0.7))
) %>%
  dplyr::mutate(
    # 医疗费用与年龄、BMI、吸烟状态相关
    cost = 1000 + 50 * age + 100 * bmi + 
           ifelse(smoker == "Yes", 5000, 0) + 
           rnorm(n_patients, 0, 1000)
  )

# 查看数据
head(medical_cost)
##   age  bmi smoker      cost
## 1  24 25.5     No  5467.349
## 2  31 28.5    Yes  9982.663
## 3  63 24.6     No  6578.476
## 4  71 25.3    Yes 11175.188
## 5  67 18.3     No  6179.478
## 6  39 29.5    Yes 12297.854
# 可视化:费用与年龄的关系
ggplot2::ggplot(medical_cost, ggplot2::aes(x = age, y = cost, color = smoker)) +
  ggplot2::geom_point(size = 2, alpha = 0.6) +
  ggplot2::geom_smooth(method = "lm", se = FALSE) +
  ggplot2::labs(
    title = "回归问题示例:医疗费用预测",
    subtitle = "预测连续的费用数值(元)",
    x = "年龄(岁)", y = "医疗费用(元)", color = "是否吸烟"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

聚类问题(Clustering)

定义:将相似的样本归为同一组,发现数据的内在结构。

特点: 1. 无监督学习(无标签) 2. 目标是发现数据的自然分组 3. “正确”的聚类结果可能不唯一

# 聚类问题示例:患者分型
# 使用模拟数据演示

set.seed(789)
n_patients <- 150

# 生成三种不同类型的患者群体
patient_types <- data.frame(
  # 年轻低风险组
  age = c(rnorm(50, 35, 5), rnorm(50, 55, 8), rnorm(50, 70, 6)),
  risk_score = c(rnorm(50, 20, 5), rnorm(50, 50, 8), rnorm(50, 75, 7)),
  true_group = rep(c("低风险组", "中风险组", "高风险组"), each = 50)
)

# 进行K-means聚类(假装不知道真实分组)
kmeans_patients <- kmeans(patient_types[, c("age", "risk_score")], centers = 3, nstart = 10)
patient_types$cluster <- as.factor(kmeans_patients$cluster)

# 可视化聚类结果
ggplot2::ggplot(patient_types, ggplot2::aes(x = age, y = risk_score)) +
  ggplot2::geom_point(ggplot2::aes(color = cluster), size = 3, alpha = 0.7) +
  ggplot2::labs(
    title = "聚类问题示例:患者风险分型",
    subtitle = "无监督地将患者分为不同风险群体",
    x = "年龄(岁)", y = "风险评分", color = "聚类结果"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

三种问题类型对比

问题类型 输出类型 目标 典型算法 评估指标 医学示例
分类问题 离散类别 预测样本所属类别 逻辑回归、决策树、SVM 准确率、F1分数 疾病诊断、病理分型
回归问题 连续数值 预测连续值 线性回归、随机森林回归 MSE、R² 费用预测、生存时间预测
聚类问题 无(发现结构) 将相似样本分组 K-means、层次聚类 轮廓系数、SSE 患者分型、基因表达聚类

意义与使用场景

核心意义:正确识别问题类型是选择合适算法的前提,不同问题类型对应不同的建模方法和评估指标。

举例:心血管疾病风险预测项目

为什么需要区分问题类型?

同一个医学问题可能被建模为不同类型:预测”是否患病”是分类问题,预测”患病概率”仍是分类(但输出概率),预测”生存年数”是回归问题,发现”患者亚群”是聚类问题。明确问题类型后才能选择正确的方法。

研究目的是什么?

  1. 分类:判断患者是否属于高风险组,用于筛查决策
  2. 回归:预测具体的风险评分或生存时间,用于个体化预后
  3. 聚类:发现具有相似风险特征的患者亚群,用于精准医疗

可能得到什么样的结果?

分类模型可能输出”高风险/低风险”标签,准确率85%;回归模型可能输出”5年生存概率78%“,R²=0.65;聚类模型可能发现3个患者亚群,各有不同的预后特征。

需要注意什么?

有时问题类型可以相互转化:回归问题可以通过设置阈值转为分类问题;分类问题的概率输出可以用于排序。选择哪种类型取决于业务需求和数据特点。


1.5 特征(feature)与标签(label)

核心概念

在监督学习中,数据通常表示为特征-标签对:

概念 定义 符号表示 说明
特征 描述样本的属性 \(X\)\(x\) 输入变量
标签 样本的目标值 \(y\) 输出变量(监督学习特有)

特征(Feature)

定义:描述样本的可测量属性或特性。

特征类型

特征类型 数据类型 示例 编码方式
数值特征 连续/离散 年龄、血压 直接使用
类别特征 名义/有序 性别、疾病分期 独热编码/标签编码
文本特征 字符串 病历描述 TF-IDF/词嵌入
图像特征 像素矩阵 X光片 CNN提取
# 特征示例:心脏病预测数据集
# 创建模拟的心脏病预测数据

set.seed(111)
n <- 100

heart_data <- data.frame(
  # 数值特征
  age = round(runif(n, 30, 80)),                    # 年龄(连续)
  resting_bp = round(rnorm(n, 130, 20)),            # 静息血压(连续)
  cholesterol = round(rnorm(n, 200, 40)),           # 胆固醇(连续)
  max_heart_rate = round(220 - runif(n, 30, 80) + rnorm(n, 0, 10)),  # 最大心率
  
  # 类别特征
  sex = sample(c("Male", "Female"), n, replace = TRUE),           # 性别(名义)
  chest_pain = sample(c("Typical", "Atypical", "Non-anginal", "Asymptomatic"), 
                      n, replace = TRUE),                          # 胸痛类型(名义)
  fasting_blood_sugar = sample(c(">120", "<=120"), n, replace = TRUE),  # 空腹血糖(二值)
  
  # 标签
  heart_disease = sample(c("Yes", "No"), n, replace = TRUE, prob = c(0.4, 0.6))
)

# 查看数据结构
head(heart_data)
##   age resting_bp cholesterol max_heart_rate    sex   chest_pain
## 1  60        134         196            183   Male Asymptomatic
## 2  66        161         125            153   Male  Non-anginal
## 3  49        148         173            147 Female      Typical
## 4  56        137         208            157   Male Asymptomatic
## 5  49        134          96            162   Male Asymptomatic
## 6  51        113         196            154 Female  Non-anginal
##   fasting_blood_sugar heart_disease
## 1                >120            No
## 2                >120            No
## 3                >120            No
## 4                >120           Yes
## 5                >120            No
## 6               <=120            No
# 特征类型总结
特征名称 特征类型 数据类型 取值范围
age 数值-连续 整数 30-80岁
resting_bp 数值-连续 整数 约80-180 mmHg
cholesterol 数值-连续 整数 约100-300 mg/dL
sex 类别-名义 字符 Male/Female
chest_pain 类别-名义 字符 4种类型
fasting_blood_sugar 类别-二值 字符 >120/<=120

标签(Label)

定义:监督学习中样本的目标值或正确答案。

标签类型

任务类型 标签类型 示例
二分类 二值 患病/健康
多分类 多个离散值 肿瘤I/II/III期
回归 连续数值 医疗费用、生存时间
# 标签分布可视化
ggplot2::ggplot(heart_data, ggplot2::aes(x = heart_disease, fill = heart_disease)) +
  ggplot2::geom_bar(width = 0.6) +
  ggplot2::geom_text(stat = "count", ggplot2::aes(label = after_stat(count)), vjust = -0.5) +
  ggplot2::labs(
    title = "标签分布:心脏病诊断结果",
    x = "诊断结果", y = "样本数量", fill = "诊断结果"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "none")

特征矩阵与标签向量

在机器学习中,数据通常表示为:

  • 特征矩阵 \(X\)\(n \times p\) 矩阵,\(n\) 个样本,\(p\) 个特征
  • 标签向量 \(y\)\(n \times 1\) 向量

\[X = \begin{bmatrix} x_{11} & x_{12} & \cdots & x_{1p} \\ x_{21} & x_{22} & \cdots & x_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ x_{n1} & x_{n2} & \cdots & x_{np} \end{bmatrix}, \quad y = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}\]

# 特征矩阵与标签向量的表示
# 提取数值特征作为特征矩阵
X <- heart_data %>%
  dplyr::select(age, resting_bp, cholesterol, max_heart_rate) %>%
  as.matrix()

# 提取标签向量
y <- heart_data$heart_disease

# 显示维度
cat("特征矩阵 X 的维度:", dim(X), "\n")
## 特征矩阵 X 的维度: 100 4
cat("  - 样本数量 n =", nrow(X), "\n")
##   - 样本数量 n = 100
cat("  - 特征数量 p =", ncol(X), "\n")
##   - 特征数量 p = 4
cat("\n标签向量 y 的长度:", length(y), "\n")
## 
## 标签向量 y 的长度: 100
cat("标签类别:", unique(y), "\n")
## 标签类别: No Yes
# 特征矩阵前几行
cat("\n特征矩阵 X 前5行:\n")
## 
## 特征矩阵 X 前5行:
print(head(X, 5))
##      age resting_bp cholesterol max_heart_rate
## [1,]  60        134         196            183
## [2,]  66        161         125            153
## [3,]  49        148         173            147
## [4,]  56        137         208            157
## [5,]  49        134          96            162

意义与使用场景

核心意义:特征是机器学习模型的”原材料”,标签是监督学习的”答案”。特征的质量直接决定模型的上限。

举例:基于电子病历预测患者再入院风险

为什么特征选择很重要?

再入院风险受多种因素影响:年龄、诊断、用药、既往史、社会因素等。选择哪些作为特征,决定了模型能学到什么。如果遗漏重要特征(如既往再入院史),模型性能会受限;如果包含无关特征(如病历编号),可能引入噪声。

研究目的是什么?

  1. 识别高风险患者:在出院前预测30天内再入院概率
  2. 指导干预措施:对高风险患者进行重点随访
  3. 资源配置:合理安排出院计划和随访资源

可能得到什么样的结果?

模型可能发现:年龄大于65岁、有多次既往住院史、出院时用药超过5种的患者,再入院风险显著升高。这些特征组合可以生成风险评分。

需要注意什么?

特征工程是机器学习中最耗时的环节之一。需要注意:1. 特征的完整性和准确性;2. 特征之间的相关性;3. 特征的可解释性;4. 避免数据泄露(如使用了未来信息)。


1.6 训练集、验证集、测试集划分

三种数据集的作用

在机器学习中,数据通常划分为三个独立的部分:

数据集 作用 比例建议 使用场景
训练集 训练模型参数 60-80% 模型学习
验证集 调整超参数、模型选择 10-20% 模型调优
测试集 评估最终性能 10-20% 最终评估

为什么要划分?

核心问题:模型在训练数据上的表现好,不代表在新数据上也好。

目标:评估模型的泛化能力——在未见数据上的表现。

三种数据集的分工

  1. 训练集:用于学习模型参数(如权重)
  2. 验证集:用于选择模型、调整超参数(如正则化系数)
  3. 测试集:最终评估,只能使用一次!
# 数据划分示例
# 使用iris数据集演示

set.seed(42)

# 获取样本数量
n <- nrow(iris)

# 随机打乱样本顺序
indices <- sample(1:n)

# 计算划分点
train_end <- round(0.7 * n)
val_end <- round(0.85 * n)

# 划分数据集
train_idx <- indices[1:train_end]
val_idx <- indices[(train_end + 1):val_end]
test_idx <- indices[(val_end + 1):n]

# 创建三个数据集
train_data <- iris[train_idx, ]
val_data <- iris[val_idx, ]
test_data <- iris[test_idx, ]

# 显示划分结果
cat("训练集样本数:", nrow(train_data), "\n")
## 训练集样本数: 105
cat("验证集样本数:", nrow(val_data), "\n")
## 验证集样本数: 23
cat("测试集样本数:", nrow(test_data), "\n")
## 测试集样本数: 22
# 可视化数据划分
split_viz <- data.frame(
  dataset = rep(c("训练集", "验证集", "测试集"), 
                times = c(nrow(train_data), nrow(val_data), nrow(test_data))),
  index = c(train_idx, val_idx, test_idx)
)

split_viz$dataset <- factor(split_viz$dataset, levels = c("训练集", "验证集", "测试集"))

ggplot2::ggplot(split_viz, ggplot2::aes(x = index, y = 1, color = dataset)) +
  ggplot2::geom_point(size = 2, alpha = 0.7) +
  ggplot2::labs(
    title = "数据集划分可视化",
    subtitle = "随机打乱后按比例划分",
    x = "样本索引", y = "", color = "数据集"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(
    axis.text.y = ggplot2::element_blank(),
    axis.ticks.y = ggplot2::element_blank(),
    legend.position = "bottom"
  )

使用caret包进行数据划分

# 使用caret包进行更规范的数据划分
library(caret)

# 创建训练集和测试集(70%-30%划分)
set.seed(42)
train_index <- caret::createDataPartition(iris$Species, p = 0.7, list = FALSE)

train_set <- iris[train_index, ]
test_set <- iris[-train_index, ]

# 在训练集内部再划分出验证集
set.seed(42)
val_index <- caret::createDataPartition(train_set$Species, p = 0.8, list = FALSE)

train_final <- train_set[val_index, ]
val_set <- train_set[-val_index, ]

# 显示最终划分
cat("训练集样本数:", nrow(train_final), "\n")
## 训练集样本数: 84
cat("验证集样本数:", nrow(val_set), "\n")
## 验证集样本数: 21
cat("测试集样本数:", nrow(test_set), "\n")
## 测试集样本数: 45

分层抽样

当类别不平衡时,需要使用分层抽样确保各数据集中类别比例一致。

# 分层抽样示例
# 创建不平衡数据集

set.seed(123)
imbalanced_data <- data.frame(
  feature = rnorm(1000),
  label = c(rep("Positive", 100), rep("Negative", 900))  # 1:9不平衡
)

# 普通随机划分 vs 分层划分
set.seed(42)

# 普通随机划分
random_split <- sample(1:nrow(imbalanced_data), 0.7 * nrow(imbalanced_data))
train_random <- imbalanced_data[random_split, ]

# 分层划分(使用caret)
stratified_index <- caret::createDataPartition(imbalanced_data$label, p = 0.7, list = FALSE)
train_stratified <- imbalanced_data[stratified_index, ]

# 比较两种划分方式的类别比例
cat("原始数据Positive比例:", sum(imbalanced_data$label == "Positive") / nrow(imbalanced_data) * 100, "%\n")
## 原始数据Positive比例: 10 %
cat("随机划分Positive比例:", sum(train_random$label == "Positive") / nrow(train_random) * 100, "%\n")
## 随机划分Positive比例: 10.14286 %
cat("分层划分Positive比例:", sum(train_stratified$label == "Positive") / nrow(train_stratified) * 100, "%\n")
## 分层划分Positive比例: 10 %

意义与使用场景

核心意义:合理的数据划分是评估模型泛化能力的基础,避免”用考试题训练”的作弊行为。

举例:构建糖尿病预测模型

为什么要划分三个数据集?

假设我们有1000名患者的数据。如果全部用于训练,模型可能”记住”了这些患者的特征,在新患者上表现很差。划分后:700人用于训练,150人用于调参,150人用于最终测试。测试集的150人模拟”新患者”,评估模型的真实泛化能力。

研究目的是什么?

  1. 训练集:让模型学习糖尿病的预测规律
  2. 验证集:选择最佳模型参数(如决策树深度)
  3. 测试集:报告模型的最终性能

可能得到什么样的结果?

训练集准确率85%,验证集准确率82%,测试集准确率80%。如果测试集准确率远低于训练集,说明模型过拟合,需要简化模型或增加数据。

需要注意什么?

  1. 测试集只能使用一次,不能反复调整模型后再测试
  2. 时间序列数据不能随机划分,要按时间顺序
  3. 医学数据要考虑患者级别的划分(同一患者的数据不能同时出现在训练集和测试集)

1.7 过拟合与欠拟合

核心概念

过拟合(Overfitting):模型在训练数据上表现很好,但在新数据上表现差。模型”记住”了训练数据的噪声和细节,而非真正的规律。

欠拟合(Underfitting):模型在训练数据和新数据上都表现不好。模型太简单,无法捕捉数据中的规律。

# 过拟合与欠拟合示意图
# 使用多项式回归演示

set.seed(42)
n <- 30

# 生成真实数据(带噪声的二次函数)
x <- seq(-3, 3, length.out = n)
y_true <- 0.5 * x^2 - 1
y <- y_true + rnorm(n, sd = 0.5)

# 创建数据框
poly_data <- data.frame(x = x, y = y, y_true = y_true)

# 拟合不同复杂度的模型
# 欠拟合:线性模型(太简单)
model_underfit <- lm(y ~ x, data = poly_data)

# 适度拟合:二次多项式
model_good <- lm(y ~ poly(x, degree = 2, raw = TRUE), data = poly_data)

# 过拟合:高阶多项式(太复杂)
model_overfit <- lm(y ~ poly(x, degree = 15, raw = TRUE), data = poly_data)

# 预测
poly_data$pred_underfit <- predict(model_underfit)
poly_data$pred_good <- predict(model_good)
poly_data$pred_overfit <- predict(model_overfit)

# 可视化
p1 <- ggplot2::ggplot(poly_data, ggplot2::aes(x = x, y = y)) +
  ggplot2::geom_point(size = 2, alpha = 0.7) +
  ggplot2::geom_line(ggplot2::aes(y = pred_underfit), color = "blue", linewidth = 1) +
  ggplot2::labs(title = "欠拟合(线性模型)", subtitle = "模型太简单,无法捕捉规律", x = "", y = "") +
  ggplot2::theme_minimal() +
  ggplot2::ylim(-3, 5)

p2 <- ggplot2::ggplot(poly_data, ggplot2::aes(x = x, y = y)) +
  ggplot2::geom_point(size = 2, alpha = 0.7) +
  ggplot2::geom_line(ggplot2::aes(y = pred_good), color = "green", linewidth = 1) +
  ggplot2::labs(title = "适度拟合(二次多项式)", subtitle = "模型复杂度适中", x = "", y = "") +
  ggplot2::theme_minimal() +
  ggplot2::ylim(-3, 5)

p3 <- ggplot2::ggplot(poly_data, ggplot2::aes(x = x, y = y)) +
  ggplot2::geom_point(size = 2, alpha = 0.7) +
  ggplot2::geom_line(ggplot2::aes(y = pred_overfit), color = "red", linewidth = 1) +
  ggplot2::labs(title = "过拟合(15阶多项式)", subtitle = "模型太复杂,拟合了噪声", x = "", y = "") +
  ggplot2::theme_minimal() +
  ggplot2::ylim(-3, 5)

p1 + p2 + p3 + patchwork::plot_annotation(
  title = "欠拟合 vs 适度拟合 vs 过拟合",
  subtitle = "蓝点为训练数据,曲线为模型预测"
)

过拟合与欠拟合的特征

特征 欠拟合 过拟合
训练误差
测试误差
原因 模型太简单 模型太复杂
表现 学不到规律 记住噪声
解决方法 增加模型复杂度 正则化、增加数据、简化模型
# 计算训练误差和测试误差
# 生成新数据作为测试集
set.seed(999)
x_test <- seq(-3, 3, length.out = 100)
y_test_true <- 0.5 * x_test^2 - 1
y_test <- y_test_true + rnorm(100, sd = 0.5)
test_data <- data.frame(x = x_test, y = y_test)

# 预测测试集
test_data$pred_underfit <- predict(model_underfit, newdata = test_data)
test_data$pred_good <- predict(model_good, newdata = test_data)
test_data$pred_overfit <- predict(model_overfit, newdata = test_data)

# 计算均方误差
calc_mse <- function(actual, predicted) {
  mean((actual - predicted)^2)
}

cat("欠拟合(线性)- 训练MSE:", calc_mse(poly_data$y, poly_data$pred_underfit), "测试MSE:", calc_mse(test_data$y, test_data$pred_underfit), "\n")
## 欠拟合(线性)- 训练MSE: 2.525686 测试MSE: 2.151376
cat("适度拟合(二次)- 训练MSE:", calc_mse(poly_data$y, poly_data$pred_good), "测试MSE:", calc_mse(test_data$y, test_data$pred_good), "\n")
## 适度拟合(二次)- 训练MSE: 0.3484976 测试MSE: 0.2811028
cat("过拟合(15阶)- 训练MSE:", calc_mse(poly_data$y, poly_data$pred_overfit), "测试MSE:", calc_mse(test_data$y, test_data$pred_overfit), "\n")
## 过拟合(15阶)- 训练MSE: 0.2188181 测试MSE: 0.4459991

如何判断过拟合/欠拟合?

# 模型复杂度与误差的关系图
# 模拟不同复杂度下的训练误差和测试误差

complexity <- 1:15
train_errors <- c()
test_errors <- c()

for (deg in complexity) {
  model <- lm(y ~ poly(x, degree = deg, raw = TRUE), data = poly_data)
  train_errors <- c(train_errors, calc_mse(poly_data$y, predict(model)))
  test_errors <- c(test_errors, calc_mse(test_data$y, predict(model, newdata = test_data)))
}

error_data <- data.frame(
  complexity = rep(complexity, 2),
  error = c(train_errors, test_errors),
  type = rep(c("训练误差", "测试误差"), each = length(complexity))
)

ggplot2::ggplot(error_data, ggplot2::aes(x = complexity, y = error, color = type)) +
  ggplot2::geom_line(linewidth = 1) +
  ggplot2::geom_point(size = 2) +
  ggplot2::geom_vline(xintercept = 2, linetype = "dashed", color = "gray50") +
  ggplot2::annotate("text", x = 2.5, y = max(error_data$error) * 0.9, 
                    label = "最优复杂度", hjust = 0) +
  ggplot2::labs(
    title = "模型复杂度与误差的关系",
    subtitle = "训练误差持续下降,测试误差先降后升(过拟合)",
    x = "模型复杂度(多项式阶数)", y = "均方误差(MSE)", color = ""
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

解决方法

问题 解决方法
欠拟合 增加特征、使用更复杂模型、减少正则化
过拟合 增加训练数据、正则化、简化模型、早停、Dropout

意义与使用场景

核心意义:过拟合和欠拟合是机器学习中最常见的问题,找到合适的模型复杂度是建模的关键。

举例:构建心脏病预测模型

过拟合的场景

使用包含100个特征的深度神经网络,训练集准确率99%,但测试集只有65%。模型”记住”了训练集中每个患者的特殊情况,包括噪声和巧合。

欠拟合的场景

使用简单的逻辑回归,只用年龄一个特征,训练集和测试集准确率都只有60%。模型太简单,无法捕捉疾病与多种因素之间的复杂关系。

如何解决?

过拟合:增加数据、使用正则化、简化模型(如减少层数)、特征选择。欠拟合:增加特征、使用更复杂的模型(如随机森林)、减少正则化强度。


1.8 偏差-方差权衡(bias-variance tradeoff)

数学基础

期望预测误差可以分解为三部分:

\[E[(y - \hat{f}(x))^2] = \underbrace{Bias[\hat{f}(x)]^2}_{偏差} + \underbrace{Var[\hat{f}(x)]}_{方差} + \underbrace{\sigma^2}_{不可约误差}\]

其中: - 偏差(Bias):模型预测的期望值与真实值的差异 \[Bias[\hat{f}(x)] = E[\hat{f}(x)] - f(x)\] - 方差(Variance):模型预测的波动程度 \[Var[\hat{f}(x)] = E[(\hat{f}(x) - E[\hat{f}(x)])^2]\] - 不可约误差(Irreducible Error):数据本身的噪声,无法消除

偏差与方差的直观理解

概念 含义 高偏差表现 高方差表现
偏差 预测的准确性 预测值偏离真实值 -
方差 预测的稳定性 - 不同训练集得到差异大的模型

靶心类比

  • 低偏差低方差:子弹集中在靶心(理想状态)
  • 高偏差低方差:子弹集中但偏离靶心(欠拟合)
  • 低偏差高方差:子弹分散但平均在靶心附近(过拟合)
  • 高偏差高方差:子弹分散且偏离靶心(最差情况)
# 偏差-方差权衡示意图
# 使用靶心图演示

set.seed(42)

# 创建四种情况的模拟数据
n_points <- 50

# 低偏差低方差
low_bias_low_var <- data.frame(
  x = rnorm(n_points, mean = 0, sd = 0.3),
  y = rnorm(n_points, mean = 0, sd = 0.3),
  type = "低偏差 低方差\n(理想)"
)

# 高偏差低方差
high_bias_low_var <- data.frame(
  x = rnorm(n_points, mean = 2, sd = 0.3),
  y = rnorm(n_points, mean = 2, sd = 0.3),
  type = "高偏差 低方差\n(欠拟合)"
)

# 低偏差高方差
low_bias_high_var <- data.frame(
  x = rnorm(n_points, mean = 0, sd = 1.5),
  y = rnorm(n_points, mean = 0, sd = 1.5),
  type = "低偏差 高方差\n(过拟合)"
)

# 高偏差高方差
high_bias_high_var <- data.frame(
  x = rnorm(n_points, mean = 2, sd = 1.5),
  y = rnorm(n_points, mean = 2, sd = 1.5),
  type = "高偏差 高方差\n(最差)"
)

# 合并数据
target_data <- rbind(low_bias_low_var, high_bias_low_var, 
                     low_bias_high_var, high_bias_high_var)
target_data$type <- factor(target_data$type, levels = unique(target_data$type))

# 绘制靶心图
ggplot2::ggplot(target_data, ggplot2::aes(x = x, y = y)) +
  ggplot2::geom_point(alpha = 0.6, size = 2) +
  ggplot2::geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  ggplot2::geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
  ggplot2::annotate("point", x = 0, y = 0, shape = 21, size = 5, 
                    fill = "red", color = "black", stroke = 1) +
  ggplot2::facet_wrap(~ type) +
  ggplot2::labs(
    title = "偏差-方差权衡:靶心图",
    subtitle = "红点为靶心(真实值),蓝点为预测值",
    x = "", y = ""
  ) +
  ggplot2::theme_minimal() +
  ggplot2::coord_fixed(xlim = c(-4, 4), ylim = c(-4, 4))

模型复杂度与偏差-方差的关系

# 模型复杂度与偏差/方差的关系
complexity <- 1:20
bias_squared <- 1 / complexity + 0.1  # 偏差随复杂度降低
variance <- complexity / 20  # 方差随复杂度增加
total_error <- bias_squared + variance + 0.1  # 总误差

bv_data <- data.frame(
  complexity = rep(complexity, 3),
  value = c(bias_squared, variance, total_error),
  component = rep(c("偏差²", "方差", "总误差"), each = length(complexity))
)

ggplot2::ggplot(bv_data, ggplot2::aes(x = complexity, y = value, color = component)) +
  ggplot2::geom_line(linewidth = 1) +
  ggplot2::geom_point(size = 2) +
  ggplot2::geom_vline(xintercept = which.min(total_error), linetype = "dashed", color = "gray50") +
  ggplot2::annotate("text", x = which.min(total_error) + 0.5, y = max(bv_data$value) * 0.9,
                    label = "最优复杂度", hjust = 0) +
  ggplot2::labs(
    title = "偏差-方差权衡",
    subtitle = "总误差 = 偏差² + 方差 + 不可约误差",
    x = "模型复杂度", y = "误差", color = "组成部分"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

不同算法的偏差-方差特点

算法类型 偏差 方差 特点
线性回归 简单稳定,但可能欠拟合
决策树(深) 复杂灵活,但容易过拟合
随机森林 集成降低方差
Boosting 集成降低偏差

意义与使用场景

核心意义:偏差-方差权衡是机器学习的核心矛盾,需要在模型复杂度上找到平衡点。

举例:选择疾病预测模型

高偏差模型(如简单逻辑回归)

只用年龄预测糖尿病风险。模型稳定(低方差),但预测不够准确(高偏差),因为忽略了BMI、血糖等重要因素。

高方差模型(如深度决策树)

使用包含100个特征的深度决策树。在训练集上完美预测,但换一批患者后预测结果变化很大(高方差),因为模型过度拟合了训练数据的特殊模式。

如何权衡?

选择随机森林或梯度提升树,通过集成学习在偏差和方差之间取得平衡。同时使用交叉验证选择最优超参数。


1.9 模型泛化能力

核心概念

泛化能力(Generalization):模型在未见数据上的表现能力。

泛化误差(Generalization Error):模型在新数据上的期望误差。

\[\text{泛化误差} = E_{(x,y) \sim P} [L(y, \hat{f}(x))]\]

其中 \(P\) 是数据的真实分布,\(L\) 是损失函数。

泛化的核心问题

训练误差 vs 泛化误差

  • 训练误差:模型在训练数据上的误差
  • 泛化误差:模型在所有可能数据上的期望误差

目标:最小化泛化误差,而非训练误差。

# 泛化能力示意图
# 使用不同大小的训练集演示

set.seed(42)

# 生成数据
n_total <- 1000
x_all <- runif(n_total, 0, 10)
y_all <- 2 * x_all + 1 + rnorm(n_total, sd = 2)  # 真实关系: y = 2x + 1 + 噪声
all_data <- data.frame(x = x_all, y = y_all)

# 测试集(固定)
test_idx <- sample(1:n_total, 200)
test_data <- all_data[test_idx, ]
train_pool <- all_data[-test_idx, ]

# 不同训练集大小下的训练误差和测试误差
train_sizes <- c(20, 50, 100, 200, 400, 600, 800)
results <- data.frame()

for (n_train in train_sizes) {
  # 多次实验取平均
  train_errors <- c()
  test_errors <- c()
  
  for (i in 1:10) {
    idx <- sample(1:nrow(train_pool), n_train)
    train_data <- train_pool[idx, ]
    
    model <- lm(y ~ x, data = train_data)
    
    train_errors <- c(train_errors, mean((train_data$y - predict(model, train_data))^2))
    test_errors <- c(test_errors, mean((test_data$y - predict(model, test_data))^2))
  }
  
  results <- rbind(results, data.frame(
    train_size = n_train,
    train_error = mean(train_errors),
    test_error = mean(test_errors)
  ))
}

# 可视化
results_long <- tidyr::pivot_longer(results, cols = c(train_error, test_error),
                                     names_to = "type", values_to = "error")

ggplot2::ggplot(results_long, ggplot2::aes(x = train_size, y = error, color = type)) +
  ggplot2::geom_line(linewidth = 1) +
  ggplot2::geom_point(size = 2) +
  ggplot2::labs(
    title = "训练数据量与泛化能力的关系",
    subtitle = "随着训练数据增加,训练误差和测试误差趋于接近",
    x = "训练样本数量", y = "均方误差(MSE)", color = ""
  ) +
  ggplot2::scale_color_discrete(labels = c("训练误差", "测试误差(泛化误差估计)")) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

提高泛化能力的方法

方法 原理 适用场景
增加训练数据 减少过拟合 数据不足时
正则化 限制模型复杂度 过拟合时
交叉验证 更准确估计泛化误差 模型选择时
特征选择 减少噪声特征 特征过多时
集成学习 降低方差 单模型不稳定时

泛化误差界

根据统计学习理论,泛化误差有理论上界:

\[\text{泛化误差} \leq \text{训练误差} + \sqrt{\frac{d \log(n/d) + \log(1/\delta)}{n}}\]

其中: - \(d\):模型复杂度(如VC维) - \(n\):训练样本数 - \(\delta\):置信度

启示: 1. 模型越复杂,泛化误差上界越大 2. 训练数据越多,泛化误差上界越小

意义与使用场景

核心意义:泛化能力是机器学习模型的终极目标,模型的价值在于对未知数据的预测能力。

举例:医学诊断模型的临床应用

训练场景:在A医院的历史数据上训练肺癌诊断模型,训练集准确率95%。

泛化问题:当模型部署到B医院时,准确率可能下降到80%。原因包括:A、B医院的患者人群不同、设备差异、诊断标准差异等。

如何提高泛化能力?

  1. 多中心数据训练:使用多家医院的数据
  2. 数据增强:增加数据的多样性
  3. 正则化:防止模型过度拟合特定医院的特点
  4. 域适应:针对新医院的数据进行调整

1.10 学习曲线(learning curve)

核心概念

学习曲线:展示模型性能随训练数据量变化的曲线。

用途: 1. 诊断模型是否过拟合/欠拟合 2. 判断增加数据是否有帮助 3. 确定最优训练数据量

学习曲线的典型形态

# 学习曲线示例
# 使用caret包绘制学习曲线

set.seed(42)

# 加载数据
data(Sonar, package = "mlbench")

# 定义训练控制
train_control <- caret::trainControl(
  method = "cv",
  number = 5
)

# 计算不同训练集大小下的性能
train_sizes <- seq(0.1, 1, by = 0.1)
n <- nrow(Sonar)

results <- data.frame()

for (prop in train_sizes) {
  n_train <- round(prop * n)
  
  # 多次采样取平均
  accs <- c()
  for (i in 1:5) {
    idx <- sample(1:n, n_train)
    train_subset <- Sonar[idx, ]
    
    model <- caret::train(
      Class ~ .,
      data = train_subset,
      method = "rpart",
      trControl = train_control,
      tuneLength = 3
    )
    
    accs <- c(accs, max(model$results$Accuracy))
  }
  
  results <- rbind(results, data.frame(
    train_size = n_train,
    accuracy = mean(accs),
    sd = sd(accs)
  ))
}

# 绘制学习曲线
ggplot2::ggplot(results, ggplot2::aes(x = train_size, y = accuracy)) +
  ggplot2::geom_line(color = "steelblue", linewidth = 1) +
  ggplot2::geom_point(color = "steelblue", size = 2) +
  ggplot2::geom_ribbon(ggplot2::aes(ymin = accuracy - sd, ymax = accuracy + sd),
                       alpha = 0.2, fill = "steelblue") +
  ggplot2::labs(
    title = "学习曲线示例",
    subtitle = "模型性能随训练数据量增加而提升",
    x = "训练样本数量", y = "交叉验证准确率"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::ylim(0.5, 1)

不同情况下的学习曲线

# 模拟不同情况的学习曲线

train_sizes <- 1:100

# 欠拟合:性能低且平稳
underfit_curve <- 0.6 + 0.05 * (1 - exp(-train_sizes / 20))

# 过拟合:训练性能高,验证性能低
overfit_train <- 0.95 + 0.04 * (1 - exp(-train_sizes / 30))
overfit_val <- 0.5 + 0.3 * (1 - exp(-train_sizes / 50))

# 适度拟合:两条曲线接近
good_fit_train <- 0.85 + 0.1 * (1 - exp(-train_sizes / 20))
good_fit_val <- 0.8 + 0.12 * (1 - exp(-train_sizes / 25))

# 创建数据框
lc_data <- data.frame(
  train_size = rep(train_sizes, 6),
  accuracy = c(underfit_curve, underfit_curve,
               overfit_train, overfit_val,
               good_fit_train, good_fit_val),
  type = rep(c("欠拟合", "过拟合", "适度拟合"), each = 200),
  curve = rep(rep(c("训练集", "验证集"), each = 100), 3)
)

lc_data$type <- factor(lc_data$type, levels = c("欠拟合", "适度拟合", "过拟合"))

# 绘制
ggplot2::ggplot(lc_data, ggplot2::aes(x = train_size, y = accuracy, color = curve)) +
  ggplot2::geom_line(linewidth = 1) +
  ggplot2::facet_wrap(~ type) +
  ggplot2::labs(
    title = "不同情况下的学习曲线",
    subtitle = "训练集曲线(红)vs 验证集曲线(蓝)",
    x = "训练样本数量", y = "准确率", color = ""
  ) +
  ggplot2::scale_color_manual(values = c("训练集" = "red", "验证集" = "blue")) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

如何解读学习曲线?

曲线特征 诊断 建议
训练和验证曲线都低且接近 欠拟合 增加模型复杂度
训练曲线高,验证曲线低 过拟合 增加数据、正则化
曲线仍在上升 数据不足 继续收集数据
曲线已平稳 数据充足 无需更多数据

1.11 交叉验证(K折、留一、分层K折)

核心概念

交叉验证(Cross-Validation):将数据分成多份,轮流用其中一份做验证、其余做训练,综合评估模型性能。

目的:更充分地利用数据,获得更稳定的性能估计。

K折交叉验证(K-Fold CV)

方法:将数据分成K份,进行K次训练和验证,每次用其中1份验证、其余K-1份训练。

数学表达

\[CV_{(K)} = \frac{1}{K} \sum_{i=1}^{K} E_i\]

其中 \(E_i\) 是第 \(i\) 折的验证误差。

# K折交叉验证示意图
# 使用5折交叉验证演示

set.seed(42)
n <- 50
data <- data.frame(
  id = 1:n,
  value = rnorm(n)
)

# 创建5折划分
k <- 5
folds <- sample(rep(1:k, length.out = n))

# 可视化交叉验证过程
cv_viz <- data.frame(
  id = rep(1:n, k),
  fold = rep(1:k, each = n),
  role = rep(folds, k)
)

for (i in 1:k) {
  cv_viz$role[(i-1)*n + 1:n] <- ifelse(folds == i, "验证集", "训练集")
}

cv_viz$fold <- factor(cv_viz$fold)

ggplot2::ggplot(cv_viz, ggplot2::aes(x = id, y = fold, fill = role)) +
  ggplot2::geom_tile(color = "white", linewidth = 0.5) +
  ggplot2::scale_fill_manual(values = c("训练集" = "lightblue", "验证集" = "salmon")) +
  ggplot2::labs(
    title = "5折交叉验证示意图",
    subtitle = "每一行代表一次验证,红色为验证集,蓝色为训练集",
    x = "样本编号", y = "折数", fill = "角色"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

使用caret包进行交叉验证

# 使用caret包进行K折交叉验证
set.seed(42)

# 定义交叉验证控制
cv_control <- caret::trainControl(
  method = "cv",        # K折交叉验证
  number = 5,           # 5折
  verboseIter = FALSE   # 不显示每次迭代的进度
)

# 训练模型(使用5折交叉验证)
model_cv <- caret::train(
  Species ~ .,
  data = iris,
  method = "rpart",
  trControl = cv_control,
  tuneLength = 3
)

# 查看交叉验证结果
print(model_cv)
## CART 
## 
## 150 samples
##   4 predictor
##   3 classes: 'setosa', 'versicolor', 'virginica' 
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 120, 120, 120, 120, 120 
## Resampling results across tuning parameters:
## 
##   cp    Accuracy   Kappa
##   0.00  0.9400000  0.91 
##   0.44  0.7600000  0.64 
##   0.50  0.3333333  0.00 
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was cp = 0.

留一交叉验证(Leave-One-Out CV, LOOCV)

方法:每次留一个样本做验证,其余全部做训练。适用于小样本数据。

特点: 1. 估计最准确(偏差最小) 2. 计算成本高(需要训练n次) 3. 方差可能较大

# LOOCV示例
# 使用caret包
set.seed(42)

loocv_control <- caret::trainControl(
  method = "LOOCV"  # 留一交叉验证
)

# 训练模型(使用LOOCV)
model_loocv <- caret::train(
  Species ~ .,
  data = iris,
  method = "rpart",
  trControl = loocv_control
)

# 查看结果
cat("LOOCV准确率:", model_loocv$results$Accuracy, "\n")
## LOOCV准确率: 0.9333333 0.3333333 0

分层K折交叉验证(Stratified K-Fold)

方法:确保每折中各类别的比例与整体数据一致。

适用场景:类别不平衡的分类问题。

# 分层K折交叉验证
set.seed(42)

# 创建不平衡数据集
imbalanced_iris <- iris %>%
  dplyr::filter(!(Species == "setosa" & row_number() > 10))  # setosa只保留10个

# 普通K折 vs 分层K折
# 使用caret的createFolds函数
normal_folds <- caret::createFolds(imbalanced_iris$Species, k = 5)
stratified_folds <- caret::createFolds(imbalanced_iris$Species, k = 5, returnTrain = TRUE)

# 比较每折中各类别的比例
compare_folds <- data.frame()
for (i in 1:5) {
  test_idx <- normal_folds[[i]]
  compare_folds <- rbind(compare_folds, data.frame(
    折数 = i,
    方法 = "普通K折",
    setosa比例 = mean(imbalanced_iris$Species[test_idx] == "setosa")
  ))
}

# 分层K折的比例应该更一致
for (i in 1:5) {
  test_idx <- normal_folds[[i]]
  cat("折", i, "setosa比例:", mean(imbalanced_iris$Species[test_idx] == "setosa"), "\n")
}
## 折 1 setosa比例: 0.09090909 
## 折 2 setosa比例: 0.09090909 
## 折 3 setosa比例: 0.09090909 
## 折 4 setosa比例: 0.09090909 
## 折 5 setosa比例: 0.09090909

不同交叉验证方法对比

方法 训练次数 偏差 方差 适用场景
K折(K=5) 5 一般场景
K折(K=10) 10 中高 常用选择
LOOCV n 最低 小样本
分层K折 K 类别不平衡

意义与使用场景

核心意义:交叉验证是评估模型泛化能力的标准方法,比单次划分更稳定可靠。

举例:医学研究中评估疾病预测模型

为什么需要交叉验证?

医学数据通常样本量有限(如200例患者)。如果只做一次训练-测试划分,测试集可能恰好包含”简单”或”困难”的样本,导致评估结果不稳定。交叉验证让每个样本都有机会被验证,评估更全面。

如何选择交叉验证方法?

  1. 样本量充足(>1000):5折或10折交叉验证
  2. 样本量中等(100-1000):10折交叉验证
  3. 样本量小(<100):留一交叉验证
  4. 类别不平衡:分层K折交叉验证

1.12 评价指标概述

为什么需要多种评价指标?

不同任务、不同场景需要关注不同的性能维度:

场景 关注重点 适合指标
疾病筛查 不漏诊 召回率
确诊检查 不误诊 精确率
综合评估 平衡 F1分数
模型比较 整体性能 AUC

分类问题常用指标

指标 公式 含义
准确率 \(\frac{TP+TN}{TP+TN+FP+FN}\) 预测正确的比例
精确率 \(\frac{TP}{TP+FP}\) 预测为正的样本中真正为正的比例
召回率 \(\frac{TP}{TP+FN}\) 真正为正的样本被正确预测的比例
F1分数 \(2 \times \frac{精确率 \times 召回率}{精确率 + 召回率}\) 精确率和召回率的调和平均

回归问题常用指标

指标 公式 含义
MAE \(\frac{1}{n}\sum|y_i - \hat{y}_i|\) 平均绝对误差
MSE \(\frac{1}{n}\sum(y_i - \hat{y}_i)^2\) 均方误差
RMSE \(\sqrt{MSE}\) 均方根误差
\(1 - \frac{SS_{res}}{SS_{tot}}\) 决定系数
# 评价指标计算示例
# 使用模拟数据演示

set.seed(42)
n <- 100

# 模拟真实值和预测值
actual <- sample(c(0, 1), n, replace = TRUE, prob = c(0.7, 0.3))
predicted_prob <- runif(n)
predicted <- ifelse(predicted_prob > 0.5, 1, 0)

# 计算混淆矩阵
conf_matrix <- table(实际 = actual, 预测 = predicted)
print(conf_matrix)
##     预测
## 实际  0  1
##    0 31 35
##    1 13 21
# 手动计算评价指标
TP <- conf_matrix[2, 2]
TN <- conf_matrix[1, 1]
FP <- conf_matrix[1, 2]
FN <- conf_matrix[2, 1]

accuracy <- (TP + TN) / sum(conf_matrix)
precision <- TP / (TP + FP)
recall <- TP / (TP + FN)
f1 <- 2 * precision * recall / (precision + recall)

cat("准确率:", accuracy, "\n")
## 准确率: 0.52
cat("精确率:", precision, "\n")
## 精确率: 0.375
cat("召回率:", recall, "\n")
## 召回率: 0.6176471
cat("F1分数:", f1, "\n")
## F1分数: 0.4666667

使用caret计算评价指标

# 使用caret包计算评价指标
set.seed(42)

# 训练模型
model <- caret::train(
  Species ~ .,
  data = iris,
  method = "rpart",
  trControl = caret::trainControl(method = "cv", number = 5)
)

# 查看详细结果
# 确保因子水平一致
obs <- factor(model$pred$obs, levels = levels(iris$Species))
pred <- factor(model$pred$pred, levels = levels(iris$Species))
conf_mat <- caret::confusionMatrix(pred, obs)
print(conf_mat)
## Confusion Matrix and Statistics
## 
##             Reference
## Prediction   setosa versicolor virginica
##   setosa          0          0         0
##   versicolor      0          0         0
##   virginica       0          0         0
## 
## Overall Statistics
##                                   
##                Accuracy : NaN     
##                  95% CI : (NA, NA)
##     No Information Rate : NA      
##     P-Value [Acc > NIR] : NA      
##                                   
##                   Kappa : NaN     
##                                   
##  Mcnemar's Test P-Value : NA      
## 
## Statistics by Class:
## 
##                      Class: setosa Class: versicolor Class: virginica
## Sensitivity                     NA                NA               NA
## Specificity                     NA                NA               NA
## Pos Pred Value                  NA                NA               NA
## Neg Pred Value                  NA                NA               NA
## Prevalence                     NaN               NaN              NaN
## Detection Rate                 NaN               NaN              NaN
## Detection Prevalence           NaN               NaN              NaN
## Balanced Accuracy               NA                NA               NA

意义与使用场景

核心意义:选择正确的评价指标是模型评估的基础,不同指标反映模型的不同侧面。

举例:癌症筛查模型的评价指标选择

场景:开发一个乳腺癌筛查模型,用于大规模人群筛查。

为什么召回率比精确率更重要?

在筛查场景中,漏诊(假阴性)的代价远高于误诊(假阳性)。漏诊可能导致患者错过最佳治疗时机,而误诊只是需要进一步检查确认。因此,应该优先优化召回率,确保尽可能多地识别出真正的患者。

如果精确率和召回率都很重要呢?

使用F1分数作为综合指标,或者根据实际成本设定不同的权重(如加权F分数)。

需要注意什么?

在类别不平衡的数据中,准确率可能具有误导性。例如,如果95%的样本是健康人,一个”全预测为健康”的模型准确率也有95%,但完全没有实用价值。


第二章 数据预处理与特征工程

2.1 特征缩放(标准化、归一化、鲁棒缩放)

为什么需要特征缩放?

不同特征的量纲和数值范围可能差异巨大,这会导致:

  1. 距离敏感算法受影响:KNN、SVM等依赖距离计算的算法会被大量纲特征主导
  2. 梯度下降收敛慢:神经网络等优化算法收敛速度受影响
  3. 正则化效果偏差:L1/L2正则化对不同量纲特征惩罚不均

三种常用缩放方法

方法 公式 结果范围 适用场景
标准化(Z-score) \(\frac{x - \mu}{\sigma}\) 均值0,标准差1 正态分布数据
归一化(Min-Max) \(\frac{x - x_{min}}{x_{max} - x_{min}}\) [0, 1] 有明确边界的数据
鲁棒缩放 \(\frac{x - Q_2}{Q_3 - Q_1}\) 无固定范围 有异常值的数据
# 特征缩放示例
set.seed(42)

# 创建模拟医疗数据
medical_data <- data.frame(
  age = round(rnorm(100, mean = 50, sd = 15)),        # 年龄:20-80岁
  blood_pressure = round(rnorm(100, mean = 120, sd = 20)),  # 血压:80-180 mmHg
  cholesterol = round(rnorm(100, mean = 200, sd = 50)),     # 胆固醇:100-400 mg/dL
  income = round(rnorm(100, mean = 50000, sd = 20000))      # 收入:1万-10万元
)

# 查看原始数据范围
summary(medical_data)
##       age        blood_pressure   cholesterol        income     
##  Min.   : 5.00   Min.   : 80.0   Min.   : 65.0   Min.   :16350  
##  1st Qu.:40.75   1st Qu.:108.0   1st Qu.:164.8   1st Qu.:39346  
##  Median :51.00   Median :118.5   Median :199.0   Median :49086  
##  Mean   :50.51   Mean   :118.2   Mean   :199.5   Mean   :50659  
##  3rd Qu.:60.00   3rd Qu.:129.2   3rd Qu.:232.8   3rd Qu.:63496  
##  Max.   :84.00   Max.   :174.0   Max.   :323.0   Max.   :98443

标准化(Z-score Normalization)

公式

\[x_{scaled} = \frac{x - \mu}{\sigma}\]

特点: 1. 转换后均值为0,标准差为1 2. 保留异常值信息 3. 适用于近似正态分布的数据

# 标准化实现
scale_standard <- function(x) {
  (x - mean(x)) / sd(x)
}

# 应用标准化
medical_data_std <- as.data.frame(lapply(medical_data, scale_standard))

# 查看标准化后的数据
summary(medical_data_std)
##       age           blood_pressure      cholesterol            income        
##  Min.   :-2.90227   Min.   :-2.11674   Min.   :-2.643720   Min.   :-1.95782  
##  1st Qu.:-0.62242   1st Qu.:-0.56723   1st Qu.:-0.682459   1st Qu.:-0.64557  
##  Median : 0.03125   Median : 0.01383   Median :-0.009044   Median :-0.08975  
##  Mean   : 0.00000   Mean   : 0.00000   Mean   : 0.000000   Mean   : 0.00000  
##  3rd Qu.: 0.60520   3rd Qu.: 0.60874   3rd Qu.: 0.654540   3rd Qu.: 0.73253  
##  Max.   : 2.13573   Max.   : 3.08519   Max.   : 2.429013   Max.   : 2.72680
# 使用caret包的preProcess函数
preproc_std <- caret::preProcess(medical_data, method = c("center", "scale"))
medical_data_std2 <- predict(preproc_std, medical_data)

# 可视化对比
viz_data <- data.frame(
  特征 = rep(names(medical_data), each = nrow(medical_data)),
  原始值 = unlist(medical_data),
  标准化后 = unlist(medical_data_std)
)

viz_data$特征 <- factor(viz_data$特征, levels = names(medical_data))

p1 <- ggplot2::ggplot(viz_data, ggplot2::aes(x = 特征, y = 原始值, fill = 特征)) +
  ggplot2::geom_boxplot() +
  ggplot2::labs(title = "原始数据分布", y = "数值") +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "none", axis.text.x = ggplot2::element_text(angle = 45, hjust = 1))

p2 <- ggplot2::ggplot(viz_data, ggplot2::aes(x = 特征, y = 标准化后, fill = 特征)) +
  ggplot2::geom_boxplot() +
  ggplot2::labs(title = "标准化后数据分布", y = "标准化数值") +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "none", axis.text.x = ggplot2::element_text(angle = 45, hjust = 1))

p1 + p2 + patchwork::plot_annotation(
  title = "特征缩放前后对比",
  subtitle = "标准化消除了量纲差异"
)

归一化(Min-Max Normalization)

公式

\[x_{scaled} = \frac{x - x_{min}}{x_{max} - x_{min}}\]

特点: 1. 转换后数据在[0, 1]范围内 2. 对异常值敏感 3. 适用于有明确边界的数据

# 归一化实现
scale_minmax <- function(x) {
  (x - min(x)) / (max(x) - min(x))
}

# 应用归一化
medical_data_norm <- as.data.frame(lapply(medical_data, scale_minmax))

# 查看归一化后的数据范围
summary(medical_data_norm)
##       age         blood_pressure    cholesterol         income      
##  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:0.4525   1st Qu.:0.2979   1st Qu.:0.3866   1st Qu.:0.2801  
##  Median :0.5823   Median :0.4096   Median :0.5194   Median :0.3988  
##  Mean   :0.5761   Mean   :0.4069   Mean   :0.5212   Mean   :0.4179  
##  3rd Qu.:0.6962   3rd Qu.:0.5239   3rd Qu.:0.6502   3rd Qu.:0.5743  
##  Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000
# 使用caret包
preproc_norm <- caret::preProcess(medical_data, method = "range")
medical_data_norm2 <- predict(preproc_norm, medical_data)

鲁棒缩放(Robust Scaling)

公式

\[x_{scaled} = \frac{x - Q_2}{IQR} = \frac{x - Q_2}{Q_3 - Q_1}\]

特点: 1. 使用中位数和四分位距,对异常值不敏感 2. 适用于有异常值的数据 3. 结果没有固定范围

# 鲁棒缩放实现
scale_robust <- function(x) {
  (x - median(x)) / IQR(x)
}

# 添加异常值演示
medical_data_outlier <- medical_data
medical_data_outlier$cholesterol[1] <- 1000  # 添加异常值

# 对比标准化和鲁棒缩放对异常值的处理
std_with_outlier <- scale_standard(medical_data_outlier$cholesterol)
robust_with_outlier <- scale_robust(medical_data_outlier$cholesterol)

cat("标准化 - 均值:", mean(std_with_outlier), "标准差:", sd(std_with_outlier), "\n")
## 标准化 - 均值: -8.632418e-17 标准差: 1
cat("鲁棒缩放 - 均值:", mean(robust_with_outlier), "标准差:", sd(robust_with_outlier), "\n")
## 鲁棒缩放 - 均值: 0.1125088 标准差: 1.331805

意义与使用场景

核心意义:特征缩放是数据预处理的基础步骤,确保不同特征对模型的贡献均衡。

举例:基于患者生理指标预测疾病风险

为什么需要特征缩放?

假设模型使用年龄(20-80岁)、血压(80-180 mmHg)、血糖(3-30 mmol/L)三个特征。如果不缩放,血压的数值范围最大,在KNN等距离敏感算法中会主导距离计算,导致模型主要根据血压做决策,忽略其他特征。

如何选择缩放方法?

  1. 数据近似正态分布:标准化
  2. 数据有明确边界(如图像像素值0-255):归一化
  3. 数据有异常值:鲁棒缩放
  4. 不确定:先尝试标准化

需要注意什么?

缩放参数(均值、标准差、最小值、最大值)必须只从训练集计算,然后应用到测试集。不能在整个数据集上计算缩放参数,否则会造成数据泄露。


2.2 编码分类变量(独热编码、标签编码、目标编码)

为什么需要编码?

大多数机器学习算法只能处理数值输入,分类变量需要转换为数值形式。

三种常用编码方法

方法 原理 适用场景 缺点
独热编码 每个类别创建一个二进制列 名义变量、类别少 维度爆炸
标签编码 每个类别映射为一个整数 有序变量 引入虚假顺序
目标编码 用目标变量的均值编码 高基数变量、树模型 容易过拟合
# 创建模拟医疗数据
set.seed(42)
patient_data <- data.frame(
  patient_id = 1:20,
  gender = sample(c("Male", "Female"), 20, replace = TRUE),
  blood_type = sample(c("A", "B", "AB", "O"), 20, replace = TRUE),
  disease_stage = sample(c("I", "II", "III", "IV"), 20, replace = TRUE),
  hospital = sample(paste0("H", 1:10), 20, replace = TRUE),  # 高基数变量
  outcome = sample(c(0, 1), 20, replace = TRUE)
)

head(patient_data)
##   patient_id gender blood_type disease_stage hospital outcome
## 1          1   Male         AB           III       H5       1
## 2          2   Male         AB           III      H10       0
## 3          3   Male          A            II       H2       1
## 4          4   Male          A            II       H6       0
## 5          5 Female         AB            IV       H6       0
## 6          6 Female          O            IV       H2       1

独热编码(One-Hot Encoding)

原理:为每个类别创建一个二进制列,存在为1,不存在为0。

# 独热编码实现
# 使用caret包的dummyVars函数

# 对性别进行独热编码
dummy_model <- caret::dummyVars(~ gender, data = patient_data)
gender_encoded <- predict(dummy_model, patient_data)

head(gender_encoded)
##   genderFemale genderMale
## 1            0          1
## 2            0          1
## 3            0          1
## 4            0          1
## 5            1          0
## 6            1          0
# 对血型进行独热编码
dummy_model2 <- caret::dummyVars(~ blood_type, data = patient_data)
blood_type_encoded <- predict(dummy_model2, patient_data)

head(blood_type_encoded)
##   blood_typeA blood_typeAB blood_typeB blood_typeO
## 1           0            1           0           0
## 2           0            1           0           0
## 3           1            0           0           0
## 4           1            0           0           0
## 5           0            1           0           0
## 6           0            0           0           1
# 使用model.matrix函数(R基础方法)
gender_onehot <- model.matrix(~ gender - 1, data = patient_data)
head(gender_onehot)
##   genderFemale genderMale
## 1            0          1
## 2            0          1
## 3            0          1
## 4            0          1
## 5            1          0
## 6            1          0

标签编码(Label Encoding)

原理:将每个类别映射为一个整数。

适用场景:有序分类变量(如疾病分期I < II < III < IV)。

# 标签编码实现
# 使用因子转换

# 疾病分期是有序变量
patient_data$disease_stage_ordered <- factor(
  patient_data$disease_stage,
  levels = c("I", "II", "III", "IV"),
  ordered = TRUE
)

# 查看编码
patient_data$disease_stage_encoded <- as.numeric(patient_data$disease_stage_ordered)

head(patient_data[, c("disease_stage", "disease_stage_ordered", "disease_stage_encoded")])
##   disease_stage disease_stage_ordered disease_stage_encoded
## 1           III                   III                     3
## 2           III                   III                     3
## 3            II                    II                     2
## 4            II                    II                     2
## 5            IV                    IV                     4
## 6            IV                    IV                     4
# 注意:对于名义变量,标签编码会引入虚假顺序
# 例如血型的标签编码
patient_data$blood_type_encoded <- as.numeric(factor(patient_data$blood_type))

# 这会引入虚假的顺序关系:A < B < AB < O,但实际上血型没有顺序

目标编码(Target Encoding)

原理:用该类别对应的目标变量均值来编码。

公式

\[x_{encoded} = \frac{\sum_{i \in C} y_i + \alpha \cdot \bar{y}}{n_C + \alpha}\]

其中 \(C\) 是当前类别,\(n_C\) 是该类别的样本数,\(\bar{y}\) 是全局目标均值,\(\alpha\) 是平滑参数。

# 目标编码实现
# 对医院变量进行目标编码

# 计算每个医院的目标变量均值
hospital_target_mean <- patient_data %>%
  dplyr::group_by(hospital) %>%
  dplyr::summarise(
    target_mean = mean(outcome),
    count = n(),
    .groups = "drop"
  )

# 全局目标均值
global_mean <- mean(patient_data$outcome)

# 平滑参数
alpha <- 5

# 计算平滑后的目标编码
hospital_target_mean <- hospital_target_mean %>%
  dplyr::mutate(
    target_encoded = (target_mean * count + global_mean * alpha) / (count + alpha)
  )

# 合并到原数据
patient_data <- patient_data %>%
  dplyr::left_join(
    hospital_target_mean %>% dplyr::select(hospital, target_encoded),
    by = "hospital"
  )

print(hospital_target_mean)
## # A tibble: 9 × 4
##   hospital target_mean count target_encoded
##   <chr>          <dbl> <int>          <dbl>
## 1 H1             1         2          0.679
## 2 H10            0.667     3          0.594
## 3 H2             0.667     3          0.594
## 4 H3             0         1          0.458
## 5 H4             1         1          0.625
## 6 H5             0.333     3          0.469
## 7 H6             0         4          0.306
## 8 H7             1         2          0.679
## 9 H8             1         1          0.625

编码方法选择指南

变量类型 推荐方法 示例 注意事项
名义变量(类别少) 独热编码 性别、血型 避免虚拟变量陷阱
名义变量(类别多) 目标编码/嵌入 医院、邮编 需要交叉验证防过拟合
有序变量 标签编码 疾病分期 确保顺序正确
二值变量 标签编码(0/1) 是否吸烟 直接映射即可

意义与使用场景

核心意义:正确的编码方式保留变量的信息,避免引入虚假关系。

举例:疾病分期预测模型

如何编码疾病分期?

疾病分期(I、II、III、IV)是有序变量,分期越高表示病情越重。应该使用标签编码,保留这种顺序关系。如果使用独热编码,模型无法学习到”IV期比I期更严重”的信息。

如何编码医院名称?

医院名称是名义变量,但可能有几十甚至上百家医院。独热编码会导致维度爆炸。可以使用目标编码,用每家医院的历史治愈率来编码,既压缩了维度,又保留了医院质量信息。

需要注意什么?

目标编码必须使用交叉验证方式计算,只用训练集的目标均值,否则会造成数据泄露(使用了测试集的标签信息)。


2.3 缺失值处理(均值/中位数填补、模型预测填补)

缺失值的类型

类型 定义 示例 处理策略
MCAR 完全随机缺失 数据录入随机遗漏 可删除或填补
MAR 随机缺失 男性更可能不报告体重 可用其他变量预测
MNAR 非随机缺失 重症患者未完成问卷 需要专门方法
# 创建带缺失值的数据
set.seed(42)
data_with_na <- data.frame(
  patient_id = 1:100,
  age = round(rnorm(100, 50, 15)),
  bmi = round(rnorm(100, 25, 4), 1),
  blood_pressure = round(rnorm(100, 120, 20)),
  cholesterol = round(rnorm(100, 200, 50))
)

# 随机引入缺失值
data_with_na$bmi[sample(1:100, 15)] <- NA
data_with_na$blood_pressure[sample(1:100, 10)] <- NA
data_with_na$cholesterol[sample(1:100, 20)] <- NA

# 查看缺失情况
missing_summary <- data.frame(
  变量 = names(data_with_na)[-1],
  缺失数 = sapply(data_with_na[-1], function(x) sum(is.na(x))),
  缺失比例 = sapply(data_with_na[-1], function(x) mean(is.na(x)) * 100)
)

print(missing_summary)
##                          变量 缺失数 缺失比例
## age                       age      0        0
## bmi                       bmi     15       15
## blood_pressure blood_pressure     10       10
## cholesterol       cholesterol     20       20
# 可视化缺失模式
visdat <- data_with_na[-1]
visdat$rown <- 1:nrow(visdat)
visdat_long <- tidyr::pivot_longer(visdat, cols = -rown, names_to = "variable", values_to = "value")
visdat_long$is_na <- is.na(visdat_long$value)

ggplot2::ggplot(visdat_long, ggplot2::aes(x = variable, y = rown, fill = is_na)) +
  ggplot2::geom_tile() +
  ggplot2::scale_fill_manual(values = c("FALSE" = "lightblue", "TRUE" = "red")) +
  ggplot2::labs(
    title = "缺失值模式可视化",
    x = "变量", y = "样本", fill = "是否缺失"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

简单填补方法

# 均值填补
data_mean_imputed <- data_with_na
data_mean_imputed$bmi[is.na(data_mean_imputed$bmi)] <- mean(data_mean_imputed$bmi, na.rm = TRUE)
data_mean_imputed$blood_pressure[is.na(data_mean_imputed$blood_pressure)] <- mean(data_mean_imputed$blood_pressure, na.rm = TRUE)
data_mean_imputed$cholesterol[is.na(data_mean_imputed$cholesterol)] <- mean(data_mean_imputed$cholesterol, na.rm = TRUE)

# 中位数填补(对异常值更稳健)
data_median_imputed <- data_with_na
data_median_imputed$bmi[is.na(data_median_imputed$bmi)] <- median(data_median_imputed$bmi, na.rm = TRUE)
data_median_imputed$blood_pressure[is.na(data_median_imputed$blood_pressure)] <- median(data_median_imputed$blood_pressure, na.rm = TRUE)
data_median_imputed$cholesterol[is.na(data_median_imputed$cholesterol)] <- median(data_median_imputed$cholesterol, na.rm = TRUE)

# 使用caret包的preProcess
preproc_impute <- caret::preProcess(data_with_na[-1], method = "medianImpute")
data_caret_imputed <- predict(preproc_impute, data_with_na[-1])

# 比较填补结果
impute_compare <- data.frame(
  变量 = c("bmi", "blood_pressure", "cholesterol"),
  原始均值 = c(mean(data_with_na$bmi, na.rm = TRUE),
                mean(data_with_na$blood_pressure, na.rm = TRUE),
                mean(data_with_na$cholesterol, na.rm = TRUE)),
  均值填补后均值 = c(mean(data_mean_imputed$bmi),
                      mean(data_mean_imputed$blood_pressure),
                      mean(data_mean_imputed$cholesterol)),
  中位数填补后均值 = c(mean(data_median_imputed$bmi),
                        mean(data_median_imputed$blood_pressure),
                        mean(data_median_imputed$cholesterol))
)

print(impute_compare)
##             变量  原始均值 均值填补后均值 中位数填补后均值
## 1            bmi  24.66353       24.66353           24.669
## 2 blood_pressure 121.02222      121.02222          121.020
## 3    cholesterol 202.73750      202.73750          201.690

K近邻填补(KNN Imputation)

原理:利用K个最相似样本的值进行填补。

# KNN填补
preproc_knn <- caret::preProcess(data_with_na[-1], method = "knnImpute")
data_knn_imputed <- predict(preproc_knn, data_with_na[-1])

# 注意:KNN填补需要先标准化数据
# caret会自动处理

# 查看填补结果
head(data_knn_imputed)
##          age        bmi blood_pressure cholesterol
## 1  1.3066901  1.4338767     -2.0204150 -0.06333831
## 2 -0.5427005  1.2663831      0.2944158  0.81587837
## 3  0.2863367 -0.8440364      1.0824434 -0.01706374
## 4  0.5414250  2.1596823      1.9689744 -0.46129954
## 5  0.3501088 -0.6597935     -1.4293944 -0.22529927
## 6 -0.1600679  0.2055902     -1.1831358 -0.13275015

模型预测填补

原理:用其他变量作为特征,训练模型预测缺失值。

# 使用mice包进行多重插补(如果可用)
# 这里演示简单的模型预测填补

# 用其他变量预测bmi的缺失值
bmi_complete <- data_with_na[!is.na(data_with_na$bmi), ]
bmi_missing <- data_with_na[is.na(data_with_na$bmi), ]

# 训练线性回归模型
bmi_model <- lm(bmi ~ age + blood_pressure + cholesterol, data = bmi_complete)

# 预测缺失值
predicted_bmi <- predict(bmi_model, newdata = bmi_missing)

# 填补
data_model_imputed <- data_with_na
data_model_imputed$bmi[is.na(data_model_imputed$bmi)] <- predicted_bmi

# 可视化填补效果
ggplot2::ggplot() +
  ggplot2::geom_histogram(ggplot2::aes(x = data_with_na$bmi, fill = "原始(含缺失)"),
                          alpha = 0.5, bins = 20, na.rm = TRUE) +
  ggplot2::geom_histogram(ggplot2::aes(x = data_model_imputed$bmi, fill = "模型填补后"),
                          alpha = 0.5, bins = 20) +
  ggplot2::labs(
    title = "模型预测填补效果",
    x = "BMI", y = "频数", fill = ""
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

缺失值处理策略选择

缺失比例 推荐方法 原因
< 5% 简单填补/删除 影响小
5-20% KNN/模型填补 保留更多信息
> 20% 考虑删除变量 信息损失严重

意义与使用场景

核心意义:合理的缺失值处理可以保留数据信息,避免引入偏差。

举例:临床试验数据分析

场景:一项药物临床试验,部分患者未完成所有随访检查。

为什么不能简单删除缺失样本?

如果删除所有有缺失值的样本,可能导致:1. 样本量大幅减少;2. 删除的样本可能有特殊性(如病情重的患者更可能失访),引入选择偏差。

如何选择填补方法?

如果缺失是随机的(MCAR),可以用均值填补。如果缺失与其他变量相关(MAR),应该用模型预测填补。如果缺失与自身值相关(MNAR),需要更复杂的方法如模式混合模型。


2.4 特征选择(过滤法、包装法、嵌入法)

为什么需要特征选择?

  1. 降低维度:减少计算成本
  2. 提高性能:去除噪声特征
  3. 增强可解释性:识别关键因素
  4. 防止过拟合:减少模型复杂度

三种特征选择方法

方法 原理 优点 缺点
过滤法 基于统计指标筛选 快速、独立于模型 忽略特征间关系
包装法 用模型性能评估特征组合 考虑特征间关系 计算成本高
嵌入法 在模型训练中自动选择 高效、模型特定 依赖特定模型
# 使用Sonar数据集演示特征选择
data(Sonar, package = "mlbench")

# 查看数据维度
cat("数据维度:", dim(Sonar), "\n")
## 数据维度: 208 61
cat("特征数量:", ncol(Sonar) - 1, "\n")
## 特征数量: 60
# 划分训练集和测试集
set.seed(42)
train_index <- caret::createDataPartition(Sonar$Class, p = 0.7, list = FALSE)
train_data <- Sonar[train_index, ]
test_data <- Sonar[-train_index, ]

过滤法(Filter Methods)

原理:基于统计指标(如相关性、方差、信息增益)评估每个特征的重要性。

# 方差阈值法:删除方差过小的特征
feature_vars <- apply(train_data[, -ncol(train_data)], 2, var)
low_var_features <- names(feature_vars[feature_vars < 0.01])

cat("低方差特征数量:", length(low_var_features), "\n")
## 低方差特征数量: 22
# 相关性过滤:删除与目标变量相关性低的特征
# 计算每个特征与目标变量的相关性
correlations <- sapply(names(train_data)[-ncol(train_data)], function(col) {
  cor(as.numeric(train_data[[col]]), as.numeric(train_data$Class))
})

# 选择相关性最高的前20个特征
top_features <- names(sort(abs(correlations), decreasing = TRUE))[1:20]

# 可视化特征相关性
cor_df <- data.frame(
  feature = names(correlations),
  correlation = abs(correlations)
)

cor_df <- cor_df[order(cor_df$correlation, decreasing = TRUE), ]
cor_df$feature <- factor(cor_df$feature, levels = cor_df$feature)

ggplot2::ggplot(head(cor_df, 30), ggplot2::aes(x = feature, y = correlation)) +
  ggplot2::geom_bar(stat = "identity", fill = "steelblue") +
  ggplot2::labs(
    title = "特征与目标变量的相关性(前30)",
    x = "特征", y = "相关系数绝对值"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1))

包装法(Wrapper Methods)

原理:通过训练模型评估特征子集的性能,搜索最优特征组合。

# 递归特征消除(Recursive Feature Elimination, RFE)
set.seed(42)

# 定义控制参数
rfe_control <- caret::rfeControl(
  functions = caret::rfFuncs,  # 使用随机森林
  method = "cv",               # 交叉验证
  number = 5                   # 5折
)

# 执行RFE
rfe_result <- caret::rfe(
  x = train_data[, -ncol(train_data)],
  y = train_data$Class,
  sizes = c(5, 10, 20, 30, 40, 50),
  rfeControl = rfe_control
)

# 查看结果
print(rfe_result)
## 
## Recursive feature selection
## 
## Outer resampling method: Cross-Validated (5 fold) 
## 
## Resampling performance over subset size:
## 
##  Variables Accuracy  Kappa AccuracySD KappaSD Selected
##          5   0.7875 0.5714    0.08620 0.17187         
##         10   0.7870 0.5691    0.07636 0.15293         
##         20   0.8208 0.6385    0.05967 0.12223         
##         30   0.8489 0.6955    0.04050 0.08250         
##         40   0.8553 0.7081    0.05310 0.10718         
##         50   0.8351 0.6672    0.03162 0.06428         
##         60   0.8694 0.7364    0.05197 0.10351        *
## 
## The top 5 variables (out of 60):
##    V11, V12, V9, V52, V36
# 最优特征数量
cat("\n最优特征数量:", rfe_result$bestSubset, "\n")
## 
## 最优特征数量: 60
cat("选中的特征:\n")
## 选中的特征:
print(rfe_result$optVariables)
##  [1] "V11" "V12" "V9"  "V52" "V36" "V48" "V10" "V13" "V47" "V28" "V46" "V21"
## [13] "V49" "V27" "V37" "V51" "V20" "V16" "V17" "V18" "V2"  "V45" "V14" "V5" 
## [25] "V1"  "V4"  "V34" "V23" "V8"  "V31" "V15" "V3"  "V33" "V19" "V44" "V59"
## [37] "V42" "V6"  "V43" "V24" "V26" "V25" "V54" "V35" "V53" "V39" "V30" "V32"
## [49] "V60" "V58" "V22" "V40" "V38" "V55" "V7"  "V56" "V29" "V41" "V57" "V50"

嵌入法(Embedded Methods)

原理:在模型训练过程中自动进行特征选择。

# 使用带正则化的模型(如LASSO)
# 使用glmnet包

# 准备数据
x_train <- as.matrix(train_data[, -ncol(train_data)])
y_train <- ifelse(train_data$Class == "M", 1, 0)

# 使用caret训练LASSO模型
set.seed(42)
lasso_model <- caret::train(
  x = x_train,
  y = train_data$Class,
  method = "glmnet",
  trControl = caret::trainControl(method = "cv", number = 5),
  tuneGrid = data.frame(alpha = 1, lambda = seq(0.001, 0.1, length.out = 20))
)

# 查看模型
print(lasso_model)
## glmnet 
## 
## 146 samples
##  60 predictor
##   2 classes: 'M', 'R' 
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 117, 118, 116, 117, 116 
## Resampling results across tuning parameters:
## 
##   lambda       Accuracy   Kappa    
##   0.001000000  0.7386864  0.4749746
##   0.006210526  0.7796223  0.5568986
##   0.011421053  0.7660755  0.5296068
##   0.016631579  0.7724959  0.5426466
##   0.021842105  0.7582102  0.5130910
##   0.027052632  0.7786700  0.5554115
##   0.032263158  0.7786700  0.5554115
##   0.037473684  0.7791461  0.5550158
##   0.042684211  0.7860427  0.5685974
##   0.047894737  0.7998522  0.5952532
##   0.053105263  0.8003448  0.5965009
##   0.058315789  0.8005747  0.5962016
##   0.063526316  0.8008046  0.5973741
##   0.068736842  0.7939080  0.5825084
##   0.073947368  0.7867652  0.5688028
##   0.079157895  0.7870115  0.5695442
##   0.084368421  0.7870115  0.5695442
##   0.089578947  0.7732020  0.5410120
##   0.094789474  0.7798686  0.5549569
##   0.100000000  0.7658292  0.5260499
## 
## Tuning parameter 'alpha' was held constant at a value of 1
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were alpha = 1 and lambda = 0.06352632.
# 获取非零系数的特征
lasso_coefs <- coef(lasso_model$finalModel, lasso_model$bestTune$lambda)
selected_features <- rownames(lasso_coefs)[lasso_coefs[, 1] != 0][-1]  # 排除截距

cat("\nLASSO选中的特征数量:", length(selected_features), "\n")
## 
## LASSO选中的特征数量: 9
cat("选中的特征:\n")
## 选中的特征:
print(selected_features)
## [1] "V11" "V12" "V21" "V28" "V36" "V45" "V48" "V49" "V52"

特征选择方法对比

方法 计算速度 考虑特征交互 推荐场景
过滤法 初步筛选
包装法(RFE) 精度要求高
嵌入法(LASSO) 高维数据

意义与使用场景

核心意义:特征选择可以提高模型性能、降低过拟合风险、增强可解释性。

举例:基因表达数据分析

场景:分析1000个基因的表达水平与疾病的关系,样本只有100例。

为什么特征选择至关重要?

特征数(1000)远大于样本数(100),如果不进行特征选择,模型会严重过拟合。需要筛选出与疾病真正相关的少数基因。

如何选择方法?

  1. 先用过滤法快速筛选(如单变量检验p值)
  2. 再用嵌入法精细选择(如LASSO)
  3. 最后用包装法验证(如RFE)

2.5 特征提取(PCA、LDA)

特征选择 vs 特征提取

方法 原理 特点
特征选择 从原始特征中选择子集 保留原始特征含义
特征提取 创建新的特征组合 新特征是原始特征的线性组合

主成分分析(PCA)

原理:找到数据方差最大的方向,将数据投影到低维空间。

数学表达

寻找投影方向 \(w\),使得投影后的方差最大:

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

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

# PCA示例:使用iris数据集
data(iris)

# 只取数值特征
iris_numeric <- iris[, 1:4]

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

# 查看主成分
print(summary(pca_result))
## Importance of components:
##                           PC1    PC2     PC3     PC4
## Standard deviation     1.7084 0.9560 0.38309 0.14393
## Proportion of Variance 0.7296 0.2285 0.03669 0.00518
## Cumulative Proportion  0.7296 0.9581 0.99482 1.00000
# 方差解释比例
var_explained <- pca_result$sdev^2 / sum(pca_result$sdev^2)

var_df <- data.frame(
  PC = paste0("PC", 1:4),
  方差解释比例 = var_explained,
  累积解释比例 = cumsum(var_explained)
)

print(var_df)
##    PC 方差解释比例 累积解释比例
## 1 PC1  0.729624454    0.7296245
## 2 PC2  0.228507618    0.9581321
## 3 PC3  0.036689219    0.9948213
## 4 PC4  0.005178709    1.0000000
# 碎石图
ggplot2::ggplot(var_df, ggplot2::aes(x = PC, y = 方差解释比例)) +
  ggplot2::geom_bar(stat = "identity", fill = "steelblue") +
  ggplot2::geom_line(ggplot2::aes(group = 1), color = "red", linewidth = 1) +
  ggplot2::geom_point(color = "red", size = 2) +
  ggplot2::labs(
    title = "PCA碎石图",
    x = "主成分", y = "方差解释比例"
  ) +
  ggplot2::theme_minimal()

# 可视化前两个主成分
pca_scores <- as.data.frame(pca_result$x)
pca_scores$Species <- iris$Species

ggplot2::ggplot(pca_scores, ggplot2::aes(x = PC1, y = PC2, color = Species)) +
  ggplot2::geom_point(size = 3, alpha = 0.7) +
  ggplot2::labs(
    title = "PCA降维可视化",
    subtitle = "前两个主成分解释了95.8%的方差",
    x = paste0("PC1 (", round(var_explained[1] * 100, 1), "%)"),
    y = paste0("PC2 (", round(var_explained[2] * 100, 1), "%)")
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

线性判别分析(LDA)

原理:找到能最好区分类别的投影方向,最大化类间方差、最小化类内方差。

数学表达

\[\max_w \frac{w^T S_B w}{w^T S_W w}\]

其中 \(S_B\) 是类间散度矩阵,\(S_W\) 是类内散度矩阵。

# LDA示例
# 使用MASS包的lda函数

lda_result <- MASS::lda(Species ~ ., data = iris)

# 查看结果
print(lda_result)
## Call:
## lda(Species ~ ., data = iris)
## 
## Prior probabilities of groups:
##     setosa versicolor  virginica 
##  0.3333333  0.3333333  0.3333333 
## 
## Group means:
##            Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa            5.006       3.428        1.462       0.246
## versicolor        5.936       2.770        4.260       1.326
## virginica         6.588       2.974        5.552       2.026
## 
## Coefficients of linear discriminants:
##                     LD1         LD2
## Sepal.Length  0.8293776 -0.02410215
## Sepal.Width   1.5344731 -2.16452123
## Petal.Length -2.2012117  0.93192121
## Petal.Width  -2.8104603 -2.83918785
## 
## Proportion of trace:
##    LD1    LD2 
## 0.9912 0.0088
# 预测并可视化
lda_pred <- predict(lda_result, iris)
lda_scores <- as.data.frame(lda_pred$x)
lda_scores$Species <- iris$Species

ggplot2::ggplot(lda_scores, ggplot2::aes(x = LD1, y = LD2, color = Species)) +
  ggplot2::geom_point(size = 3, alpha = 0.7) +
  ggplot2::labs(
    title = "LDA降维可视化",
    subtitle = "LDA寻找最佳分类方向",
    x = "LD1", y = "LD2"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

PCA vs LDA 对比

# 对比PCA和LDA的降维效果
p1 <- ggplot2::ggplot(pca_scores, ggplot2::aes(x = PC1, y = PC2, color = Species)) +
  ggplot2::geom_point(size = 2, alpha = 0.7) +
  ggplot2::labs(title = "PCA:无监督降维", x = "PC1", y = "PC2") +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "none")

p2 <- ggplot2::ggplot(lda_scores, ggplot2::aes(x = LD1, y = LD2, color = Species)) +
  ggplot2::geom_point(size = 2, alpha = 0.7) +
  ggplot2::labs(title = "LDA:有监督降维", x = "LD1", y = "LD2") +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "none")

p1 + p2 + patchwork::plot_annotation(
  title = "PCA vs LDA 降维效果对比",
  subtitle = "LDA利用类别信息,分离效果更好"
)

意义与使用场景

核心意义:特征提取可以在保留信息的同时降低维度,便于可视化和建模。

举例:基因表达数据分析

场景:分析10000个基因的表达谱,需要降维后进行可视化或建模。

选择PCA还是LDA?

如果目的是数据探索、可视化,选择PCA(无监督)。如果目的是分类,选择LDA(有监督,利用类别信息)。但LDA最多只能提取K-1个成分(K是类别数),对于二分类问题只有1个成分。


2.6 特征交互与多项式特征

为什么需要特征交互?

有时特征之间存在交互效应,单个特征的影响取决于另一个特征的值。

医学示例:药物A的疗效可能取决于患者是否同时服用药物B。

# 创建模拟数据演示特征交互
set.seed(42)
n <- 200

interaction_data <- data.frame(
  age = round(runif(n, 30, 80)),
  bmi = round(rnorm(n, 25, 4), 1),
  smoking = sample(c(0, 1), n, replace = TRUE),
  exercise = sample(c(0, 1), n, replace = TRUE)
)

# 生成目标变量(包含交互效应)
# 吸烟对风险的影响随年龄增加而增强
interaction_data$risk_score <- 
  0.5 * interaction_data$age + 
  2 * interaction_data$bmi + 
  20 * interaction_data$smoking +
  0.3 * interaction_data$age * interaction_data$smoking +  # 交互效应
  rnorm(n, 0, 5)

# 查看数据
head(interaction_data)
##   age  bmi smoking exercise risk_score
## 1  76 29.8       1        1  145.54570
## 2  77 29.2       1        1  144.57387
## 3  44 21.0       1        1   97.18772
## 4  72 32.4       1        0  143.08005
## 5  62 22.3       1        1  110.59923
## 6  56 25.4       0        1   77.80938

创建交互特征

# 手动创建交互特征
interaction_data$age_smoking <- interaction_data$age * interaction_data$smoking
interaction_data$bmi_exercise <- interaction_data$bmi * interaction_data$exercise

# 使用模型验证交互效应的重要性
model_no_interaction <- lm(risk_score ~ age + bmi + smoking + exercise, 
                           data = interaction_data)
model_with_interaction <- lm(risk_score ~ age + bmi + smoking + exercise + 
                               age:smoking + bmi:exercise, 
                             data = interaction_data)

# 比较模型
anova_result <- anova(model_no_interaction, model_with_interaction)
print(anova_result)
## Analysis of Variance Table
## 
## Model 1: risk_score ~ age + bmi + smoking + exercise
## Model 2: risk_score ~ age + bmi + smoking + exercise + age:smoking + bmi:exercise
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1    195 5809.0                                  
## 2    193 4988.8  2    820.14 15.864 4.181e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 查看交互项系数
summary(model_with_interaction)
## 
## Call:
## lm(formula = risk_score ~ age + bmi + smoking + exercise + age:smoking + 
##     bmi:exercise, data = interaction_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.3919  -3.5348  -0.0625   3.7397  17.4802 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -2.48097    3.80750  -0.652    0.515    
## age           0.54573    0.03428  15.919  < 2e-16 ***
## bmi           1.98366    0.13167  15.065  < 2e-16 ***
## smoking      22.58550    2.88465   7.830 3.17e-13 ***
## exercise      5.33206    4.77227   1.117    0.265    
## age:smoking   0.26817    0.04969   5.397 1.97e-07 ***
## bmi:exercise -0.23562    0.18994  -1.240    0.216    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.084 on 193 degrees of freedom
## Multiple R-squared:  0.9562, Adjusted R-squared:  0.9548 
## F-statistic: 702.1 on 6 and 193 DF,  p-value: < 2.2e-16

多项式特征

原理:添加特征的幂次项,捕捉非线性关系。

# 多项式特征示例
set.seed(42)
n <- 100

# 生成非线性关系数据
poly_data <- data.frame(
  x = runif(n, -5, 5)
)
poly_data$y <- 0.5 * poly_data$x^2 - 2 * poly_data$x + 3 + rnorm(n, 0, 2)

# 线性模型
model_linear <- lm(y ~ x, data = poly_data)

# 多项式模型(二次)
model_poly2 <- lm(y ~ poly(x, degree = 2, raw = TRUE), data = poly_data)

# 多项式模型(三次)
model_poly3 <- lm(y ~ poly(x, degree = 3, raw = TRUE), data = poly_data)

# 预测
poly_data$pred_linear <- predict(model_linear)
poly_data$pred_poly2 <- predict(model_poly2)
poly_data$pred_poly3 <- predict(model_poly3)

# 可视化
ggplot2::ggplot(poly_data, ggplot2::aes(x = x, y = y)) +
  ggplot2::geom_point(alpha = 0.5) +
  ggplot2::geom_line(ggplot2::aes(y = pred_linear, color = "线性"), linewidth = 1) +
  ggplot2::geom_line(ggplot2::aes(y = pred_poly2, color = "二次多项式"), linewidth = 1) +
  ggplot2::geom_line(ggplot2::aes(y = pred_poly3, color = "三次多项式"), linewidth = 1) +
  ggplot2::labs(
    title = "多项式特征捕捉非线性关系",
    x = "x", y = "y", color = "模型"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

意义与使用场景

核心意义:特征交互和多项式特征可以捕捉复杂的非线性关系。

举例:药物相互作用研究

场景:研究两种药物联合使用的效果。

为什么需要交互特征?

药物A单独使用可能效果一般,药物B单独使用也效果一般,但两者联合使用时效果显著增强(协同效应)。这种交互效应需要通过交互特征来捕捉。

需要注意什么?

添加过多交互特征和多项式特征会导致:1. 维度爆炸;2. 过拟合风险增加。应该基于领域知识选择有意义的交互特征,或使用正则化方法控制复杂度。


2.7 处理不平衡数据(过采样、欠采样、SMOTE)

类别不平衡问题

在医学数据中,类别不平衡非常常见,如罕见疾病的诊断。

问题:模型可能倾向于预测多数类,导致少数类的预测性能很差。

# 创建不平衡数据集
set.seed(42)
n <- 1000

imbalanced_data <- data.frame(
  feature1 = c(rnorm(100, mean = 2, sd = 1), rnorm(900, mean = 0, sd = 1)),
  feature2 = c(rnorm(100, mean = 2, sd = 1), rnorm(900, mean = 0, sd = 1)),
  class = factor(c(rep("Positive", 100), rep("Negative", 900)))
)

# 查看类别分布
table(imbalanced_data$class)
## 
## Negative Positive 
##      900      100
# 可视化
ggplot2::ggplot(imbalanced_data, ggplot2::aes(x = feature1, y = feature2, color = class)) +
  ggplot2::geom_point(size = 1.5, alpha = 0.6) +
  ggplot2::labs(
    title = "不平衡数据集可视化",
    subtitle = "Positive: 100 (10%), Negative: 900 (90%)",
    x = "Feature 1", y = "Feature 2"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

过采样(Oversampling)

原理:增加少数类样本,使两类数量接近。

# 过采样:复制少数类样本
positive_samples <- imbalanced_data[imbalanced_data$class == "Positive", ]
negative_samples <- imbalanced_data[imbalanced_data$class == "Negative", ]

# 简单复制过采样
oversampled <- rbind(
  positive_samples[sample(1:nrow(positive_samples), 800, replace = TRUE), ],
  negative_samples
)

# 查看平衡后的分布
table(oversampled$class)
## 
## Negative Positive 
##      900      800

欠采样(Undersampling)

原理:减少多数类样本,使两类数量接近。

# 欠采样:随机抽取多数类样本
undersampled <- rbind(
  positive_samples,
  negative_samples[sample(1:nrow(negative_samples), 100), ]
)

# 查看平衡后的分布
table(undersampled$class)
## 
## Negative Positive 
##      100      100

SMOTE(Synthetic Minority Over-sampling Technique)

原理:在少数类样本之间插值生成新样本。

# SMOTE原理示意
set.seed(42)

# 选择两个少数类样本
sample1 <- as.numeric(positive_samples[1, 1:2])
sample2 <- as.numeric(positive_samples[2, 1:2])

# 在两点之间随机插值生成新样本
alpha <- runif(1)
new_sample <- sample1 + alpha * (sample2 - sample1)

cat("原始样本1:", sample1, "\n")
## 原始样本1: 3.370958 4.325058
cat("原始样本2:", sample2, "\n")
## 原始样本2: 1.435302 2.524122
cat("生成的新样本:", new_sample, "\n")
## 生成的新样本: 1.600208 2.677551
# 可视化SMOTE
smote_viz <- data.frame(
  x = c(sample1[1], sample2[1], new_sample[1]),
  y = c(sample1[2], sample2[2], new_sample[2]),
  type = c("原始样本1", "原始样本2", "SMOTE生成")
)

ggplot2::ggplot(smote_viz, ggplot2::aes(x = x, y = y, color = type)) +
  ggplot2::geom_point(size = 3) +
  ggplot2::geom_line(data = data.frame(x = c(sample1[1], sample2[1]),
                                        y = c(sample1[2], sample2[2])),
                     ggplot2::aes(x = x, y = y), color = "gray50", linetype = "dashed") +
  ggplot2::labs(
    title = "SMOTE原理示意",
    subtitle = "在少数类样本之间插值生成新样本",
    x = "Feature 1", y = "Feature 2", color = ""
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

不同方法的对比

方法 优点 缺点 适用场景
过采样 不丢失信息 容易过拟合 数据量小
欠采样 计算效率高 丢失信息 数据量大
SMOTE 生成新样本 可能生成噪声 常用选择

意义与使用场景

核心意义:处理类别不平衡可以提高模型对少数类的识别能力。

举例:罕见疾病诊断

场景:某种罕见疾病的发病率只有1%,需要开发诊断模型。

为什么需要处理不平衡?

如果不处理,模型可能学会”总是预测健康”,准确率高达99%,但完全无法识别患者。处理不平衡后,模型会更关注少数类(患者),提高召回率。

如何选择方法?

数据量小(<1000):过采样或SMOTE。数据量大(>10000):欠采样。也可以结合使用:先欠采样多数类,再用SMOTE生成少数类。


2.8 时间特征工程(滞后、滚动统计)

时间序列特征

对于时间序列数据,需要创建能够捕捉时间依赖性的特征。

# 创建模拟时间序列数据
set.seed(42)
n_days <- 365

time_data <- data.frame(
  date = seq(as.Date("2023-01-01"), by = "day", length.out = n_days),
  patients = round(100 + 20 * sin(2 * pi * 1:n_days / 365) +  # 年周期
                   5 * sin(2 * pi * 1:n_days / 7) +            # 周周期
                   rnorm(n_days, 0, 10))                       # 随机波动
)

# 查看数据
head(time_data)
##         date patients
## 1 2023-01-01      118
## 2 2023-01-02      100
## 3 2023-01-03      107
## 4 2023-01-04      106
## 5 2023-01-05      101
## 6 2023-01-06       97
# 可视化
ggplot2::ggplot(time_data, ggplot2::aes(x = date, y = patients)) +
  ggplot2::geom_line(color = "steelblue") +
  ggplot2::labs(
    title = "医院每日门诊量时间序列",
    x = "日期", y = "门诊人数"
  ) +
  ggplot2::theme_minimal()

滞后特征(Lag Features)

原理:使用过去时间点的值作为当前时间点的特征。

# 创建滞后特征
time_data <- time_data %>%
  dplyr::mutate(
    lag_1 = dplyr::lag(patients, 1),   # 前1天
    lag_7 = dplyr::lag(patients, 7),   # 前7天(上周同一天)
    lag_30 = dplyr::lag(patients, 30)  # 前30天
  )

# 查看滞后特征
head(time_data[30:35, ], 10)
##          date patients lag_1 lag_7 lag_30
## 30 2023-01-30      108   118   111     NA
## 31 2023-01-31      117   108   122    118
## 32 2023-02-01      115   117   125    100
## 33 2023-02-02      116   115    99    107
## 34 2023-02-03      101   116   102    106
## 35 2023-02-04      116   101    92    101

滚动统计特征

原理:计算过去一段时间窗口内的统计量。

# 创建滚动统计特征
time_data <- time_data %>%
  dplyr::mutate(
    # 7天滚动均值
    rolling_mean_7 = zoo::rollmean(patients, k = 7, fill = NA, align = "right"),
    # 7天滚动标准差
    rolling_sd_7 = zoo::rollapply(patients, width = 7, FUN = sd, fill = NA, align = "right"),
    # 7天滚动最大值
    rolling_max_7 = zoo::rollapply(patients, width = 7, FUN = max, fill = NA, align = "right")
  )

# 可视化滚动统计
ggplot2::ggplot(time_data[30:n_days, ], ggplot2::aes(x = date)) +
  ggplot2::geom_line(ggplot2::aes(y = patients, color = "原始值"), alpha = 0.5) +
  ggplot2::geom_line(ggplot2::aes(y = rolling_mean_7, color = "7天滚动均值"), linewidth = 1) +
  ggplot2::labs(
    title = "滚动统计特征",
    x = "日期", y = "门诊人数", color = ""
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

时间分解特征

# 提取时间分解特征
time_data <- time_data %>%
  dplyr::mutate(
    year = as.numeric(format(date, "%Y")),
    month = as.numeric(format(date, "%m")),
    day = as.numeric(format(date, "%d")),
    weekday = as.numeric(format(date, "%w")),  # 0-6, 周日为0
    is_weekend = weekday %in% c(0, 6)
  )

# 按星期几统计
weekday_stats <- time_data %>%
  dplyr::group_by(weekday) %>%
  dplyr::summarise(
    mean_patients = mean(patients),
    sd_patients = sd(patients),
    .groups = "drop"
  )

weekday_stats$weekday_name <- c("周日", "周一", "周二", "周三", "周四", "周五", "周六")

ggplot2::ggplot(weekday_stats, ggplot2::aes(x = weekday_name, y = mean_patients)) +
  ggplot2::geom_bar(stat = "identity", fill = "steelblue") +
  ggplot2::geom_errorbar(ggplot2::aes(ymin = mean_patients - sd_patients,
                                       ymax = mean_patients + sd_patients),
                         width = 0.2) +
  ggplot2::labs(
    title = "不同星期几的门诊量分布",
    x = "", y = "平均门诊人数"
  ) +
  ggplot2::theme_minimal()

意义与使用场景

核心意义:时间特征工程可以捕捉时间序列数据中的模式和依赖性。

举例:预测医院每日门诊量

场景:医院需要预测未来一周的每日门诊量,以便合理安排医护人员。

如何构建特征?

  1. 滞后特征:前一天的门诊量、上周同一天的门诊量
  2. 滚动统计:过去7天的平均门诊量
  3. 时间分解:星期几、是否节假日、月份
  4. 趋势特征:门诊量的变化趋势

2.9 数据分箱(等频、等距、基于决策树)

什么是数据分箱?

数据分箱(Binning/Discretization):将连续变量转换为离散的区间(箱)。

用途: 1. 处理非线性关系 2. 减少噪声影响 3. 提高可解释性

# 创建模拟数据
set.seed(42)
binning_data <- data.frame(
  age = round(runif(500, 18, 90)),
  income = round(rnorm(500, mean = 50000, sd = 20000))
)

# 查看原始数据分布
p1 <- ggplot2::ggplot(binning_data, ggplot2::aes(x = age)) +
  ggplot2::geom_histogram(bins = 30, fill = "steelblue", color = "white") +
  ggplot2::labs(title = "年龄分布(原始)", x = "年龄", y = "频数") +
  ggplot2::theme_minimal()

p2 <- ggplot2::ggplot(binning_data, ggplot2::aes(x = income)) +
  ggplot2::geom_histogram(bins = 30, fill = "steelblue", color = "white") +
  ggplot2::labs(title = "收入分布(原始)", x = "收入", y = "频数") +
  ggplot2::theme_minimal()

p1 + p2

等距分箱(Equal-width Binning)

原理:将数据范围等分为若干区间。

# 等距分箱:将年龄分为5个等宽区间
binning_data$age_equal_width <- cut(
  binning_data$age,
  breaks = 5,
  labels = c("青年", "中青年", "中年", "中老年", "老年")
)

# 查看分箱结果
table(binning_data$age_equal_width)
## 
##   青年 中青年   中年 中老年   老年 
##    114     92    110     82    102
# 可视化
ggplot2::ggplot(binning_data, ggplot2::aes(x = age_equal_width)) +
  ggplot2::geom_bar(fill = "steelblue") +
  ggplot2::labs(
    title = "等距分箱结果",
    x = "年龄区间", y = "频数"
  ) +
  ggplot2::theme_minimal()

等频分箱(Equal-frequency Binning)

原理:每个箱包含大致相同数量的样本。

# 等频分箱:将年龄分为5个等频区间
binning_data$age_equal_freq <- cut(
  binning_data$age,
  breaks = quantile(binning_data$age, probs = seq(0, 1, 0.2)),
  include.lowest = TRUE,
  labels = c("Q1", "Q2", "Q3", "Q4", "Q5")
)

# 查看分箱结果
table(binning_data$age_equal_freq)
## 
##  Q1  Q2  Q3  Q4  Q5 
## 102 104  97 101  96
# 可视化
ggplot2::ggplot(binning_data, ggplot2::aes(x = age_equal_freq)) +
  ggplot2::geom_bar(fill = "steelblue") +
  ggplot2::labs(
    title = "等频分箱结果",
    x = "年龄区间(按分位数)", y = "频数"
  ) +
  ggplot2::theme_minimal()

基于决策树的分箱

原理:使用决策树找到最优分割点,使目标变量的区分度最大。

# 创建带目标变量的数据
set.seed(42)
binning_data$risk <- ifelse(
  binning_data$age < 40, 
  sample(c(0, 1), sum(binning_data$age < 40), replace = TRUE, prob = c(0.8, 0.2)),
  ifelse(
    binning_data$age < 60,
    sample(c(0, 1), sum(binning_data$age >= 40 & binning_data$age < 60), replace = TRUE, prob = c(0.5, 0.5)),
    sample(c(0, 1), sum(binning_data$age >= 60), replace = TRUE, prob = c(0.3, 0.7))
  )
)

# 使用决策树找到最优分割点
tree_model <- rpart::rpart(
  risk ~ age,
  data = binning_data,
  method = "class",
  control = rpart::rpart.control(maxdepth = 3, minbucket = 30)
)

# 查看决策树
print(tree_model)
## n= 500 
## 
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
## 
## 1) root 500 242 1 (0.4840000 0.5160000)  
##   2) age< 39.5 151  25 0 (0.8344371 0.1655629) *
##   3) age>=39.5 349 116 1 (0.3323782 0.6676218) *
# 可视化决策树
rpart.plot::rpart.plot(tree_model, type = 3, extra = 101)

# 根据决策树结果进行分箱
binning_data$age_tree_bin <- cut(
  binning_data$age,
  breaks = c(-Inf, 39.5, 59.5, Inf),
  labels = c("低风险", "中风险", "高风险")
)

# 查看分箱后的风险分布
risk_by_bin <- binning_data %>%
  dplyr::group_by(age_tree_bin) %>%
  dplyr::summarise(
    n = n(),
    risk_rate = mean(risk),
    .groups = "drop"
  )

print(risk_by_bin)
## # A tibble: 3 × 3
##   age_tree_bin     n risk_rate
##   <fct>        <int>     <dbl>
## 1 低风险         151     0.166
## 2 中风险         152     0.612
## 3 高风险         197     0.711

分箱方法对比

方法 优点 缺点 适用场景
等距分箱 简单直观 对异常值敏感 均匀分布数据
等频分箱 每箱样本均衡 边界可能不直观 偏态分布数据
决策树分箱 考虑目标变量 需要有标签 有监督学习

意义与使用场景

核心意义:数据分箱可以将复杂的非线性关系简化为分段线性关系,提高模型的可解释性。

举例:年龄与疾病风险的关系

场景:分析年龄与心血管疾病风险的关系。

为什么需要分箱?

年龄与风险的关系可能是非线性的:青年风险低,中年风险逐渐上升,老年风险高但可能趋于平稳。分箱后可以更好地捕捉这种非线性关系,同时结果更易于解释(如”60岁以上为高风险人群”)。

如何选择分箱方法?

如果目的是探索性分析,用等距分箱便于可视化。如果目的是建模,用决策树分箱可以找到最优分割点。如果数据偏态严重,用等频分箱更稳健。


第三章 模型评估与验证

3.1 分类评估指标:准确率、精确率、召回率、F1分数

混淆矩阵基础

分类模型的预测结果可以用混淆矩阵(Confusion Matrix)来表示:

预测为正 预测为负
实际为正 TP(真正例) FN(假负例)
实际为负 FP(假正例) TN(真负例)

四大核心指标

准确率(Accuracy)

\[Accuracy = \frac{TP + TN}{TP + TN + FP + FN}\]

精确率(Precision)

\[Precision = \frac{TP}{TP + FP}\]

召回率(Recall/Sensitivity)

\[Recall = \frac{TP}{TP + FN}\]

F1分数

\[F1 = 2 \times \frac{Precision \times Recall}{Precision + Recall}\]

# 创建模拟预测结果
set.seed(42)
n <- 200

# 模拟真实标签和预测概率
actual <- sample(c(0, 1), n, replace = TRUE, prob = c(0.7, 0.3))
pred_prob <- runif(n)
predicted <- ifelse(pred_prob > 0.5, 1, 0)

# 计算混淆矩阵
conf_matrix <- table(实际 = actual, 预测 = predicted)
print(conf_matrix)
##     预测
## 实际  0  1
##    0 77 52
##    1 39 32
# 提取混淆矩阵元素
TN <- conf_matrix[1, 1]
FP <- conf_matrix[1, 2]
FN <- conf_matrix[2, 1]
TP <- conf_matrix[2, 2]

# 计算各指标
accuracy <- (TP + TN) / sum(conf_matrix)
precision <- TP / (TP + FP)
recall <- TP / (TP + FN)
f1 <- 2 * precision * recall / (precision + recall)

# 结果汇总
metrics_df <- data.frame(
  指标 = c("准确率", "精确率", "召回率", "F1分数"),
  公式 = c("(TP+TN)/Total", "TP/(TP+FP)", "TP/(TP+FN)", "2*P*R/(P+R)"),
  值 = c(accuracy, precision, recall, f1)
)

print(metrics_df)
##     指标          公式        值
## 1 准确率 (TP+TN)/Total 0.5450000
## 2 精确率    TP/(TP+FP) 0.3809524
## 3 召回率    TP/(TP+FN) 0.4507042
## 4 F1分数   2*P*R/(P+R) 0.4129032

使用caret计算混淆矩阵

# 使用caret包计算完整的混淆矩阵
confusion_result <- caret::confusionMatrix(
  factor(predicted, levels = c(0, 1)),
  factor(actual, levels = c(0, 1)),
  positive = "1"
)

print(confusion_result)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  0  1
##          0 77 39
##          1 52 32
##                                           
##                Accuracy : 0.545           
##                  95% CI : (0.4733, 0.6154)
##     No Information Rate : 0.645           
##     P-Value [Acc > NIR] : 0.9986          
##                                           
##                   Kappa : 0.0457          
##                                           
##  Mcnemar's Test P-Value : 0.2084          
##                                           
##             Sensitivity : 0.4507          
##             Specificity : 0.5969          
##          Pos Pred Value : 0.3810          
##          Neg Pred Value : 0.6638          
##              Prevalence : 0.3550          
##          Detection Rate : 0.1600          
##    Detection Prevalence : 0.4200          
##       Balanced Accuracy : 0.5238          
##                                           
##        'Positive' Class : 1               
## 

指标选择指南

场景 优先指标 原因
类别平衡 准确率 整体性能
疾病筛查 召回率 不漏诊
确诊检查 精确率 不误诊
综合评估 F1分数 平衡P和R

意义与使用场景

核心意义:不同指标反映模型不同侧面的性能,需根据实际需求选择。

举例:癌症筛查模型

场景:开发一个乳腺癌筛查模型,用于大规模人群筛查。

为什么召回率更重要?

在筛查场景中,漏诊(假阴性)的代价远高于误诊(假阳性)。漏诊可能导致患者错过最佳治疗时机,而误诊只是需要进一步检查确认。因此,应该优先优化召回率,确保尽可能多地识别出真正的患者。


3.2 混淆矩阵(confusion matrix)

混淆矩阵的完整结构

对于二分类问题,混淆矩阵包含四个基本元素:

# 可视化混淆矩阵
conf_df <- data.frame(
  预测 = rep(c("预测为正", "预测为负"), each = 2),
  实际 = rep(c("实际为正", "实际为负"), 2),
  数值 = c(TP, FN, FP, TN)
)

conf_df$预测 <- factor(conf_df$预测, levels = c("预测为正", "预测为负"))
conf_df$实际 <- factor(conf_df$实际, levels = c("实际为正", "实际为负"))

ggplot2::ggplot(conf_df, ggplot2::aes(x = 预测, y = 实际, fill = 数值)) +
  ggplot2::geom_tile(color = "white", linewidth = 1) +
  ggplot2::geom_text(ggplot2::aes(label = 数值), size = 6) +
  ggplot2::scale_fill_gradient(low = "lightblue", high = "steelblue") +
  ggplot2::labs(
    title = "混淆矩阵可视化",
    x = "预测值", y = "实际值"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "none")

多分类混淆矩阵

# 多分类混淆矩阵示例
set.seed(42)
multi_actual <- sample(c("A", "B", "C"), 200, replace = TRUE)
multi_pred <- sample(c("A", "B", "C"), 200, replace = TRUE, 
                     prob = c(0.4, 0.35, 0.25))

# 计算多分类混淆矩阵
multi_conf <- table(实际 = multi_actual, 预测 = multi_pred)
print(multi_conf)
##     预测
## 实际  A  B  C
##    A 28 31 16
##    B 37 22 19
##    C 24 11 12
# 可视化
multi_conf_df <- as.data.frame(multi_conf)

ggplot2::ggplot(multi_conf_df, ggplot2::aes(x = 预测, y = 实际, fill = Freq)) +
  ggplot2::geom_tile(color = "white", linewidth = 1) +
  ggplot2::geom_text(ggplot2::aes(label = Freq), size = 5) +
  ggplot2::scale_fill_gradient(low = "white", high = "steelblue") +
  ggplot2::labs(
    title = "多分类混淆矩阵",
    x = "预测类别", y = "实际类别"
  ) +
  ggplot2::theme_minimal()

从混淆矩阵衍生更多指标

# 更多衍生指标
specificity <- TN / (TN + FP)  # 特异度
npv <- TN / (TN + FN)          # 阴性预测值
fpr <- FP / (TN + FP)          # 假正率
fnr <- FN / (TP + FN)          # 假负率

derived_metrics <- data.frame(
  指标 = c("特异度(Specificity)", "阴性预测值(NPV)", 
           "假正率(FPR)", "假负率(FNR)"),
  公式 = c("TN/(TN+FP)", "TN/(TN+FN)", "FP/(TN+FP)", "FN/(TP+FN)"),
  值 = c(specificity, npv, fpr, fnr)
)

print(derived_metrics)
##                    指标       公式        值
## 1 特异度(Specificity) TN/(TN+FP) 0.5968992
## 2     阴性预测值(NPV) TN/(TN+FN) 0.6637931
## 3         假正率(FPR) FP/(TN+FP) 0.4031008
## 4         假负率(FNR) FN/(TP+FN) 0.5492958

3.3 ROC曲线与AUC

ROC曲线原理

ROC曲线(Receiver Operating Characteristic Curve):展示不同阈值下真阳性率(TPR)和假阳性率(FPR)的关系。

  • 横轴:假阳性率 FPR = FP / (FP + TN)
  • 纵轴:真阳性率 TPR = TP / (TP + FN) = 召回率

AUC(Area Under Curve):ROC曲线下的面积,衡量模型整体分类能力。

# 使用pROC包绘制ROC曲线
set.seed(42)
n <- 200

# 模拟数据
actual <- sample(c(0, 1), n, replace = TRUE, prob = c(0.7, 0.3))
# 好的模型:预测概率与真实标签相关
pred_prob <- ifelse(actual == 1, 
                    rbeta(n, 5, 2),   # 正类:概率偏高
                    rbeta(n, 2, 5))   # 负类:概率偏低

# 计算ROC曲线
roc_obj <- pROC::roc(actual, pred_prob, quiet = TRUE)

# 绘制ROC曲线
ggplot2::ggplot(data.frame(
  specificity = roc_obj$specificities,
  sensitivity = roc_obj$sensitivities
), ggplot2::aes(x = 1 - specificity, y = sensitivity)) +
  ggplot2::geom_line(color = "steelblue", linewidth = 1) +
  ggplot2::geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "gray50") +
  ggplot2::annotate("text", x = 0.75, y = 0.25, 
                    label = paste("AUC =", round(pROC::auc(roc_obj), 3)),
                    size = 5) +
  ggplot2::labs(
    title = "ROC曲线",
    subtitle = "曲线越靠近左上角,模型性能越好",
    x = "假阳性率(1 - 特异度)", y = "真阳性率(敏感度)"
  ) +
  ggplot2::theme_minimal()

AUC的含义

AUC值 模型性能
0.5 随机猜测
0.5-0.7 较差
0.7-0.8 一般
0.8-0.9 良好
0.9-1.0 优秀
# 比较不同模型的ROC曲线
set.seed(42)

# 模拟三个不同质量的模型
# 好模型
pred_good <- ifelse(actual == 1, rbeta(n, 6, 1.5), rbeta(n, 1.5, 6))
# 中等模型
pred_medium <- ifelse(actual == 1, rbeta(n, 3, 2), rbeta(n, 2, 3))
# 差模型
pred_poor <- ifelse(actual == 1, rbeta(n, 2, 2), rbeta(n, 2, 2))

# 计算ROC
roc_good <- pROC::roc(actual, pred_good, quiet = TRUE)
roc_medium <- pROC::roc(actual, pred_medium, quiet = TRUE)
roc_poor <- pROC::roc(actual, pred_poor, quiet = TRUE)

# 创建数据框
roc_compare <- rbind(
  data.frame(FPR = 1 - roc_good$specificities, TPR = roc_good$sensitivities, 
             Model = paste("好模型 (AUC =", round(pROC::auc(roc_good), 2), ")")),
  data.frame(FPR = 1 - roc_medium$specificities, TPR = roc_medium$sensitivities,
             Model = paste("中等模型 (AUC =", round(pROC::auc(roc_medium), 2), ")")),
  data.frame(FPR = 1 - roc_poor$specificities, TPR = roc_poor$sensitivities,
             Model = paste("差模型 (AUC =", round(pROC::auc(roc_poor), 2), ")"))
)

ggplot2::ggplot(roc_compare, ggplot2::aes(x = FPR, y = TPR, color = Model)) +
  ggplot2::geom_line(linewidth = 1) +
  ggplot2::geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "gray50") +
  ggplot2::labs(
    title = "不同模型的ROC曲线对比",
    x = "假阳性率", y = "真阳性率"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

意义与使用场景

核心意义:ROC曲线和AUC提供了模型在不同阈值下的整体性能评估,不依赖于特定阈值。

举例:疾病诊断模型评估

场景:评估一个糖尿病风险预测模型。

为什么用AUC而不是准确率?

准确率依赖于阈值选择。如果阈值设为0.5,准确率可能是80%;如果阈值设为0.3,准确率可能变为75%。AUC综合了所有可能阈值下的性能,更全面地评估模型质量。


3.4 PR曲线与平均精度(AP)

PR曲线原理

PR曲线(Precision-Recall Curve):展示不同阈值下精确率和召回率的关系。

  • 横轴:召回率(Recall)
  • 纵轴:精确率(Precision)

适用场景:类别不平衡时,PR曲线比ROC曲线更能反映模型性能。

# 绘制PR曲线
set.seed(42)

# 创建不平衡数据
actual_imb <- sample(c(0, 1), 500, replace = TRUE, prob = c(0.9, 0.1))
pred_imb <- ifelse(actual_imb == 1, rbeta(500, 5, 2), rbeta(500, 2, 5))

# 计算不同阈值下的精确率和召回率
thresholds <- seq(0, 1, by = 0.01)
pr_data <- data.frame()

for (thresh in thresholds) {
  pred_class <- ifelse(pred_imb >= thresh, 1, 0)
  
  TP <- sum(pred_class == 1 & actual_imb == 1)
  FP <- sum(pred_class == 1 & actual_imb == 0)
  FN <- sum(pred_class == 0 & actual_imb == 1)
  
  precision <- ifelse(TP + FP == 0, 1, TP / (TP + FP))
  recall <- ifelse(TP + FN == 0, 0, TP / (TP + FN))
  
  pr_data <- rbind(pr_data, data.frame(
    threshold = thresh,
    precision = precision,
    recall = recall
  ))
}

# 绘制PR曲线
ggplot2::ggplot(pr_data, ggplot2::aes(x = recall, y = precision)) +
  ggplot2::geom_line(color = "steelblue", linewidth = 1) +
  ggplot2::geom_hline(yintercept = mean(actual_imb), linetype = "dashed", 
                       color = "gray50") +
  ggplot2::labs(
    title = "PR曲线",
    subtitle = "虚线表示随机分类器的基线(正类比例)",
    x = "召回率(Recall)", y = "精确率(Precision)"
  ) +
  ggplot2::theme_minimal()

ROC vs PR曲线对比

# 对比ROC和PR曲线在不平衡数据上的表现
# 计算ROC
roc_imb <- pROC::roc(actual_imb, pred_imb, quiet = TRUE)

# 创建对比图
p1 <- ggplot2::ggplot(data.frame(
  FPR = 1 - roc_imb$specificities,
  TPR = roc_imb$sensitivities
), ggplot2::aes(x = FPR, y = TPR)) +
  ggplot2::geom_line(color = "steelblue", linewidth = 1) +
  ggplot2::geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "gray50") +
  ggplot2::labs(
    title = paste("ROC曲线 (AUC =", round(pROC::auc(roc_imb), 2), ")"),
    x = "假阳性率", y = "真阳性率"
  ) +
  ggplot2::theme_minimal()

p2 <- ggplot2::ggplot(pr_data, ggplot2::aes(x = recall, y = precision)) +
  ggplot2::geom_line(color = "steelblue", linewidth = 1) +
  ggplot2::geom_hline(yintercept = mean(actual_imb), linetype = "dashed", color = "gray50") +
  ggplot2::labs(
    title = "PR曲线",
    x = "召回率", y = "精确率"
  ) +
  ggplot2::theme_minimal()

p1 + p2 + patchwork::plot_annotation(
  title = "类别不平衡数据:ROC vs PR曲线",
  subtitle = paste("正类比例:", mean(actual_imb) * 100, "%")
)

意义与使用场景

核心意义:PR曲线在类别不平衡时更能反映模型对少数类的识别能力。

举例:罕见疾病诊断

场景:诊断一种发病率为1%的罕见疾病。

为什么用PR曲线?

在极度不平衡的数据中,ROC曲线可能过于乐观。一个”全预测为健康”的模型,ROC曲线看起来还不错(因为TN很高),但PR曲线会清楚地显示模型对患者的识别能力很差。


3.5 回归评估指标

常用回归指标

指标 公式 特点
MAE \(\frac{1}{n}\sum|y_i - \hat{y}_i|\) 对异常值稳健
MSE \(\frac{1}{n}\sum(y_i - \hat{y}_i)^2\) 对大误差敏感
RMSE \(\sqrt{MSE}\) 与原数据同量纲
\(1 - \frac{SS_{res}}{SS_{tot}}\) 解释方差比例
# 回归评估指标示例
set.seed(42)
n <- 100

# 模拟数据
x <- runif(n, 0, 10)
y_true <- 2 * x + 3
y_pred <- y_true + rnorm(n, 0, 2)  # 添加预测误差

# 计算各指标
residuals <- y_true - y_pred

MAE <- mean(abs(residuals))
MSE <- mean(residuals^2)
RMSE <- sqrt(MSE)
SS_res <- sum(residuals^2)
SS_tot <- sum((y_true - mean(y_true))^2)
R2 <- 1 - SS_res / SS_tot

# 调整R²(假设有1个自变量)
n <- length(y_true)
p <- 1
R2_adj <- 1 - (1 - R2) * (n - 1) / (n - p - 1)

reg_metrics <- data.frame(
  指标 = c("MAE", "MSE", "RMSE", "R²", "调整R²"),
  值 = c(MAE, MSE, RMSE, R2, R2_adj)
)

print(reg_metrics)
##     指标        值
## 1    MAE 1.4575998
## 2    MSE 3.4293694
## 3   RMSE 1.8518557
## 4     R² 0.9050144
## 5 调整R² 0.9040451

可视化预测效果

# 预测值 vs 实际值
plot_data <- data.frame(
  实际值 = y_true,
  预测值 = y_pred
)

ggplot2::ggplot(plot_data, ggplot2::aes(x = 实际值, y = 预测值)) +
  ggplot2::geom_point(alpha = 0.6) +
  ggplot2::geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +
  ggplot2::geom_smooth(method = "lm", se = FALSE, color = "steelblue") +
  ggplot2::labs(
    title = "预测值 vs 实际值",
    subtitle = "虚线为理想预测线(y=x)",
    x = "实际值", y = "预测值"
  ) +
  ggplot2::theme_minimal()

# 残差图
resid_data <- data.frame(
  预测值 = y_pred,
  残差 = residuals
)

ggplot2::ggplot(resid_data, ggplot2::aes(x = 预测值, y = 残差)) +
  ggplot2::geom_point(alpha = 0.6) +
  ggplot2::geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  ggplot2::labs(
    title = "残差图",
    subtitle = "检查残差是否随机分布",
    x = "预测值", y = "残差"
  ) +
  ggplot2::theme_minimal()

意义与使用场景

核心意义:回归指标衡量预测值与实际值的偏差程度,R²反映模型解释力。

举例:医疗费用预测

场景:预测患者的年度医疗费用。

如何解读指标?

MAE = 2000元:平均预测误差为2000元。RMSE = 3000元:考虑到大误差的影响,误差的”典型值”为3000元。R² = 0.75:模型解释了75%的费用变异。


3.6 多分类评估

宏平均、微平均、加权平均

对于多分类问题,需要将各类别的指标进行汇总:

方法 计算方式 特点
宏平均 各类别指标的算术平均 平等对待每个类别
微平均 汇总所有TP/FP/FN后计算 受多数类影响大
加权平均 按类别样本数加权平均 考虑类别不平衡
# 多分类评估示例
set.seed(42)

# 模拟三分类数据
n <- 300
multi_actual <- sample(c("A", "B", "C"), n, replace = TRUE, prob = c(0.5, 0.3, 0.2))

# 模拟预测(有一定准确率)
multi_pred <- sapply(multi_actual, function(x) {
  if (x == "A") sample(c("A", "B", "C"), 1, prob = c(0.7, 0.2, 0.1))
  else if (x == "B") sample(c("A", "B", "C"), 1, prob = c(0.2, 0.6, 0.2))
  else sample(c("A", "B", "C"), 1, prob = c(0.1, 0.2, 0.7))
})

# 混淆矩阵
multi_conf <- table(实际 = multi_actual, 预测 = multi_pred)
print(multi_conf)
##     预测
## 实际   A   B   C
##    A 112  21  18
##    B  24  60   7
##    C   4  10  44
# 计算每个类别的指标
calc_class_metrics <- function(conf_matrix, class_name) {
  TP <- conf_matrix[class_name, class_name]
  FP <- sum(conf_matrix[, class_name]) - TP
  FN <- sum(conf_matrix[class_name, ]) - TP
  
  precision <- TP / (TP + FP)
  recall <- TP / (TP + FN)
  f1 <- 2 * precision * recall / (precision + recall)
  
  return(c(Precision = precision, Recall = recall, F1 = f1))
}

class_metrics <- data.frame()
for (cls in c("A", "B", "C")) {
  metrics <- calc_class_metrics(multi_conf, cls)
  class_metrics <- rbind(class_metrics, data.frame(
    类别 = cls,
    Precision = metrics["Precision"],
    Recall = metrics["Recall"],
    F1 = metrics["F1"]
  ))
}

# 计算宏平均
macro_avg <- data.frame(
  类别 = "宏平均",
  Precision = mean(class_metrics$Precision),
  Recall = mean(class_metrics$Recall),
  F1 = mean(class_metrics$F1)
)

# 计算加权平均
class_counts <- table(multi_actual)
weights <- class_counts / sum(class_counts)
weighted_avg <- data.frame(
  类别 = "加权平均",
  Precision = sum(class_metrics$Precision * weights),
  Recall = sum(class_metrics$Recall * weights),
  F1 = sum(class_metrics$F1 * weights)
)

# 汇总结果
all_metrics <- rbind(class_metrics, macro_avg, weighted_avg)

print(all_metrics)
##                类别 Precision    Recall        F1
## Precision         A 0.8000000 0.7417219 0.7697595
## Precision1        B 0.6593407 0.6593407 0.6593407
## Precision2        C 0.6376812 0.7586207 0.6929134
## 1            宏平均 0.6990073 0.7198944 0.7073378
## 11         加权平均 0.7259517 0.7200000 0.7214088

使用caret计算多分类指标

# 使用caret计算多分类混淆矩阵
multi_confusion <- caret::confusionMatrix(
  factor(multi_pred, levels = c("A", "B", "C")),
  factor(multi_actual, levels = c("A", "B", "C"))
)

print(multi_confusion$byClass)
##          Sensitivity Specificity Pos Pred Value Neg Pred Value Precision
## Class: A   0.7417219   0.8120805      0.8000000      0.7562500 0.8000000
## Class: B   0.6593407   0.8516746      0.6593407      0.8516746 0.6593407
## Class: C   0.7586207   0.8966942      0.6376812      0.9393939 0.6376812
##             Recall        F1 Prevalence Detection Rate Detection Prevalence
## Class: A 0.7417219 0.7697595  0.5033333      0.3733333            0.4666667
## Class: B 0.6593407 0.6593407  0.3033333      0.2000000            0.3033333
## Class: C 0.7586207 0.6929134  0.1933333      0.1466667            0.2300000
##          Balanced Accuracy
## Class: A         0.7769012
## Class: B         0.7555077
## Class: C         0.8276575

3.7 交叉验证的变体

时间序列交叉验证

对于时间序列数据,不能随机划分,必须按时间顺序。

# 时间序列交叉验证示意
n <- 100
time_series <- 1:n

# 创建时间序列CV的划分
ts_splits <- data.frame()

for (i in 1:5) {
  train_end <- 20 * i
  test_start <- train_end + 1
  test_end <- min(train_end + 20, n)
  
  ts_splits <- rbind(ts_splits, data.frame(
    fold = i,
    type = "训练集",
    start = 1,
    end = train_end
  ))
  ts_splits <- rbind(ts_splits, data.frame(
    fold = i,
    type = "测试集",
    start = test_start,
    end = test_end
  ))
}

# 可视化
ts_splits$fold <- factor(ts_splits$fold)
ts_splits$type <- factor(ts_splits$type, levels = c("训练集", "测试集"))

ggplot2::ggplot(ts_splits, ggplot2::aes(x = start, xend = end, y = fold, yend = fold, color = type)) +
  ggplot2::geom_segment(linewidth = 3) +
  ggplot2::scale_color_manual(values = c("训练集" = "steelblue", "测试集" = "salmon")) +
  ggplot2::labs(
    title = "时间序列交叉验证",
    subtitle = "训练集在前,测试集在后,保持时间顺序",
    x = "时间", y = "折数"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

分组交叉验证

当数据有分组结构时(如同一患者的多次测量),需要确保同一组的数据不会同时出现在训练集和测试集。

# 分组交叉验证示意
set.seed(42)

# 创建分组数据(5个患者,每个患者10次测量)
group_data <- data.frame(
  patient_id = rep(1:5, each = 10),
  measurement = 1:50
)

# 使用caret创建分组划分
group_folds <- caret::createFolds(group_data$patient_id, k = 5)

# 可视化分组CV
group_viz <- data.frame()
for (i in 1:5) {
  test_patients <- unique(group_data$patient_id)[group_folds[[i]]]
  test_idx <- group_data$patient_id %in% test_patients
  
  group_viz <- rbind(group_viz, data.frame(
    fold = i,
    patient = group_data$patient_id,
    type = ifelse(test_idx, "测试集", "训练集")
  ))
}

group_viz$fold <- factor(group_viz$fold)
group_viz$type <- factor(group_viz$type, levels = c("训练集", "测试集"))

ggplot2::ggplot(group_viz, ggplot2::aes(x = patient, y = fold, fill = type)) +
  ggplot2::geom_tile(color = "white", linewidth = 0.5) +
  ggplot2::scale_fill_manual(values = c("训练集" = "steelblue", "测试集" = "salmon")) +
  ggplot2::labs(
    title = "分组交叉验证",
    subtitle = "同一患者的数据要么全在训练集,要么全在测试集",
    x = "患者ID", y = "折数"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

意义与使用场景

核心意义:交叉验证变体针对特殊数据结构设计,确保评估结果可靠,避免数据泄露。

举例:预测患者再入院风险

为什么需要选择正确的交叉验证方法?

患者数据有特殊结构:同一患者可能有多次住院记录,时间序列数据有先后顺序。如果用普通K折交叉验证,同一患者的数据可能同时出现在训练集和测试集,导致过于乐观的评估结果;或者用未来的数据预测过去,不符合实际应用场景。

研究目的是什么?

  1. 时间序列CV:按时间顺序划分,用历史数据预测未来
  2. 分组CV:确保同一患者数据不跨集,评估模型对新患者的泛化能力
  3. 分层CV:保持各类别比例,适用于不平衡数据

可能得到什么样的结果?

使用正确的交叉验证方法后,模型评估结果可能比普通CV低5-10%,但这才是真实的泛化性能。例如,普通CV显示AUC=0.85,分组CV显示AUC=0.78,后者更接近实际部署效果。

需要注意什么?

选择交叉验证方法时,要考虑数据结构和实际应用场景。如果模型将用于预测新患者,必须使用分组CV;如果用于时间序列预测,必须使用时间序列CV。


3.8 模型比较的统计检验

配对t检验

比较两个模型在多次交叉验证中的性能差异。

# 模拟两个模型的交叉验证结果
set.seed(42)
k <- 10  # 10折交叉验证

# 模型A的准确率
model_a_acc <- rnorm(k, mean = 0.85, sd = 0.03)
# 模型B的准确率(略差)
model_b_acc <- rnorm(k, mean = 0.82, sd = 0.04)

# 配对t检验
t_test_result <- t.test(model_a_acc, model_b_acc, paired = TRUE)
print(t_test_result)
## 
##  Paired t-test
## 
## data:  model_a_acc and model_b_acc
## t = 2.1428, df = 9, p-value = 0.06075
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.002950963  0.108865307
## sample estimates:
## mean difference 
##      0.05295717
# 可视化比较
compare_df <- data.frame(
  fold = 1:k,
  模型A = model_a_acc,
  模型B = model_b_acc
)

compare_long <- tidyr::pivot_longer(compare_df, cols = c(模型A, 模型B),
                                     names_to = "模型", values_to = "准确率")

ggplot2::ggplot(compare_long, ggplot2::aes(x = fold, y = 准确率, color = 模型)) +
  ggplot2::geom_point(size = 3) +
  ggplot2::geom_line(ggplot2::aes(group = 模型)) +
  ggplot2::labs(
    title = "两模型交叉验证准确率比较",
    subtitle = paste("配对t检验 p值 =", round(t_test_result$p.value, 4)),
    x = "折数", y = "准确率"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

McNemar检验

比较两个分类模型在相同样本上的预测结果。

# McNemar检验示例
set.seed(42)
n <- 200

# 模拟两个模型的预测
actual <- sample(c(0, 1), n, replace = TRUE)
model1_pred <- ifelse(runif(n) < 0.8, actual, 1 - actual)  # 80%准确率
model2_pred <- ifelse(runif(n) < 0.75, actual, 1 - actual) # 75%准确率

# 创建列联表
# b01: 模型1正确,模型2错误
# b10: 模型1错误,模型2正确
b01 <- sum(model1_pred == actual & model2_pred != actual)
b10 <- sum(model1_pred != actual & model2_pred == actual)

contingency_table <- matrix(
  c(sum(model1_pred == actual & model2_pred == actual), b01,
    b10, sum(model1_pred != actual & model2_pred != actual)),
  nrow = 2,
  dimnames = list(模型2 = c("正确", "错误"), 模型1 = c("正确", "错误"))
)

print(contingency_table)
##       模型1
## 模型2 正确 错误
##   正确  121   31
##   错误   38   10
# McNemar检验
mcnemar_result <- mcnemar.test(contingency_table)
print(mcnemar_result)
## 
##  McNemar's Chi-squared test with continuity correction
## 
## data:  contingency_table
## McNemar's chi-squared = 0.52174, df = 1, p-value = 0.4701

3.9 超参数调优

贝叶斯优化

使用贝叶斯方法智能搜索最优参数。

# 贝叶斯优化示意(概念演示)
# 实际应用可使用rBayesianOptimization包

# 创建模拟的目标函数
set.seed(42)
objective_function <- function(x) {
  # 假设最优参数在x=0.7附近
  -((x - 0.7)^2) + rnorm(1, 0, 0.1)
}

# 模拟贝叶斯优化的迭代过程
iterations <- 10
results <- data.frame()

for (i in 1:iterations) {
  # 在早期随机探索,后期利用已知信息
  if (i <= 3) {
    x <- runif(1)
  } else {
    # 在已知最优点附近探索
    best_x <- results$x[which.max(results$score)]
    x <- best_x + rnorm(1, 0, 0.1)
    x <- max(0, min(1, x))  # 限制在[0,1]范围
  }
  
  score <- objective_function(x)
  results <- rbind(results, data.frame(iteration = i, x = x, score = score))
}

# 可视化优化过程
ggplot2::ggplot(results, ggplot2::aes(x = iteration, y = score)) +
  ggplot2::geom_point(size = 2, color = "steelblue") +
  ggplot2::geom_line(color = "steelblue") +
  ggplot2::labs(
    title = "贝叶斯优化过程示意",
    subtitle = "逐步逼近最优参数",
    x = "迭代次数", y = "目标函数值"
  ) +
  ggplot2::theme_minimal()

意义与使用场景

核心意义:超参数调优是提升模型性能的关键步骤,需要平衡搜索效率和计算成本。

举例:随机森林参数调优

场景:调优随机森林的树数量和最大特征数。

如何选择调优方法?

参数空间小(<100种组合):网格搜索。参数空间大:随机搜索或贝叶斯优化。计算资源有限:随机搜索。需要精确最优:贝叶斯优化。


第四章 线性模型(监督学习)

4.1 线性回归

基本原理

线性回归是机器学习中最基础的模型,假设目标变量与特征之间存在线性关系:

\[y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_p x_p + \epsilon\]

其中: - \(\beta_0\) 是截距 - \(\beta_1, \beta_2, ..., \beta_p\) 是回归系数 - \(\epsilon\) 是误差项

目标:最小化残差平方和(OLS):

\[\min_{\beta} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2\]

# 线性回归示例:医疗费用预测
set.seed(42)
n <- 200

# 模拟数据
age <- round(runif(n, 18, 80))
bmi <- rnorm(n, 25, 5)
smoker <- sample(c(0, 1), n, replace = TRUE, prob = c(0.8, 0.2))

# 医疗费用与各因素相关
medical_cost <- 1000 + 50 * age + 200 * bmi + 5000 * smoker + rnorm(n, 0, 1000)

# 创建数据框
medical_data <- data.frame(
  age = age,
  bmi = bmi,
  smoker = factor(smoker, levels = c(0, 1), labels = c("非吸烟", "吸烟")),
  cost = medical_cost
)

# 拟合线性回归模型
lm_model <- lm(cost ~ age + bmi + smoker, data = medical_data)

# 查看模型摘要
summary(lm_model)
## 
## Call:
## lm(formula = cost ~ age + bmi + smoker, data = medical_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2418.3  -631.9   -11.5   720.7  3277.6 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  980.584    441.075   2.223   0.0273 *  
## age           49.585      4.153  11.939   <2e-16 ***
## bmi          199.889     15.498  12.898   <2e-16 ***
## smoker吸烟  4922.772    191.751  25.673   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1048 on 196 degrees of freedom
## Multiple R-squared:  0.8437, Adjusted R-squared:  0.8413 
## F-statistic: 352.7 on 3 and 196 DF,  p-value: < 2.2e-16

模型解读

# 提取回归系数
coef_df <- data.frame(
  变量 = names(coef(lm_model)),
  系数 = round(coef(lm_model), 2),
  标准误 = round(summary(lm_model)$coefficients[, 2], 2),
  t值 = round(summary(lm_model)$coefficients[, 3], 2),
  p值 = round(summary(lm_model)$coefficients[, 4], 4)
)

print(coef_df)
##                    变量    系数 标准误   t值    p值
## (Intercept) (Intercept)  980.58 441.07  2.22 0.0273
## age                 age   49.59   4.15 11.94 0.0000
## bmi                 bmi  199.89  15.50 12.90 0.0000
## smoker吸烟   smoker吸烟 4922.77 191.75 25.67 0.0000

可视化回归结果

# 预测值 vs 实际值
medical_data$predicted <- predict(lm_model)

ggplot2::ggplot(medical_data, ggplot2::aes(x = cost, y = predicted, color = smoker)) +
  ggplot2::geom_point(alpha = 0.6) +
  ggplot2::geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "gray50") +
  ggplot2::labs(
    title = "线性回归:预测值 vs 实际值",
    x = "实际医疗费用", y = "预测医疗费用"
  ) +
  ggplot2::theme_minimal()

# 残差诊断
medical_data$residuals <- residuals(lm_model)

ggplot2::ggplot(medical_data, ggplot2::aes(x = predicted, y = residuals)) +
  ggplot2::geom_point(alpha = 0.6) +
  ggplot2::geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  ggplot2::geom_smooth(method = "loess", se = FALSE, color = "steelblue") +
  ggplot2::labs(
    title = "残差图",
    subtitle = "检查残差是否随机分布(无模式)",
    x = "预测值", y = "残差"
  ) +
  ggplot2::theme_minimal()

意义与使用场景

核心意义:线性回归提供了可解释性强的预测模型,系数直接反映特征对目标的影响。

举例:医疗费用预测

场景:预测患者的年度医疗费用。

如何解读系数?

age系数 = 50:年龄每增加1岁,医疗费用平均增加50元。smoker系数 = 5000:吸烟者的医疗费用平均比非吸烟者高5000元。


4.2 逻辑回归

基本原理

逻辑回归用于二分类问题,通过sigmoid函数将线性组合映射到[0,1]区间:

\[P(y=1|x) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 x_1 + ... + \beta_p x_p)}}\]

对数几率(Log-odds)

\[\log\frac{P(y=1)}{P(y=0)} = \beta_0 + \beta_1 x_1 + ... + \beta_p x_p\]

# 逻辑回归示例:疾病诊断
set.seed(42)
n <- 300

# 模拟数据
age <- round(runif(n, 30, 80))
bmi <- rnorm(n, 26, 4)
blood_pressure <- rnorm(n, 120, 15)

# 疾病风险与因素相关
log_odds <- -10 + 0.1 * age + 0.15 * bmi + 0.03 * blood_pressure
prob <- 1 / (1 + exp(-log_odds))
disease <- rbinom(n, 1, prob)

# 创建数据框
disease_data <- data.frame(
  age = age,
  bmi = bmi,
  blood_pressure = blood_pressure,
  disease = factor(disease, levels = c(0, 1), labels = c("健康", "患病"))
)

# 拟合逻辑回归模型
logit_model <- glm(disease ~ age + bmi + blood_pressure, 
                   data = disease_data, 
                   family = binomial(link = "logit"))

# 查看模型摘要
summary(logit_model)
## 
## Call:
## glm(formula = disease ~ age + bmi + blood_pressure, family = binomial(link = "logit"), 
##     data = disease_data)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -6.70130    2.59968  -2.578  0.00994 ** 
## age             0.06596    0.01692   3.897 9.73e-05 ***
## bmi             0.16000    0.06133   2.609  0.00909 ** 
## blood_pressure  0.01355    0.01424   0.952  0.34124    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 195.05  on 299  degrees of freedom
## Residual deviance: 169.26  on 296  degrees of freedom
## AIC: 177.26
## 
## Number of Fisher Scoring iterations: 6

概率输出与决策边界

# 预测概率
disease_data$prob <- predict(logit_model, type = "response")

# 默认决策边界:0.5
disease_data$pred_class <- ifelse(disease_data$prob >= 0.5, "患病", "健康")

# 可视化概率分布
ggplot2::ggplot(disease_data, ggplot2::aes(x = prob, fill = disease)) +
  ggplot2::geom_histogram(bins = 30, alpha = 0.7, position = "identity") +
  ggplot2::geom_vline(xintercept = 0.5, linetype = "dashed", color = "red") +
  ggplot2::labs(
    title = "逻辑回归:预测概率分布",
    subtitle = "虚线为默认决策边界(0.5)",
    x = "预测概率", y = "频数", fill = "实际状态"
  ) +
  ggplot2::theme_minimal()

# 决策边界可视化(以年龄和BMI为例)
ggplot2::ggplot(disease_data, ggplot2::aes(x = age, y = bmi, color = disease)) +
  ggplot2::geom_point(alpha = 0.6) +
  ggplot2::labs(
    title = "决策边界示意",
    x = "年龄", y = "BMI", color = "实际状态"
  ) +
  ggplot2::theme_minimal()

模型评估

# 计算混淆矩阵
conf_matrix <- table(实际 = disease_data$disease, 预测 = disease_data$pred_class)
print(conf_matrix)
##       预测
## 实际 患病
##   健康   30
##   患病  270
# 计算ROC和AUC
roc_obj <- pROC::roc(as.numeric(disease_data$disease) - 1, 
                     disease_data$prob, quiet = TRUE)

ggplot2::ggplot(data.frame(
  FPR = 1 - roc_obj$specificities,
  TPR = roc_obj$sensitivities
), ggplot2::aes(x = FPR, y = TPR)) +
  ggplot2::geom_line(color = "steelblue", linewidth = 1) +
  ggplot2::geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "gray50") +
  ggplot2::annotate("text", x = 0.75, y = 0.25, 
                    label = paste("AUC =", round(pROC::auc(roc_obj), 3)),
                    size = 5) +
  ggplot2::labs(
    title = "逻辑回归ROC曲线",
    x = "假阳性率", y = "真阳性率"
  ) +
  ggplot2::theme_minimal()

意义与使用场景

核心意义:逻辑回归输出概率,便于设置不同阈值平衡精确率和召回率。

举例:疾病风险评估

场景:评估患者心血管疾病风险。

为什么用逻辑回归?

逻辑回归不仅给出分类结果,还给出患病概率。可以根据临床需求调整阈值:筛查场景用低阈值(高召回率),确诊场景用高阈值(高精确率)。


4.3 正则化线性模型

为什么需要正则化?

当特征数量多或存在多重共线性时,普通线性回归容易过拟合。正则化通过添加惩罚项来约束系数。

岭回归(L2正则化)

\[\min_{\beta} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2 + \lambda \sum_{j=1}^{p} \beta_j^2\]

LASSO回归(L1正则化)

\[\min_{\beta} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2 + \lambda \sum_{j=1}^{p} |\beta_j|\]

弹性网(Elastic Net)

\[\min_{\beta} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2 + \lambda_1 \sum_{j=1}^{p} |\beta_j| + \lambda_2 \sum_{j=1}^{p} \beta_j^2\]

# 正则化示例
set.seed(42)
n <- 200
p <- 50  # 50个特征

# 模拟高维数据
X <- matrix(rnorm(n * p), n, p)
# 只有前5个特征有真实效应
true_beta <- c(3, 2, 1.5, 1, 0.5, rep(0, p - 5))
y <- X %*% true_beta + rnorm(n, 0, 1)

# 创建数据框
reg_data <- data.frame(y = y, X)
colnames(reg_data) <- c("y", paste0("X", 1:p))

# 使用caret进行正则化回归
# 岭回归(alpha = 0)
ridge_model <- caret::train(
  y ~ .,
  data = reg_data,
  method = "glmnet",
  trControl = caret::trainControl(method = "cv", number = 5),
  tuneGrid = expand.grid(alpha = 0, lambda = 10^seq(-3, 1, length.out = 20))
)

# LASSO回归(alpha = 1)
lasso_model <- caret::train(
  y ~ .,
  data = reg_data,
  method = "glmnet",
  trControl = caret::trainControl(method = "cv", number = 5),
  tuneGrid = expand.grid(alpha = 1, lambda = 10^seq(-3, 1, length.out = 20))
)

# 弹性网(alpha = 0.5)
elastic_model <- caret::train(
  y ~ .,
  data = reg_data,
  method = "glmnet",
  trControl = caret::trainControl(method = "cv", number = 5),
  tuneGrid = expand.grid(alpha = 0.5, lambda = 10^seq(-3, 1, length.out = 20))
)

# 比较模型性能
model_compare <- data.frame(
  模型 = c("岭回归", "LASSO", "弹性网"),
  最优lambda = c(ridge_model$bestTune$lambda, 
                  lasso_model$bestTune$lambda, 
                  elastic_model$bestTune$lambda),
  RMSE = c(min(ridge_model$results$RMSE), 
           min(lasso_model$results$RMSE), 
           min(elastic_model$results$RMSE))
)

print(model_compare)
##     模型 最优lambda     RMSE
## 1 岭回归  0.2069138 1.262869
## 2  LASSO  0.1274275 1.078834
## 3 弹性网  0.2069138 1.098975

LASSO的特征选择效果

# 提取LASSO系数
lasso_coef <- coef(lasso_model$finalModel, lasso_model$bestTune$lambda)
lasso_coef_df <- data.frame(
  特征 = rownames(lasso_coef),
  系数 = as.numeric(lasso_coef)
)
lasso_coef_df <- lasso_coef_df[lasso_coef_df$特征 != "(Intercept)", ]

# 可视化非零系数
nonzero_coef <- lasso_coef_df[lasso_coef_df$系数 != 0, ]
nonzero_coef <- nonzero_coef[order(abs(nonzero_coef$系数), decreasing = TRUE), ]

ggplot2::ggplot(head(nonzero_coef, 15), ggplot2::aes(x = reorder(特征, abs(系数)), y = 系数)) +
  ggplot2::geom_bar(stat = "identity", fill = "steelblue") +
  ggplot2::coord_flip() +
  ggplot2::labs(
    title = "LASSO选择的特征(非零系数)",
    x = "特征", y = "系数值"
  ) +
  ggplot2::theme_minimal()

三种正则化方法对比

方法 特点 适用场景
岭回归 系数收缩但不为0 多重共线性
LASSO 系数可收缩为0 特征选择
弹性网 结合L1和L2优点 高维+相关特征

意义与使用场景

核心意义:正则化防止过拟合,LASSO还能进行特征选择。

举例:基因表达数据分析

场景:从1000个基因中筛选与疾病相关的基因。

为什么用LASSO?

基因数量远大于样本量,普通回归会过拟合。LASSO能自动选择重要基因,将无关基因的系数压缩为0。


4.4 线性支持向量机

基本原理

支持向量机(SVM)寻找最优超平面,使两类样本之间的间隔最大化。

线性可分SVM

\[\min_{w,b} \frac{1}{2}\|w\|^2\]

约束条件:\(y_i(w^T x_i + b) \geq 1\)

# 线性SVM示例
set.seed(42)

# 创建线性可分数据
n <- 100
x1 <- c(rnorm(n/2, mean = 2, sd = 0.5), rnorm(n/2, mean = 4, sd = 0.5))
x2 <- c(rnorm(n/2, mean = 2, sd = 0.5), rnorm(n/2, mean = 4, sd = 0.5))
y <- factor(c(rep(-1, n/2), rep(1, n/2)))

svm_data <- data.frame(x1 = x1, x2 = x2, y = y)

# 训练线性SVM
svm_linear <- e1071::svm(y ~ x1 + x2, 
                          data = svm_data, 
                          type = "C-classification",
                          kernel = "linear",
                          cost = 1)

# 可视化决策边界
plot(svm_linear, svm_data, x1 ~ x2,
     main = "线性SVM决策边界",
     xlab = "特征1", ylab = "特征2")

软间隔与硬间隔

硬间隔:要求所有样本正确分类,容易过拟合。

软间隔:允许部分样本误分类,引入松弛变量\(\xi_i\)

\[\min_{w,b,\xi} \frac{1}{2}\|w\|^2 + C \sum_{i=1}^{n} \xi_i\]

C参数:控制误分类的惩罚程度。

# 比较不同C值
set.seed(42)

# 创建有噪声的数据
n <- 100
x1 <- c(rnorm(n/2, 2, 1), rnorm(n/2, 4, 1))
x2 <- c(rnorm(n/2, 2, 1), rnorm(n/2, 4, 1))
y <- factor(c(rep(-1, n/2), rep(1, n/2)))

noise_data <- data.frame(x1 = x1, x2 = x2, y = y)

# 不同C值的SVM
svm_small_c <- e1071::svm(y ~ x1 + x2, data = noise_data, 
                           kernel = "linear", cost = 0.1)
svm_large_c <- e1071::svm(y ~ x1 + x2, data = noise_data, 
                           kernel = "linear", cost = 100)

# 比较支持向量数量
sv_compare <- data.frame(
  C值 = c("0.1(软间隔)", "100(硬间隔)"),
  支持向量数 = c(length(svm_small_c$index), length(svm_large_c$index))
)

print(sv_compare)
##             C值 支持向量数
## 1 0.1(软间隔)         37
## 2 100(硬间隔)         15

4.5 感知机

基本原理

感知机(Perceptron)是最简单的神经网络,用于线性可分数据的二分类:

\[f(x) = \text{sign}(w^T x + b)\]

学习规则:对于误分类样本\((x_i, y_i)\)

\[w \leftarrow w + \eta y_i x_i\]

# 感知机示例
set.seed(42)

# 创建线性可分数据
n <- 100
perceptron_data <- data.frame(
  x1 = c(rnorm(n/2, 2, 0.5), rnorm(n/2, 4, 0.5)),
  x2 = c(rnorm(n/2, 2, 0.5), rnorm(n/2, 4, 0.5)),
  y = factor(c(rep(-1, n/2), rep(1, n/2)))
)

# 使用nnet包的单层感知机
perceptron_model <- nnet::nnet(y ~ x1 + x2, 
                                data = perceptron_data,
                                size = 0,  # 无隐藏层
                                skip = TRUE,  # 直接连接
                                maxit = 100,
                                trace = FALSE)

# 预测
perceptron_data$pred <- predict(perceptron_model, type = "class")

# 计算准确率
accuracy <- mean(perceptron_data$pred == perceptron_data$y)
print(paste("感知机准确率:", round(accuracy, 3)))
## [1] "感知机准确率: 1"

感知机的局限性

感知机只能解决线性可分问题,对于异或(XOR)等问题无法求解。

# XOR问题示例
xor_data <- data.frame(
  x1 = c(0, 0, 1, 1),
  x2 = c(0, 1, 0, 1),
  y = factor(c(0, 1, 1, 0))
)

ggplot2::ggplot(xor_data, ggplot2::aes(x = x1, y = x2, color = y)) +
  ggplot2::geom_point(size = 5) +
  ggplot2::labs(
    title = "XOR问题:线性不可分",
    subtitle = "感知机无法找到线性决策边界",
    x = "x1", y = "x2"
  ) +
  ggplot2::theme_minimal()

意义与使用场景

核心意义:感知机是神经网络的基础,理解感知机有助于理解深度学习的原理。

举例:简单医学诊断规则

为什么学习感知机?

虽然感知机在实际应用中很少单独使用,但它是理解神经网络的基础。感知机的学习规则(权重更新)是现代深度学习的核心思想。在医学领域,一些简单的线性可分问题可以用感知机快速解决,如根据体温和白细胞计数初步判断是否感染。

研究目的是什么?

  1. 理解线性分类器的工作原理
  2. 为学习多层神经网络打基础
  3. 快速处理简单线性可分问题

可能得到什么样的结果?

对于线性可分数据,感知机可以找到完美分类边界;对于线性不可分数据,感知机无法收敛,需要使用多层感知机或其他非线性方法。

需要注意什么?

感知机对初始权重敏感,可能收敛到不同的解。对于实际医学问题,建议使用更强大的方法如逻辑回归或SVM,它们能处理线性不可分数据并提供概率输出。


4.6 线性判别分析(LDA)

基本原理

线性判别分析(LDA)是一种生成式分类方法,假设各类数据服从高斯分布,寻找使类间距离最大、类内距离最小的投影方向。

Fisher判别准则

\[J(w) = \frac{w^T S_B w}{w^T S_W w}\]

其中\(S_B\)是类间散度矩阵,\(S_W\)是类内散度矩阵。

# LDA示例
set.seed(42)

# 使用iris数据集(取两个类别)
iris_binary <- iris[iris$Species %in% c("setosa", "versicolor"), ]
iris_binary$Species <- factor(iris_binary$Species)

# 训练LDA模型
lda_model <- MASS::lda(Species ~ Sepal.Length + Sepal.Width + 
                        Petal.Length + Petal.Width, 
                       data = iris_binary)

# 查看模型结果
print(lda_model)
## Call:
## lda(Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width, 
##     data = iris_binary)
## 
## Prior probabilities of groups:
##     setosa versicolor 
##        0.5        0.5 
## 
## Group means:
##            Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa            5.006       3.428        1.462       0.246
## versicolor        5.936       2.770        4.260       1.326
## 
## Coefficients of linear discriminants:
##                    LD1
## Sepal.Length -0.300458
## Sepal.Width  -1.773845
## Petal.Length  2.142260
## Petal.Width   3.035726
# 预测
lda_pred <- predict(lda_model, iris_binary)

# 混淆矩阵
table(实际 = iris_binary$Species, 预测 = lda_pred$class)
##             预测
## 实际       setosa versicolor
##   setosa         50          0
##   versicolor      0         50

LDA可视化

# LDA投影可视化
lda_values <- predict(lda_model, iris_binary)$x
iris_binary$LD1 <- lda_values[, 1]

ggplot2::ggplot(iris_binary, ggplot2::aes(x = LD1, fill = Species)) +
  ggplot2::geom_histogram(bins = 30, alpha = 0.7, position = "identity") +
  ggplot2::labs(
    title = "LDA投影结果",
    subtitle = "将多维数据投影到一维判别方向",
    x = "线性判别函数值", y = "频数"
  ) +
  ggplot2::theme_minimal()

LDA vs 逻辑回归

特点 LDA 逻辑回归
模型类型 生成式 判别式
假设 高斯分布 无分布假设
适用场景 类别分离明显 边界复杂
多分类 自然扩展 需要扩展

意义与使用场景

核心意义:LDA提供了概率框架下的分类方法,同时可用于降维。

举例:医学诊断分类

场景:根据多项指标区分良性和恶性肿瘤。

为什么用LDA?

当各类数据近似服从高斯分布且协方差矩阵相似时,LDA是最优分类器。它还提供了直观的判别函数,便于临床解释。


第五章 树模型与集成学习

5.1 决策树基础

决策树原理

决策树通过一系列规则将数据递归划分,形成树状结构:

  • 根节点:包含所有样本
  • 内部节点:特征测试
  • 叶节点:预测结果

主要算法

算法 分裂准则 特点
ID3 信息增益 只能处理分类特征
C4.5 增益率 可处理连续特征
CART 基尼系数 可用于分类和回归
# 决策树示例:疾病诊断
set.seed(42)

# 使用模拟数据
n <- 300
tree_data <- data.frame(
  age = round(runif(n, 20, 80)),
  bmi = rnorm(n, 25, 5),
  blood_pressure = rnorm(n, 120, 20),
  cholesterol = factor(sample(c("正常", "偏高", "高"), n, replace = TRUE, prob = c(0.5, 0.3, 0.2)))
)

# 生成目标变量(疾病风险)
risk_score <- 0.02 * tree_data$age + 0.05 * tree_data$bmi + 
              0.01 * tree_data$blood_pressure - 5
risk_score <- risk_score + ifelse(tree_data$cholesterol == "高", 2, 
                                   ifelse(tree_data$cholesterol == "偏高", 1, 0))
prob <- 1 / (1 + exp(-risk_score))
tree_data$disease <- factor(ifelse(runif(n) < prob, "患病", "健康"))

# 训练决策树(CART算法)
tree_model <- rpart::rpart(disease ~ age + bmi + blood_pressure + cholesterol,
                            data = tree_data,
                            method = "class",
                            control = rpart::rpart.control(maxdepth = 4))

# 查看树结构
print(tree_model)
## n= 300 
## 
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
## 
## 1) root 300 93 健康 (0.3100000 0.6900000)  
##   2) cholesterol=高,偏高 144 63 健康 (0.4375000 0.5625000)  
##     4) age>=46.5 89 41 患病 (0.5393258 0.4606742)  
##       8) bmi>=19.53831 77 32 患病 (0.5844156 0.4155844) *
##       9) bmi< 19.53831 12  3 健康 (0.2500000 0.7500000) *
##     5) age< 46.5 55 15 健康 (0.2727273 0.7272727) *
##   3) cholesterol=正常 156 30 健康 (0.1923077 0.8076923)  
##     6) bmi>=32.1859 9  3 患病 (0.6666667 0.3333333) *
##     7) bmi< 32.1859 147 24 健康 (0.1632653 0.8367347) *

决策树可视化

# 可视化决策树
rpart.plot::rpart.plot(tree_model,
                        main = "决策树:疾病诊断",
                        extra = 104,  # 显示每个节点的类别比例
                        box.palette = "GnBu")

分裂准则详解

信息增益(Information Gain)

\[IG(D, A) = H(D) - \sum_{v \in Values(A)} \frac{|D_v|}{|D|} H(D_v)\]

其中\(H(D)\)是熵:

\[H(D) = -\sum_{k=1}^{K} p_k \log_2 p_k\]

基尼系数(Gini Index)

\[Gini(D) = 1 - \sum_{k=1}^{K} p_k^2\]

# 比较不同分裂准则
# 基尼系数(默认)
tree_gini <- rpart::rpart(disease ~ age + bmi + blood_pressure + cholesterol,
                           data = tree_data, method = "class")

# 信息增益
tree_info <- rpart::rpart(disease ~ age + bmi + blood_pressure + cholesterol,
                           data = tree_data, method = "class",
                           parms = list(split = "information"))

# 比较变量重要性
importance_compare <- data.frame(
  变量 = names(tree_gini$variable.importance),
  基尼系数 = tree_gini$variable.importance,
  信息增益 = tree_info$variable.importance
)

# 标准化
importance_compare$基尼系数 <- importance_compare$基尼系数 / max(importance_compare$基尼系数)
importance_compare$信息增益 <- importance_compare$信息增益 / max(importance_compare$信息增益)

importance_long <- tidyr::pivot_longer(importance_compare, 
                                        cols = c(基尼系数, 信息增益),
                                        names_to = "准则", values_to = "重要性")

ggplot2::ggplot(importance_long, ggplot2::aes(x = 变量, y = 重要性, fill = 准则)) +
  ggplot2::geom_bar(stat = "identity", position = "dodge") +
  ggplot2::labs(
    title = "不同分裂准则的变量重要性",
    x = "变量", y = "标准化重要性"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1))


5.2 决策树的剪枝

预剪枝 vs 后剪枝

预剪枝:在构建过程中限制树的生长。

# 预剪枝参数
pre_prune <- rpart::rpart(disease ~ age + bmi + blood_pressure + cholesterol,
                           data = tree_data,
                           method = "class",
                           control = rpart::rpart.control(
                             maxdepth = 3,        # 最大深度
                             minsplit = 20,       # 分裂所需最小样本数
                             minbucket = 10       # 叶节点最小样本数
                           ))

print(paste("预剪枝树的叶节点数:", length(unique(pre_prune$where))))
## [1] "预剪枝树的叶节点数: 5"

后剪枝:先构建完整树,再根据复杂度参数剪枝。

# 后剪枝:使用cp参数
printcp(tree_model)
## 
## Classification tree:
## rpart::rpart(formula = disease ~ age + bmi + blood_pressure + 
##     cholesterol, data = tree_data, method = "class", control = rpart::rpart.control(maxdepth = 4))
## 
## Variables actually used in tree construction:
## [1] age         bmi         cholesterol
## 
## Root node error: 93/300 = 0.31
## 
## n= 300 
## 
##         CP nsplit rel error xerror     xstd
## 1 0.037634      0   1.00000 1.0000 0.086136
## 2 0.032258      3   0.86022 1.1183 0.088634
## 3 0.010000      4   0.82796 1.0860 0.088012
# 可视化cp与误差的关系
plotcp(tree_model)

# 选择最优cp
best_cp <- tree_model$cptable[which.min(tree_model$cptable[, "xerror"]), "CP"]

# 剪枝
pruned_tree <- rpart::prune(tree_model, cp = best_cp)

print(paste("剪枝前叶节点数:", length(unique(tree_model$where))))
## [1] "剪枝前叶节点数: 5"
print(paste("剪枝后叶节点数:", length(unique(pruned_tree$where))))
## [1] "剪枝后叶节点数: 1"

意义与使用场景

核心意义:剪枝防止决策树过拟合,提高泛化能力。

举例:医疗诊断决策树

场景:构建疾病诊断决策树。

为什么需要剪枝?

不剪枝的树可能过度拟合训练数据,对噪声敏感。剪枝后树更简洁,规则更通用,临床应用更可靠。


5.3 随机森林

Bagging思想

Bagging(Bootstrap Aggregating):通过自助采样训练多个模型,再聚合结果。

随机森林在Bagging基础上增加特征随机性

  • 每棵树使用bootstrap样本
  • 每次分裂只考虑随机选择的特征子集
# 随机森林示例
set.seed(42)

# 训练随机森林
rf_model <- randomForest::randomForest(
  disease ~ age + bmi + blood_pressure + cholesterol,
  data = tree_data,
  ntree = 500,      # 树的数量
  mtry = 2,         # 每次分裂考虑的特征数
  importance = TRUE
)

# 查看模型结果
print(rf_model)
## 
## Call:
##  randomForest(formula = disease ~ age + bmi + blood_pressure +      cholesterol, data = tree_data, ntree = 500, mtry = 2, importance = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 2
## 
##         OOB estimate of  error rate: 32%
## Confusion matrix:
##      患病 健康 class.error
## 患病   29   64   0.6881720
## 健康   32  175   0.1545894

OOB误差

OOB(Out-of-Bag)误差:每棵树约有1/3的样本未被用于训练,可用作验证集。

# OOB误差随树数量的变化
oob_error <- data.frame(
  ntree = 1:500,
  OOB误差 = rf_model$err.rate[, "OOB"]
)

ggplot2::ggplot(oob_error, ggplot2::aes(x = ntree, y = OOB误差)) +
  ggplot2::geom_line(color = "steelblue") +
  ggplot2::labs(
    title = "随机森林OOB误差",
    subtitle = "误差随树数量增加趋于稳定",
    x = "树的数量", y = "OOB误差率"
  ) +
  ggplot2::theme_minimal()

变量重要性

# 变量重要性
importance_df <- data.frame(
  变量 = rownames(rf_model$importance),
  平均精度下降 = rf_model$importance[, "MeanDecreaseAccuracy"],
  基尼系数下降 = rf_model$importance[, "MeanDecreaseGini"]
)

# 可视化
importance_long <- tidyr::pivot_longer(importance_df,
                                        cols = c(平均精度下降, 基尼系数下降),
                                        names_to = "指标", values_to = "值")

ggplot2::ggplot(importance_long, ggplot2::aes(x = reorder(变量, 值), y = 值, fill = 指标)) +
  ggplot2::geom_bar(stat = "identity", position = "dodge") +
  ggplot2::coord_flip() +
  ggplot2::labs(
    title = "随机森林变量重要性",
    x = "变量", y = "重要性"
  ) +
  ggplot2::theme_minimal()

意义与使用场景

核心意义:随机森林通过集成多个决策树,降低方差,提高泛化能力。

举例:疾病风险预测

场景:预测患者的心血管疾病风险。

为什么用随机森林?

随机森林能自动处理特征交互,不需要特征缩放,对缺失值和异常值鲁棒。OOB误差提供了无需额外验证集的性能评估。


5.4 梯度提升决策树(GBDT)

GBDT原理

梯度提升:通过迭代训练决策树,每棵树拟合前一棵树的残差。

\[F_m(x) = F_{m-1}(x) + \eta \cdot h_m(x)\]

其中\(h_m(x)\)是第\(m\)棵树,\(\eta\)是学习率。

# GBDT示例
set.seed(42)

# 为gbm准备数据(将因子转换为数值)
gbm_data <- tree_data
gbm_data$cholesterol_num <- as.numeric(gbm_data$cholesterol)

# 训练GBDT模型
gbm_model <- gbm::gbm(
  as.numeric(disease == "患病") ~ age + bmi + blood_pressure + cholesterol_num,
  data = gbm_data,
  distribution = "bernoulli",
  n.trees = 500,
  interaction.depth = 4,
  shrinkage = 0.01,
  cv.folds = 5,
  n.cores = 1
)
## CV: 1 
## CV: 2 
## CV: 3 
## CV: 4 
## CV: 5
# 查看最优树数量
best_iter <- gbm::gbm.perf(gbm_model, method = "cv")

print(paste("最优迭代次数:", best_iter))
## [1] "最优迭代次数: 122"

GBDT变量重要性

# 变量重要性
gbm_summary <- summary(gbm_model, n.trees = best_iter, plotit = FALSE)

ggplot2::ggplot(gbm_summary, ggplot2::aes(x = reorder(var, rel.inf), y = rel.inf)) +
  ggplot2::geom_bar(stat = "identity", fill = "steelblue") +
  ggplot2::coord_flip() +
  ggplot2::labs(
    title = "GBDT变量重要性",
    x = "变量", y = "相对影响"
  ) +
  ggplot2::theme_minimal()


5.5 XGBoost

XGBoost原理

XGBoost在GBDT基础上引入:

  1. 二阶泰勒展开:更精确的优化
  2. 正则化:控制模型复杂度
  3. 并行化:加速训练

目标函数

\[Obj = \sum_{i=1}^{n} L(y_i, \hat{y}_i) + \sum_{k=1}^{K} \Omega(f_k)\]

# XGBoost示例
set.seed(42)

# 准备数据
library(xgboost)

# 转换为数值矩阵
tree_data_num <- tree_data
tree_data_num$cholesterol <- as.numeric(factor(tree_data_num$cholesterol)) - 1
tree_data_num$disease <- as.numeric(tree_data_num$disease == "患病")

# 创建xgb矩阵
xgb_matrix <- xgboost::xgb.DMatrix(
  data = as.matrix(tree_data_num[, c("age", "bmi", "blood_pressure", "cholesterol")]),
  label = tree_data_num$disease
)

# 训练XGBoost
xgb_model <- xgboost::xgb.train(
  data = xgb_matrix,
  objective = "binary:logistic",
  nrounds = 100,
  max_depth = 4,
  eta = 0.1,
  subsample = 0.8,
  colsample_bytree = 0.8,
  verbose = 0
)

# 预测
xgb_pred <- predict(xgb_model, xgb_matrix)
xgb_pred_class <- ifelse(xgb_pred > 0.5, "患病", "健康")

# 准确率
accuracy <- mean(xgb_pred_class == tree_data$disease)
print(paste("XGBoost准确率:", round(accuracy, 3)))
## [1] "XGBoost准确率: 0.913"

XGBoost变量重要性

# 变量重要性
xgb_importance <- xgboost::xgb.importance(
  model = xgb_model,
  feature_names = c("age", "bmi", "blood_pressure", "cholesterol")
)

ggplot2::ggplot(xgb_importance, ggplot2::aes(x = reorder(Feature, Gain), y = Gain)) +
  ggplot2::geom_bar(stat = "identity", fill = "steelblue") +
  ggplot2::coord_flip() +
  ggplot2::labs(
    title = "XGBoost变量重要性",
    x = "特征", y = "增益"
  ) +
  ggplot2::theme_minimal()


5.6 LightGBM与CatBoost

LightGBM特点

LightGBM采用:

  1. 单边梯度采样(GOSS):保留大梯度样本
  2. 互斥特征捆绑(EFB):减少特征数量
  3. 叶子生长策略:按叶子分裂,更快

CatBoost特点

CatBoost特点:

  1. 类别特征处理:自动编码类别变量
  2. 对称树结构:减少过拟合
  3. 有序提升:避免目标泄露

集成方法对比

模型 优势 劣势
随机森林 简单、鲁棒 性能上限较低
GBDT 灵活、强大 需要调参
XGBoost 高效、正则化 参数较多
LightGBM 训练快、内存少 可能过拟合
CatBoost 类别特征处理 训练较慢

5.7 堆叠与混合集成

Stacking

堆叠(Stacking):用多个基模型的预测结果作为新特征,训练元模型。

# Stacking示例
set.seed(42)

# 划分训练集和测试集
train_idx <- sample(1:nrow(tree_data), 0.7 * nrow(tree_data))
train_data <- tree_data[train_idx, ]
test_data <- tree_data[-train_idx, ]

# 训练多个基模型
# 模型1:决策树
model_tree <- rpart::rpart(disease ~ ., data = train_data, method = "class")

# 模型2:随机森林
model_rf <- randomForest::randomForest(disease ~ ., data = train_data)

# 模型3:逻辑回归
model_logit <- glm(disease ~ age + bmi + blood_pressure + 
                    as.numeric(cholesterol),
                   data = train_data, family = binomial)

# 获取基模型预测概率
pred_tree <- predict(model_tree, test_data, type = "prob")[, 2]
pred_rf <- predict(model_rf, test_data, type = "prob")[, 2]
pred_logit <- predict(model_logit, test_data, type = "response")

# 创建元特征
meta_features <- data.frame(
  tree_prob = pred_tree,
  rf_prob = pred_rf,
  logit_prob = pred_logit,
  disease = test_data$disease
)

# 训练元模型(逻辑回归)
meta_model <- glm(disease ~ tree_prob + rf_prob + logit_prob,
                  data = meta_features, family = binomial)

# 最终预测
final_pred <- predict(meta_model, meta_features, type = "response")
final_pred_class <- ifelse(final_pred > 0.5, "患病", "健康")

# 评估
stacking_accuracy <- mean(final_pred_class == test_data$disease)
print(paste("Stacking准确率:", round(stacking_accuracy, 3)))
## [1] "Stacking准确率: 0.2"

Blending

混合(Blending):与Stacking类似,但使用保留验证集而非交叉验证。

投票分类器

# 投票分类器
# 硬投票
hard_vote <- ifelse(
  (predict(model_tree, test_data, type = "class") == "患病") +
  (predict(model_rf, test_data, type = "class") == "患病") +
  (ifelse(pred_logit > 0.5, "患病", "健康") == "患病") >= 2,
  "患病", "健康"
)

# 软投票
soft_vote_prob <- (pred_tree + pred_rf + pred_logit) / 3
soft_vote <- ifelse(soft_vote_prob > 0.5, "患病", "健康")

# 比较
voting_compare <- data.frame(
  方法 = c("决策树", "随机森林", "逻辑回归", "硬投票", "软投票"),
  准确率 = c(
    mean(predict(model_tree, test_data, type = "class") == test_data$disease),
    mean(predict(model_rf, test_data, type = "class") == test_data$disease),
    mean(ifelse(pred_logit > 0.5, "患病", "健康") == test_data$disease),
    mean(hard_vote == test_data$disease),
    mean(soft_vote == test_data$disease)
  )
)

print(voting_compare)
##       方法    准确率
## 1   决策树 0.7555556
## 2 随机森林 0.6888889
## 3 逻辑回归 0.2333333
## 4   硬投票 0.7111111
## 5   软投票 0.2666667

意义与使用场景

核心意义:集成学习通过组合多个模型,提高预测性能和稳定性。

举例:疾病诊断模型集成

场景:构建高准确率的疾病诊断系统。

如何选择集成方法?

模型差异大:用Stacking。模型性能相近:用投票。需要快速部署:用软投票。追求最优性能:用Stacking+调参。


第六章 支持向量机与核方法

6.1 最大间隔分类器

基本原理

最大间隔分类器寻找能够正确划分两类样本且间隔最大的超平面。

间隔:超平面到最近样本点的距离。

支持向量:距离超平面最近的样本点,决定超平面位置。

优化目标

\[\max_{w,b} \gamma \quad \text{s.t.} \quad y_i(w^T x_i + b) \geq \gamma\]

等价于:

\[\min_{w,b} \frac{1}{2}\|w\|^2 \quad \text{s.t.} \quad y_i(w^T x_i + b) \geq 1\]

# 最大间隔分类器示意
set.seed(42)

# 创建线性可分数据
n <- 50
x1 <- c(rnorm(n/2, mean = 1, sd = 0.5), rnorm(n/2, mean = 3, sd = 0.5))
x2 <- c(rnorm(n/2, mean = 1, sd = 0.5), rnorm(n/2, mean = 3, sd = 0.5))
y <- factor(c(rep(-1, n/2), rep(1, n/2)))

svm_data <- data.frame(x1 = x1, x2 = x2, y = y)

# 训练硬间隔SVM
hard_svm <- e1071::svm(y ~ x1 + x2, 
                        data = svm_data, 
                        type = "C-classification",
                        kernel = "linear",
                        cost = 1e5)  # 大cost近似硬间隔

# 可视化
plot(hard_svm, svm_data, x1 ~ x2,
     main = "最大间隔分类器",
     xlab = "特征1", ylab = "特征2",
     svSymbol = 16)  # 支持向量用方形标记

支持向量的作用

# 查看支持向量
print(paste("支持向量数量:", length(hard_svm$index)))
## [1] "支持向量数量: 3"
print(paste("总样本数:", nrow(svm_data)))
## [1] "总样本数: 50"
# 支持向量索引
print("支持向量索引:")
## [1] "支持向量索引:"
print(hard_svm$index)
## [1]  3 12 45

6.2 核函数

为什么需要核函数?

当数据线性不可分时,可以通过映射\(\phi(x)\)将数据映射到高维空间,使其线性可分。

核技巧:不需要显式计算\(\phi(x)\),直接计算核函数:

\[K(x_i, x_j) = \phi(x_i)^T \phi(x_j)\]

常用核函数

核函数 公式 特点
线性核 \(K(x, z) = x^T z\) 线性可分数据
多项式核 \(K(x, z) = (\gamma x^T z + r)^d\) 多项式边界
RBF核 \(K(x, z) = \exp(-\gamma\|x-z\|^2)\) 通用、最常用
Sigmoid核 \(K(x, z) = \tanh(\gamma x^T z + r)\) 类神经网络
# 不同核函数比较
set.seed(42)

# 创建非线性可分数据(环形数据)
n <- 200
theta <- runif(n, 0, 2 * pi)
r <- ifelse(sample(c(0, 1), n, replace = TRUE) == 0, 
            runif(n, 0, 0.5), runif(n, 1, 1.5))
x1 <- r * cos(theta)
x2 <- r * sin(theta)
y <- factor(ifelse(r < 0.7, 0, 1))

nonlinear_data <- data.frame(x1 = x1, x2 = x2, y = y)

# 可视化原始数据
ggplot2::ggplot(nonlinear_data, ggplot2::aes(x = x1, y = x2, color = y)) +
  ggplot2::geom_point(alpha = 0.7) +
  ggplot2::labs(
    title = "非线性可分数据",
    x = "特征1", y = "特征2"
  ) +
  ggplot2::theme_minimal()

线性核SVM

# 线性核
svm_linear <- e1071::svm(y ~ x1 + x2, 
                          data = nonlinear_data, 
                          kernel = "linear",
                          cost = 1)

# 预测准确率
pred_linear <- predict(svm_linear, nonlinear_data)
acc_linear <- mean(pred_linear == nonlinear_data$y)
print(paste("线性核准确率:", round(acc_linear, 3)))
## [1] "线性核准确率: 0.71"

多项式核SVM

# 多项式核
svm_poly <- e1071::svm(y ~ x1 + x2, 
                        data = nonlinear_data, 
                        kernel = "polynomial",
                        degree = 3,
                        coef0 = 1,
                        cost = 1)

pred_poly <- predict(svm_poly, nonlinear_data)
acc_poly <- mean(pred_poly == nonlinear_data$y)
print(paste("多项式核准确率:", round(acc_poly, 3)))
## [1] "多项式核准确率: 1"

RBF核SVM

# RBF核(径向基核)
svm_rbf <- e1071::svm(y ~ x1 + x2, 
                       data = nonlinear_data, 
                       kernel = "radial",
                       gamma = 1,
                       cost = 1)

pred_rbf <- predict(svm_rbf, nonlinear_data)
acc_rbf <- mean(pred_rbf == nonlinear_data$y)
print(paste("RBF核准确率:", round(acc_rbf, 3)))
## [1] "RBF核准确率: 1"

核函数对比

# 可视化不同核函数的决策边界
plot(svm_linear, nonlinear_data, x1 ~ x2, main = "线性核")

plot(svm_poly, nonlinear_data, x1 ~ x2, main = "多项式核")

plot(svm_rbf, nonlinear_data, x1 ~ x2, main = "RBF核")


6.3 核技巧

核技巧的本质

核技巧的核心思想:在高维空间中进行内积运算,而不需要显式地计算高维映射。

Mercer定理:只要\(K(x, z)\)是对称正定的,就存在映射\(\phi\)使得\(K(x, z) = \phi(x)^T \phi(z)\)

# 核技巧示意:2D到3D的映射
# 原始数据
x <- c(1, 2, 3)
z <- c(2, 3, 4)

# 多项式核 K(x,z) = (x^T z)^2
# 等价于映射 phi(x) = (x1^2, sqrt(2)*x1*x2, x2^2)

# 直接计算核函数
kernel_direct <- (sum(x * z))^2
print(paste("直接计算核函数:", kernel_direct))
## [1] "直接计算核函数: 400"
# 通过映射计算
phi_x <- c(x[1]^2, sqrt(2) * x[1] * x[2], x[2]^2)
phi_z <- c(z[1]^2, sqrt(2) * z[1] * z[2], z[2]^2)
kernel_mapped <- sum(phi_x * phi_z)
print(paste("通过映射计算:", kernel_mapped))
## [1] "通过映射计算: 64"

核函数的优势

  1. 计算效率:避免高维映射的计算
  2. 无限维映射:RBF核对应无限维空间
  3. 灵活性:不同核函数适应不同数据结构

意义与使用场景

核心意义:核技巧让线性方法能够处理非线性问题,是SVM等算法的核心技术。

举例:医学图像分类

为什么需要核技巧?

医学图像数据(如病理切片、MRI图像)在原始特征空间中往往线性不可分。例如,良性和恶性肿瘤在特征空间中可能形成复杂的非线性边界。核技巧通过隐式映射到高维空间,使原本线性不可分的数据变得线性可分,而无需显式计算高维映射。

研究目的是什么?

  1. 处理非线性分类问题
  2. 保持计算效率(避免高维映射)
  3. 灵活适应不同数据结构

可能得到什么样的结果?

使用RBF核的SVM可能在医学图像分类任务上达到90%以上的准确率,而线性SVM可能只有70%。核技巧的代价是需要调参(如RBF核的gamma参数)和更多计算资源。

需要注意什么?

核函数的选择和参数调优对结果影响很大。RBF核是最常用的选择,但对于特定问题,其他核函数可能更合适。核方法在大数据集上计算较慢,需要考虑近似方法或使用线性核。


6.4 支持向量回归(SVR)

SVR原理

支持向量回归将SVM扩展到回归问题,寻找在\(\epsilon\)-带内的超平面:

\[\min_{w,b} \frac{1}{2}\|w\|^2 + C \sum_{i=1}^{n} (\xi_i + \xi_i^*)\]

约束条件:\(|y_i - w^T x_i - b| \leq \epsilon + \xi_i\)

# SVR示例
set.seed(42)

# 创建回归数据
n <- 100
x <- runif(n, 0, 10)
y <- sin(x) + rnorm(n, 0, 0.2)

svr_data <- data.frame(x = x, y = y)

# 训练SVR模型
svr_model <- e1071::svm(y ~ x, 
                         data = svr_data, 
                         type = "eps-regression",
                         kernel = "radial",
                         epsilon = 0.1,
                         cost = 10)

# 预测
svr_pred <- predict(svr_model, svr_data)

# 可视化
ggplot2::ggplot(svr_data, ggplot2::aes(x = x, y = y)) +
  ggplot2::geom_point(alpha = 0.6) +
  ggplot2::geom_line(ggplot2::aes(y = svr_pred), color = "steelblue", linewidth = 1) +
  ggplot2::labs(
    title = "支持向量回归(SVR)",
    subtitle = "蓝线为SVR预测曲线",
    x = "x", y = "y"
  ) +
  ggplot2::theme_minimal()

epsilon参数的影响

# 比较不同epsilon值
svr_small_eps <- e1071::svm(y ~ x, data = svr_data, 
                             type = "eps-regression", kernel = "radial",
                             epsilon = 0.01, cost = 10)
svr_large_eps <- e1071::svm(y ~ x, data = svr_data, 
                             type = "eps-regression", kernel = "radial",
                             epsilon = 0.5, cost = 10)

# 比较支持向量数量
sv_compare <- data.frame(
  epsilon = c("0.01", "0.1", "0.5"),
  支持向量数 = c(length(svr_small_eps$index), 
                  length(svr_model$index), 
                  length(svr_large_eps$index))
)

print(sv_compare)
##   epsilon 支持向量数
## 1    0.01         95
## 2     0.1         67
## 3     0.5         10

6.5 软间隔与C参数

软间隔SVM

当数据存在噪声或异常值时,允许部分样本误分类:

\[\min_{w,b,\xi} \frac{1}{2}\|w\|^2 + C \sum_{i=1}^{n} \xi_i\]

C参数:控制误分类惩罚程度。

  • C大:严格惩罚误分类,可能过拟合
  • C小:容忍误分类,可能欠拟合
# C参数影响示例
set.seed(42)

# 创建有噪声的数据
n <- 100
x1 <- c(rnorm(n/2, 2, 1), rnorm(n/2, 4, 1))
x2 <- c(rnorm(n/2, 2, 1), rnorm(n/2, 4, 1))
y <- factor(c(rep(-1, n/2), rep(1, n/2)))

noisy_data <- data.frame(x1 = x1, x2 = x2, y = y)

# 不同C值
svm_c01 <- e1071::svm(y ~ x1 + x2, data = noisy_data, 
                       kernel = "radial", cost = 0.1)
svm_c1 <- e1071::svm(y ~ x1 + x2, data = noisy_data, 
                      kernel = "radial", cost = 1)
svm_c100 <- e1071::svm(y ~ x1 + x2, data = noisy_data, 
                        kernel = "radial", cost = 100)

# 比较支持向量数量
c_compare <- data.frame(
  C值 = c("0.1", "1", "100"),
  支持向量数 = c(length(svm_c01$index), length(svm_c1$index), length(svm_c100$index))
)

print(c_compare)
##   C值 支持向量数
## 1 0.1         63
## 2   1         29
## 3 100         17

6.6 核主成分分析(KPCA)

KPCA原理

核主成分分析将PCA扩展到非线性情况,通过核函数在高维空间进行主成分分析。

# KPCA示例
set.seed(42)

# 创建同心圆数据
n <- 200
theta <- runif(n, 0, 2 * pi)
r <- ifelse(sample(c(0, 1), n, replace = TRUE) == 0, 
            runif(n, 0.3, 0.5), runif(n, 0.8, 1))
x1 <- r * cos(theta)
x2 <- r * sin(theta)
y <- factor(ifelse(r < 0.6, 0, 1))

kpca_data <- data.frame(x1 = x1, x2 = x2, y = y)

# 使用kernlab包进行KPCA
if (requireNamespace("kernlab", quietly = TRUE)) {
  kpca_result <- kernlab::kpca(~ x1 + x2, 
                                data = kpca_data,
                                kernel = "rbfdot",
                                kpar = list(sigma = 5))
  
  # 提取主成分
  pc_scores <- kernlab::pcv(kpca_result)
  kpca_data$PC1 <- pc_scores[, 1]
  kpca_data$PC2 <- pc_scores[, 2]
  
  # 可视化KPCA结果
  ggplot2::ggplot(kpca_data, ggplot2::aes(x = PC1, y = PC2, color = y)) +
    ggplot2::geom_point(alpha = 0.7) +
    ggplot2::labs(
      title = "核主成分分析(KPCA)结果",
      x = "第一主成分", y = "第二主成分"
    ) +
    ggplot2::theme_minimal()
}


6.7 模型调优

核参数选择

RBF核的两个关键参数

  • C:正则化参数
  • gamma:核宽度参数
# 使用caret进行SVM调参
set.seed(42)

# 定义参数网格
tune_grid <- expand.grid(
  sigma = c(0.01, 0.1, 1, 10),
  C = c(0.1, 1, 10, 100)
)

# 使用caret调参
svm_tune <- caret::train(
  y ~ x1 + x2,
  data = nonlinear_data,
  method = "svmRadial",
  trControl = caret::trainControl(method = "cv", number = 5),
  tuneGrid = tune_grid
)

# 查看调参结果
print(svm_tune)
## Support Vector Machines with Radial Basis Function Kernel 
## 
## 200 samples
##   2 predictor
##   2 classes: '0', '1' 
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 161, 160, 159, 160, 160 
## Resampling results across tuning parameters:
## 
##   sigma  C      Accuracy   Kappa    
##    0.01    0.1  0.5350094  0.0000000
##    0.01    1.0  0.5350094  0.0000000
##    0.01   10.0  0.9900000  0.9795918
##    0.01  100.0  1.0000000  1.0000000
##    0.10    0.1  0.9008474  0.7952015
##    0.10    1.0  1.0000000  1.0000000
##    0.10   10.0  1.0000000  1.0000000
##    0.10  100.0  1.0000000  1.0000000
##    1.00    0.1  1.0000000  1.0000000
##    1.00    1.0  1.0000000  1.0000000
##    1.00   10.0  1.0000000  1.0000000
##    1.00  100.0  1.0000000  1.0000000
##   10.00    0.1  0.8900969  0.7823357
##   10.00    1.0  1.0000000  1.0000000
##   10.00   10.0  1.0000000  1.0000000
##   10.00  100.0  1.0000000  1.0000000
## 
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were sigma = 1 and C = 0.1.
# 可视化调参结果
ggplot2::ggplot(svm_tune$results, ggplot2::aes(x = sigma, y = Accuracy, color = factor(C))) +
  ggplot2::geom_point(size = 3) +
  ggplot2::geom_line() +
  ggplot2::labs(
    title = "SVM参数调优结果",
    x = "sigma(gamma)", y = "准确率", color = "C"
  ) +
  ggplot2::theme_minimal()

C与gamma的影响

参数 小值 大值
C 宽松边界,可能欠拟合 严格边界,可能过拟合
gamma 决策边界平滑 决策边界复杂

意义与使用场景

核心意义:SVM通过核技巧处理非线性问题,在中小规模数据上表现优异。

举例:医学图像分类

场景:区分良性和恶性肿瘤的医学图像。

为什么用SVM?

医学图像数据通常特征维度高、样本量中等。SVM的核方法能有效处理高维数据,且具有良好的泛化能力。RBF核能自动适应复杂的决策边界。


第七章 最近邻与基于实例的方法

7.1 K近邻(KNN)基础

基本原理

K近邻(K-Nearest Neighbors, KNN)是一种基于实例的学习方法:

  1. 计算待预测样本与所有训练样本的距离
  2. 选择距离最近的K个邻居
  3. 分类:K个邻居中多数类作为预测结果
  4. 回归:K个邻居的平均值作为预测结果
# KNN分类示例
set.seed(42)

# 创建模拟数据
n <- 200
knn_data <- data.frame(
  feature1 = c(rnorm(n/2, mean = 2, sd = 1), rnorm(n/2, mean = 5, sd = 1)),
  feature2 = c(rnorm(n/2, mean = 2, sd = 1), rnorm(n/2, mean = 5, sd = 1)),
  class = factor(c(rep("A", n/2), rep("B", n/2)))
)

# 划分训练集和测试集
train_idx <- sample(1:n, 0.7 * n)
train_data <- knn_data[train_idx, ]
test_data <- knn_data[-train_idx, ]

# 使用class包进行KNN分类
knn_pred <- class::knn(
  train = train_data[, c("feature1", "feature2")],
  test = test_data[, c("feature1", "feature2")],
  cl = train_data$class,
  k = 5
)

# 计算准确率
accuracy <- mean(knn_pred == test_data$class)
print(paste("KNN准确率:", round(accuracy, 3)))
## [1] "KNN准确率: 0.983"
# 混淆矩阵
table(实际 = test_data$class, 预测 = knn_pred)
##     预测
## 实际  A  B
##    A 30  1
##    B  0 29

可视化KNN决策边界

# 创建网格用于可视化决策边界
grid_x1 <- seq(min(knn_data$feature1) - 1, max(knn_data$feature1) + 1, by = 0.1)
grid_x2 <- seq(min(knn_data$feature2) - 1, max(knn_data$feature2) + 1, by = 0.1)
grid <- expand.grid(feature1 = grid_x1, feature2 = grid_x2)

# 预测网格点的类别
grid_pred <- class::knn(
  train = train_data[, c("feature1", "feature2")],
  test = grid,
  cl = train_data$class,
  k = 5
)

grid$class <- grid_pred

# 可视化
ggplot2::ggplot() +
  ggplot2::geom_tile(data = grid, ggplot2::aes(x = feature1, y = feature2, fill = class), alpha = 0.3) +
  ggplot2::geom_point(data = train_data, ggplot2::aes(x = feature1, y = feature2, color = class), size = 2) +
  ggplot2::labs(
    title = "KNN决策边界(K=5)",
    x = "特征1", y = "特征2"
  ) +
  ggplot2::theme_minimal()


7.2 距离度量

常用距离公式

欧氏距离(Euclidean Distance)

\[d(x, y) = \sqrt{\sum_{i=1}^{p}(x_i - y_i)^2}\]

曼哈顿距离(Manhattan Distance)

\[d(x, y) = \sum_{i=1}^{p}|x_i - y_i|\]

闵可夫斯基距离(Minkowski Distance)

\[d(x, y) = \left(\sum_{i=1}^{p}|x_i - y_i|^q\right)^{1/q}\]

\(q=2\)时为欧氏距离,\(q=1\)时为曼哈顿距离。

# 比较不同距离度量
# 欧氏距离(默认)
knn_euclidean <- class::knn(
  train = train_data[, c("feature1", "feature2")],
  test = test_data[, c("feature1", "feature2")],
  cl = train_data$class,
  k = 5
)

# 使用caret比较不同距离度量
knn_manhattan <- caret::knn3(
  class ~ feature1 + feature2,
  data = train_data,
  k = 5
)

# 预测
pred_manhattan <- predict(knn_manhattan, test_data[, c("feature1", "feature2")], type = "class")

# 比较
distance_compare <- data.frame(
  距离度量 = c("欧氏距离", "曼哈顿距离"),
  准确率 = c(
    mean(knn_euclidean == test_data$class),
    mean(pred_manhattan == test_data$class)
  )
)

print(distance_compare)
##     距离度量    准确率
## 1   欧氏距离 0.9833333
## 2 曼哈顿距离 0.9833333

7.3 加权KNN

加权原理

加权KNN:距离越近的邻居权重越大,减少K值选择的影响。

常用权重函数:

\[w_i = \frac{1}{d(x, x_i)^2}\]

# 加权KNN示例
# 使用kknn包
if (requireNamespace("kknn", quietly = TRUE)) {
  library(kknn)
  
  # 训练加权KNN
  weighted_knn <- kknn::kknn(
    class ~ feature1 + feature2,
    train = train_data,
    test = test_data,
    k = 5,
    distance = 2,  # 欧氏距离
    kernel = "gaussian"  # 高斯核权重
  )
  
  # 预测
  weighted_pred <- fitted(weighted_knn)
  
  # 比较加权与非加权
  compare_df <- data.frame(
    方法 = c("普通KNN", "加权KNN"),
    准确率 = c(
      mean(knn_euclidean == test_data$class),
      mean(weighted_pred == test_data$class)
    )
  )
  
  print(compare_df)
}

7.4 特征缩放的重要性

为什么需要特征缩放?

KNN基于距离计算,如果特征尺度差异大,大尺度特征会主导距离计算。

# 特征缩放示例
set.seed(42)

# 创建尺度差异大的数据
scale_data <- data.frame(
  age = round(runif(200, 20, 80)),           # 范围:20-80
  income = round(runif(200, 20000, 200000)), # 范围:20000-200000
  class = factor(sample(c("A", "B"), 200, replace = TRUE))
)

# 不缩放的KNN
train_idx <- sample(1:200, 150)
train_scale <- scale_data[train_idx, ]
test_scale <- scale_data[-train_idx, ]

knn_no_scale <- class::knn(
  train = train_scale[, c("age", "income")],
  test = test_scale[, c("age", "income")],
  cl = train_scale$class,
  k = 5
)

# 标准化后
train_scaled <- train_scale
test_scaled <- test_scale

# 标准化(使用训练集的均值和标准差)
age_mean <- mean(train_scaled$age)
age_sd <- sd(train_scaled$age)
income_mean <- mean(train_scaled$income)
income_sd <- sd(train_scaled$income)

train_scaled$age <- (train_scaled$age - age_mean) / age_sd
train_scaled$income <- (train_scaled$income - income_mean) / income_sd
test_scaled$age <- (test_scaled$age - age_mean) / age_sd
test_scaled$income <- (test_scaled$income - income_mean) / income_sd

knn_scaled <- class::knn(
  train = train_scaled[, c("age", "income")],
  test = test_scaled[, c("age", "income")],
  cl = train_scaled$class,
  k = 5
)

# 比较
scale_compare <- data.frame(
  方法 = c("未缩放", "标准化后"),
  准确率 = c(
    mean(knn_no_scale == test_scale$class),
    mean(knn_scaled == test_scale$class)
  )
)

print(scale_compare)
##       方法 准确率
## 1   未缩放   0.48
## 2 标准化后   0.50

7.5 K值选择

K值的影响

  • K小:模型复杂,容易过拟合,对噪声敏感
  • K大:模型简单,可能欠拟合,忽略局部结构
# 不同K值比较
set.seed(42)

k_values <- c(1, 3, 5, 7, 10, 15, 20, 30)
accuracies <- numeric(length(k_values))

for (i in seq_along(k_values)) {
  knn_pred <- class::knn(
    train = train_data[, c("feature1", "feature2")],
    test = test_data[, c("feature1", "feature2")],
    cl = train_data$class,
    k = k_values[i]
  )
  accuracies[i] <- mean(knn_pred == test_data$class)
}

k_compare <- data.frame(
  K值 = k_values,
  准确率 = round(accuracies, 3)
)

print(k_compare)
##   K值 准确率
## 1   1  0.983
## 2   3  0.983
## 3   5  0.983
## 4   7  0.983
## 5  10  0.983
## 6  15  0.983
## 7  20  0.967
## 8  30  0.983
# 可视化
ggplot2::ggplot(k_compare, ggplot2::aes(x = K值, y = 准确率)) +
  ggplot2::geom_line(color = "steelblue", linewidth = 1) +
  ggplot2::geom_point(size = 3, color = "steelblue") +
  ggplot2::labs(
    title = "K值对KNN准确率的影响",
    x = "K值", y = "准确率"
  ) +
  ggplot2::theme_minimal()

使用交叉验证选择K

# 使用caret进行交叉验证选择K
set.seed(42)

tune_control <- caret::trainControl(method = "cv", number = 5)

knn_tune <- caret::train(
  class ~ feature1 + feature2,
  data = train_data,
  method = "knn",
  trControl = tune_control,
  tuneGrid = data.frame(k = 1:20)
)

print(knn_tune)
## k-Nearest Neighbors 
## 
## 140 samples
##   2 predictor
##   2 classes: 'A', 'B' 
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 112, 112, 112, 113, 111 
## Resampling results across tuning parameters:
## 
##   k   Accuracy   Kappa    
##    1  0.9785532  0.9569958
##    2  0.9714103  0.9427758
##    3  0.9925926  0.9851240
##    4  0.9925926  0.9851240
##    5  0.9925926  0.9851240
##    6  0.9925926  0.9851240
##    7  0.9925926  0.9851240
##    8  0.9856960  0.9712815
##    9  0.9925926  0.9851240
##   10  0.9856960  0.9712815
##   11  0.9856960  0.9712815
##   12  0.9856960  0.9712815
##   13  0.9856960  0.9712815
##   14  0.9925926  0.9851240
##   15  0.9856960  0.9712815
##   16  0.9925926  0.9851240
##   17  0.9925926  0.9851240
##   18  0.9925926  0.9851240
##   19  0.9925926  0.9851240
##   20  0.9925926  0.9851240
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was k = 20.
# 可视化
ggplot2::ggplot(knn_tune$results, ggplot2::aes(x = k, y = Accuracy)) +
  ggplot2::geom_line(color = "steelblue", linewidth = 1) +
  ggplot2::geom_point(size = 2, color = "steelblue") +
  ggplot2::geom_errorbar(ggplot2::aes(ymin = Accuracy - AccuracySD,
                                       ymax = Accuracy + AccuracySD),
                         width = 0.5) +
  ggplot2::labs(
    title = "交叉验证选择最优K值",
    x = "K值", y = "准确率"
  ) +
  ggplot2::theme_minimal()


7.6 K-D树与球树

加速最近邻搜索

对于大规模数据,暴力搜索效率低。K-D树球树是常用的加速结构。

K-D树:递归地将空间划分为超矩形区域。

球树:使用超球体划分空间,对高维数据更有效。

# K-D树加速示意(概念演示)
# R中可以使用RANN包进行快速最近邻搜索

if (requireNamespace("RANN", quietly = TRUE)) {
  library(RANN)
  
  # 使用nn2函数快速查找最近邻
  nn_result <- RANN::nn2(
    data = train_data[, c("feature1", "feature2")],
    query = test_data[, c("feature1", "feature2")],
    k = 5
  )
  
  # nn_result$nn.idx 包含最近邻索引
  # nn_result$nn.dists 包含距离
  
  print("最近邻索引(前5个测试样本):")
  print(head(nn_result$nn.idx))
}
## [1] "最近邻索引(前5个测试样本):"
##      [,1] [,2] [,3] [,4] [,5]
## [1,]   45   18   57   76   44
## [2,]   23  127   16   52  128
## [3,]   24   34  104   96   71
## [4,]   23  127   16   52   20
## [5,]   57   77   18   45   76
## [6,]   95  106   49   35  100

7.7 最近邻回归

KNN回归

KNN也可用于回归问题,预测值为K个邻居的平均值或加权平均值。

# KNN回归示例
set.seed(42)

# 创建回归数据
n <- 200
reg_data <- data.frame(
  x = runif(n, 0, 10),
  y = sin(runif(n, 0, 10)) + rnorm(n, 0, 0.3)
)

# 划分训练集和测试集
train_idx <- sample(1:n, 150)
train_reg <- reg_data[train_idx, ]
test_reg <- reg_data[-train_idx, ]

# 使用caret进行KNN回归
knn_reg <- caret::knnreg(
  y ~ x,
  data = train_reg,
  k = 5
)

# 预测
pred_reg <- predict(knn_reg, test_reg)

# 计算RMSE
rmse <- sqrt(mean((pred_reg - test_reg$y)^2))
print(paste("KNN回归RMSE:", round(rmse, 3)))
## [1] "KNN回归RMSE: 0.896"
# 可视化
test_reg$pred <- pred_reg

ggplot2::ggplot() +
  ggplot2::geom_point(data = train_reg, ggplot2::aes(x = x, y = y), alpha = 0.3) +
  ggplot2::geom_point(data = test_reg, ggplot2::aes(x = x, y = pred), color = "steelblue", size = 2) +
  ggplot2::labs(
    title = "KNN回归",
    subtitle = "蓝点为预测值",
    x = "x", y = "y"
  ) +
  ggplot2::theme_minimal()


7.8 局部加权回归

LOWESS/LOESS

局部加权回归(Locally Weighted Regression)在每个预测点附近拟合一个局部回归模型。

# 局部加权回归示例
set.seed(42)

# 创建数据
n <- 100
lowess_data <- data.frame(
  x = runif(n, 0, 10),
  y = sin(runif(n, 0, 10)) + rnorm(n, 0, 0.2)
)

# 使用loess进行局部加权回归
loess_model <- loess(y ~ x, data = lowess_data, span = 0.5)

# 预测
x_sorted <- sort(lowess_data$x)
y_pred <- predict(loess_model, newdata = data.frame(x = x_sorted))

# 可视化
ggplot2::ggplot(lowess_data, ggplot2::aes(x = x, y = y)) +
  ggplot2::geom_point(alpha = 0.6) +
  ggplot2::geom_line(data = data.frame(x = x_sorted, y = y_pred),
                     ggplot2::aes(x = x, y = y), color = "steelblue", linewidth = 1) +
  ggplot2::labs(
    title = "局部加权回归(LOESS)",
    subtitle = "蓝线为LOESS拟合曲线",
    x = "x", y = "y"
  ) +
  ggplot2::theme_minimal()

span参数的影响

# 比较不同span值
span_values <- c(0.2, 0.5, 0.8)
colors <- c("red", "steelblue", "green")

plot_data <- data.frame()
for (i in seq_along(span_values)) {
  model <- loess(y ~ x, data = lowess_data, span = span_values[i])
  pred <- predict(model, newdata = data.frame(x = x_sorted))
  plot_data <- rbind(plot_data, data.frame(
    x = x_sorted,
    y = pred,
    span = paste("span =", span_values[i])
  ))
}

ggplot2::ggplot() +
  ggplot2::geom_point(data = lowess_data, ggplot2::aes(x = x, y = y), alpha = 0.3) +
  ggplot2::geom_line(data = plot_data, ggplot2::aes(x = x, y = y, color = span), linewidth = 1) +
  ggplot2::labs(
    title = "不同span参数的LOESS拟合",
    x = "x", y = "y", color = "参数"
  ) +
  ggplot2::theme_minimal()


7.9 意义与使用场景

核心意义

KNN和基于实例的方法是懒惰学习的代表,不需要显式的训练过程,直接存储训练数据。

优缺点

优点 缺点
简单直观 预测速度慢
无需训练 需要大量存储
可处理多分类 对特征缩放敏感
可处理非线性 高维数据效果差

适用场景

举例:医学诊断参考

场景:根据患者的症状和检查结果,参考相似病例进行诊断。

为什么用KNN?

医学诊断中,相似病例往往有相似的诊断结果。KNN可以找到历史上最相似的病例,为医生提供参考。医生可以查看这些相似病例的治疗方案和预后,辅助决策。

KNN vs 其他方法

场景 推荐方法
小数据、低维 KNN
大数据、高维 SVM、树模型
需要解释性 决策树
需要概率输出 逻辑回归

第八章 贝叶斯与概率图模型

8.1 朴素贝叶斯分类器

贝叶斯定理

贝叶斯定理描述了条件概率之间的关系:

\[P(A|B) = \frac{P(B|A) \cdot P(A)}{P(B)}\]

在分类问题中:

\[P(y|x) = \frac{P(x|y) \cdot P(y)}{P(x)}\]

朴素贝叶斯假设

朴素贝叶斯假设特征之间相互独立:

\[P(x|y) = \prod_{j=1}^{p} P(x_j|y)\]

分类决策

\[\hat{y} = \arg\max_y P(y) \prod_{j=1}^{p} P(x_j|y)\]

三种朴素贝叶斯模型

模型 特征类型 分布假设
高斯朴素贝叶斯 连续特征 高斯分布
多项式朴素贝叶斯 计数特征 多项式分布
伯努利朴素贝叶斯 二值特征 伯努利分布
# 高斯朴素贝叶斯示例
set.seed(42)

# 使用iris数据集
data(iris)

# 划分训练集和测试集
train_idx <- sample(1:nrow(iris), 0.7 * nrow(iris))
train_iris <- iris[train_idx, ]
test_iris <- iris[-train_idx, ]

# 训练高斯朴素贝叶斯
nb_gaussian <- e1071::naiveBayes(Species ~ ., data = train_iris)

# 预测
nb_pred <- predict(nb_gaussian, test_iris)

# 混淆矩阵
conf_matrix <- table(实际 = test_iris$Species, 预测 = nb_pred)
print(conf_matrix)
##             预测
## 实际       setosa versicolor virginica
##   setosa         12          0         0
##   versicolor      0         14         1
##   virginica       0          1        17
# 准确率
accuracy <- mean(nb_pred == test_iris$Species)
print(paste("高斯朴素贝叶斯准确率:", round(accuracy, 3)))
## [1] "高斯朴素贝叶斯准确率: 0.956"

查看模型参数

# 查看各类别的先验概率
print("先验概率:")
## [1] "先验概率:"
print(nb_gaussian$apriori)
## Y
##     setosa versicolor  virginica 
##         38         35         32
# 查看各特征的条件概率参数(均值和标准差)
print("条件概率参数:")
## [1] "条件概率参数:"
print(nb_gaussian$tables)
## $Sepal.Length
##             Sepal.Length
## Y                [,1]      [,2]
##   setosa     5.015789 0.3405293
##   versicolor 5.874286 0.5135476
##   virginica  6.606250 0.6983564
## 
## $Sepal.Width
##             Sepal.Width
## Y                [,1]      [,2]
##   setosa     3.426316 0.4084805
##   versicolor 2.745714 0.3229850
##   virginica  2.996875 0.3560395
## 
## $Petal.Length
##             Petal.Length
## Y                [,1]      [,2]
##   setosa     1.447368 0.1389864
##   versicolor 4.217143 0.5176222
##   virginica  5.543750 0.5529612
## 
## $Petal.Width
##             Petal.Width
## Y                 [,1]      [,2]
##   setosa     0.2473684 0.1156342
##   versicolor 1.3228571 0.2001260
##   virginica  2.0125000 0.2836968

多项式朴素贝叶斯

# 多项式朴素贝叶斯示例:文本分类
set.seed(42)

# 创建模拟文本数据
n <- 200
text_data <- data.frame(
  word1 = rpois(n, lambda = sample(c(2, 5, 1), n, replace = TRUE)),
  word2 = rpois(n, lambda = sample(c(3, 1, 4), n, replace = TRUE)),
  word3 = rpois(n, lambda = sample(c(1, 3, 2), n, replace = TRUE)),
  category = factor(sample(c("医疗", "体育", "科技"), n, replace = TRUE))
)

# 划分训练集和测试集
train_idx <- sample(1:n, 150)
train_text <- text_data[train_idx, ]
test_text <- text_data[-train_idx, ]

# 训练多项式朴素贝叶斯
nb_multinomial <- e1071::naiveBayes(category ~ ., data = train_text, type = "raw")

# 预测
pred_text <- predict(nb_multinomial, test_text)

# 准确率
accuracy <- mean(pred_text == test_text$category)
print(paste("多项式朴素贝叶斯准确率:", round(accuracy, 3)))
## [1] "多项式朴素贝叶斯准确率: 0.34"

伯努利朴素贝叶斯

# 伯努利朴素贝叶斯示例
set.seed(42)

# 创建二值特征数据
n <- 200
binary_data <- data.frame(
  symptom1 = sample(c(0, 1), n, replace = TRUE),
  symptom2 = sample(c(0, 1), n, replace = TRUE),
  symptom3 = sample(c(0, 1), n, replace = TRUE),
  symptom4 = sample(c(0, 1), n, replace = TRUE),
  disease = factor(sample(c("健康", "患病"), n, replace = TRUE, prob = c(0.7, 0.3)))
)

# 划分数据
train_idx <- sample(1:n, 150)
train_binary <- binary_data[train_idx, ]
test_binary <- binary_data[-train_idx, ]

# 训练伯努利朴素贝叶斯
nb_bernoulli <- e1071::naiveBayes(disease ~ ., data = train_binary)

# 预测
pred_binary <- predict(nb_bernoulli, test_binary)

# 准确率
accuracy <- mean(pred_binary == test_binary$disease)
print(paste("伯努利朴素贝叶斯准确率:", round(accuracy, 3)))
## [1] "伯努利朴素贝叶斯准确率: 0.7"

8.2 高斯判别分析(GDA)

GDA原理

高斯判别分析假设各类别的数据服从多元高斯分布:

\[P(x|y=k) = \frac{1}{(2\pi)^{p/2}|\Sigma_k|^{1/2}} \exp\left(-\frac{1}{2}(x-\mu_k)^T \Sigma_k^{-1} (x-\mu_k)\right)\]

# GDA示例(使用MASS包的lda函数)
set.seed(42)

# 使用iris数据
gda_model <- MASS::lda(Species ~ ., data = train_iris)

# 查看模型
print(gda_model)
## Call:
## lda(Species ~ ., data = train_iris)
## 
## Prior probabilities of groups:
##     setosa versicolor  virginica 
##  0.3619048  0.3333333  0.3047619 
## 
## Group means:
##            Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa         5.015789    3.426316     1.447368   0.2473684
## versicolor     5.874286    2.745714     4.217143   1.3228571
## virginica      6.606250    2.996875     5.543750   2.0125000
## 
## Coefficients of linear discriminants:
##                    LD1        LD2
## Sepal.Length  1.084659 -0.1357075
## Sepal.Width   1.221953 -2.0950103
## Petal.Length -2.431518  0.8233400
## Petal.Width  -2.643108 -2.4996356
## 
## Proportion of trace:
##    LD1    LD2 
## 0.9916 0.0084
# 预测
gda_pred <- predict(gda_model, test_iris)

# 混淆矩阵
table(实际 = test_iris$Species, 预测 = gda_pred$class)
##             预测
## 实际       setosa versicolor virginica
##   setosa         12          0         0
##   versicolor      0         15         0
##   virginica       0          1        17
# 准确率
accuracy <- mean(gda_pred$class == test_iris$Species)
print(paste("GDA准确率:", round(accuracy, 3)))
## [1] "GDA准确率: 0.978"

GDA vs 朴素贝叶斯

特点 GDA 朴素贝叶斯
特征独立性 不假设独立 假设独立
参数数量
数据需求
灵活性

8.3 贝叶斯网络基础

贝叶斯网络概念

贝叶斯网络是一种概率图模型,用有向无环图(DAG)表示变量之间的依赖关系。

组成: - 节点:随机变量 - :依赖关系 - 条件概率表(CPT):每个节点的条件概率分布

# 贝叶斯网络示例(概念演示)
# 使用bnlearn包(如果可用)

if (requireNamespace("bnlearn", quietly = TRUE)) {
  library(bnlearn)
  
  # 创建一个简单的贝叶斯网络结构
  # 假设:吸烟 -> 肺癌 -> 死亡;吸烟 -> 心脏病 -> 死亡
  
  # 定义网络结构
  dag <- model2network("[吸烟][肺癌|吸烟][心脏病|吸烟][死亡|肺癌:心脏病]")
  
  # 可视化网络结构
  graphviz.chart(dag, main = "贝叶斯网络结构示例")
  
} else {
  # 手动绘制贝叶斯网络概念图
  bn_nodes <- data.frame(
    node = c("吸烟", "肺癌", "心脏病", "死亡"),
    x = c(1, 0.5, 1.5, 1),
    y = c(3, 2, 2, 1)
  )
  
  ggplot2::ggplot() +
    ggplot2::geom_segment(ggplot2::aes(x = 1, xend = 0.5, y = 2.8, yend = 2.3), 
                          arrow = ggplot2::arrow(length = ggplot2::unit(0.2, "cm")), color = "gray50") +
    ggplot2::geom_segment(ggplot2::aes(x = 1, xend = 1.5, y = 2.8, yend = 2.3), 
                          arrow = ggplot2::arrow(length = ggplot2::unit(0.2, "cm")), color = "gray50") +
    ggplot2::geom_segment(ggplot2::aes(x = 0.5, xend = 1, y = 1.7, yend = 1.3), 
                          arrow = ggplot2::arrow(length = ggplot2::unit(0.2, "cm")), color = "gray50") +
    ggplot2::geom_segment(ggplot2::aes(x = 1.5, xend = 1, y = 1.7, yend = 1.3), 
                          arrow = ggplot2::arrow(length = ggplot2::unit(0.2, "cm")), color = "gray50") +
    ggplot2::geom_point(data = bn_nodes, ggplot2::aes(x = x, y = y), size = 12, color = "steelblue") +
    ggplot2::geom_text(data = bn_nodes, ggplot2::aes(x = x, y = y, label = node), vjust = -2, size = 4) +
    ggplot2::labs(
      title = "贝叶斯网络结构示例",
      subtitle = "吸烟 -> 肺癌 -> 死亡;吸烟 -> 心脏病 -> 死亡"
    ) +
    ggplot2::theme_void()
}

条件独立性

贝叶斯网络的关键性质是条件独立性

\[P(X_1, X_2, ..., X_n) = \prod_{i=1}^{n} P(X_i | Parents(X_i))\]


8.4 隐马尔可夫模型(HMM)基础概念

HMM基本要素

隐马尔可夫模型用于建模序列数据,包含:

  1. 隐藏状态序列:不可观测的状态
  2. 观测序列:可观测的输出
  3. 三个核心参数
    • 初始概率\(\pi\)
    • 转移概率\(A\)
    • 发射概率\(B\)
# HMM概念示意
# 创建一个简单的疾病状态转移模型

# 状态:健康(H)、疾病(S)
# 观测:正常检查结果(N)、异常检查结果(A)

# 转移概率矩阵
transition_matrix <- matrix(
  c(0.9, 0.1,    # H -> H, H -> S
    0.3, 0.7),   # S -> H, S -> S
  nrow = 2, byrow = TRUE
)
rownames(transition_matrix) <- c("健康", "疾病")
colnames(transition_matrix) <- c("健康", "疾病")

print("转移概率矩阵:")
## [1] "转移概率矩阵:"
print(transition_matrix)
##      健康 疾病
## 健康  0.9  0.1
## 疾病  0.3  0.7
# 发射概率矩阵
emission_matrix <- matrix(
  c(0.8, 0.2,    # H -> N, H -> A
    0.3, 0.7),   # S -> N, S -> A
  nrow = 2, byrow = TRUE
)
rownames(emission_matrix) <- c("健康", "疾病")
colnames(emission_matrix) <- c("正常", "异常")

print("发射概率矩阵:")
## [1] "发射概率矩阵:"
print(emission_matrix)
##      正常 异常
## 健康  0.8  0.2
## 疾病  0.3  0.7
# 初始概率
initial_prob <- c(健康 = 0.9, 疾病 = 0.1)
print("初始概率:")
## [1] "初始概率:"
print(initial_prob)
## 健康 疾病 
##  0.9  0.1

HMM可视化

# 可视化HMM结构
hmm_states <- data.frame(
  state = c("健康", "疾病"),
  x = c(1, 3),
  y = c(2, 2)
)

hmm_obs <- data.frame(
  obs = c("正常", "异常"),
  x = c(1, 3),
  y = c(1, 1)
)

ggplot2::ggplot() +
  ggplot2::geom_segment(ggplot2::aes(x = 1.3, xend = 2.7, y = 2, yend = 2), 
                        arrow = ggplot2::arrow(length = ggplot2::unit(0.2, "cm")), color = "gray50") +
  ggplot2::geom_segment(ggplot2::aes(x = 2.7, xend = 1.3, y = 1.85, yend = 1.85), 
                        arrow = ggplot2::arrow(length = ggplot2::unit(0.2, "cm")), color = "gray50") +
  ggplot2::geom_segment(ggplot2::aes(x = 1, xend = 1, y = 1.7, yend = 1.3), 
                        arrow = ggplot2::arrow(length = ggplot2::unit(0.2, "cm")), color = "steelblue", linetype = "dashed") +
  ggplot2::geom_segment(ggplot2::aes(x = 3, xend = 3, y = 1.7, yend = 1.3), 
                        arrow = ggplot2::arrow(length = ggplot2::unit(0.2, "cm")), color = "steelblue", linetype = "dashed") +
  ggplot2::geom_point(data = hmm_states, ggplot2::aes(x = x, y = y), 
                       size = 12, color = "steelblue", shape = 15) +
  ggplot2::geom_point(data = hmm_obs, ggplot2::aes(x = x, y = y), 
                       size = 10, color = "salmon", shape = 16) +
  ggplot2::geom_text(data = hmm_states, ggplot2::aes(x = x, y = y, label = state), 
                      vjust = -2, size = 4) +
  ggplot2::geom_text(data = hmm_obs, ggplot2::aes(x = x, y = y, label = obs), 
                      vjust = -2, size = 4) +
  ggplot2::annotate("text", x = 2, y = 2.15, label = "转移", size = 3, color = "gray30") +
  ggplot2::annotate("text", x = 1.5, y = 1.5, label = "发射", size = 3, color = "steelblue") +
  ggplot2::annotate("text", x = 3.5, y = 1.5, label = "发射", size = 3, color = "steelblue") +
  ggplot2::labs(
    title = "隐马尔可夫模型结构",
    subtitle = "方形:隐藏状态;圆形:观测值;实线:转移;虚线:发射"
  ) +
  ggplot2::theme_void()

HMM的三个基本问题

  1. 评估问题:给定模型,计算观测序列的概率(前向算法)
  2. 解码问题:给定观测序列,求最可能的状态序列(Viterbi算法)
  3. 学习问题:估计模型参数(Baum-Welch算法)

8.5 高斯过程基础

高斯过程概念

高斯过程(Gaussian Process, GP)是一种非参数贝叶斯方法,定义了函数的分布。

定义:对于任意输入点集合\(\{x_1, ..., x_n\}\),函数值\(\{f(x_1), ..., f(x_n)\}\)服从多元高斯分布。

核心要素: - 均值函数 \(m(x)\) - 协方差函数(核函数) \(k(x, x')\)

# 高斯过程回归示例
set.seed(42)

# 创建数据
n_train <- 20
x_train <- sort(runif(n_train, 0, 10))
y_train <- sin(x_train) + rnorm(n_train, 0, 0.1)

# 测试点
x_test <- seq(0, 10, length.out = 100)

# 定义RBF核函数
rbf_kernel <- function(x1, x2, sigma = 1, length_scale = 1) {
  sigma^2 * exp(-0.5 * (outer(x1, x2, "-") / length_scale)^2)
}

# 计算协方差矩阵
K <- rbf_kernel(x_train, x_train)
K_star <- rbf_kernel(x_test, x_train)
K_star_star <- rbf_kernel(x_test, x_test)

# 添加噪声
sigma_n <- 0.1
K_noise <- K + sigma_n^2 * diag(n_train)

# 高斯过程预测
L <- chol(K_noise)
alpha <- solve(t(L), solve(L, y_train))
y_pred <- K_star %*% alpha

# 计算预测方差
v <- solve(L, t(K_star))
y_var <- diag(K_star_star) - colSums(v^2)
y_sd <- sqrt(pmax(y_var, 0))

# 可视化
gp_data <- data.frame(
  x = x_test,
  y_mean = as.vector(y_pred),
  y_lower = as.vector(y_pred - 1.96 * y_sd),
  y_upper = as.vector(y_pred + 1.96 * y_sd)
)

ggplot2::ggplot() +
  ggplot2::geom_ribbon(data = gp_data, ggplot2::aes(x = x, ymin = y_lower, ymax = y_upper), 
                        fill = "steelblue", alpha = 0.3) +
  ggplot2::geom_line(data = gp_data, ggplot2::aes(x = x, y = y_mean), color = "steelblue", linewidth = 1) +
  ggplot2::geom_point(data = data.frame(x = x_train, y = y_train), 
                       ggplot2::aes(x = x, y = y), size = 2) +
  ggplot2::labs(
    title = "高斯过程回归",
    subtitle = "阴影区域为95%置信区间",
    x = "x", y = "y"
  ) +
  ggplot2::theme_minimal()

核函数的选择

# 比较不同长度尺度
length_scales <- c(0.5, 1, 2)

gp_compare <- data.frame()

for (l in length_scales) {
  K <- rbf_kernel(x_train, x_train, length_scale = l)
  K_star <- rbf_kernel(x_test, x_train, length_scale = l)
  K_star_star <- rbf_kernel(x_test, x_test, length_scale = l)
  
  K_noise <- K + sigma_n^2 * diag(n_train)
  L <- chol(K_noise)
  alpha <- solve(t(L), solve(L, y_train))
  y_pred <- K_star %*% alpha
  
  gp_compare <- rbind(gp_compare, data.frame(
    x = x_test,
    y = as.vector(y_pred),
    length_scale = paste("l =", l)
  ))
}

ggplot2::ggplot() +
  ggplot2::geom_point(data = data.frame(x = x_train, y = y_train), 
                       ggplot2::aes(x = x, y = y), size = 2) +
  ggplot2::geom_line(data = gp_compare, ggplot2::aes(x = x, y = y, color = length_scale), linewidth = 1) +
  ggplot2::labs(
    title = "不同长度尺度的高斯过程",
    x = "x", y = "y", color = "长度尺度"
  ) +
  ggplot2::theme_minimal()

高斯过程的特点

优点 缺点
提供不确定性估计 计算复杂度高\(O(n^3)\)
非参数、灵活 大数据集难以处理
自动处理噪声 核函数选择困难

8.6 意义与使用场景

核心意义

贝叶斯方法提供了概率框架下的机器学习方法,能够:

  1. 量化不确定性:输出概率分布而非点估计
  2. 融入先验知识:通过先验分布引入领域知识
  3. 处理小样本:在数据有限时仍能给出合理估计

各方法适用场景

方法 适用场景
朴素贝叶斯 文本分类、高维稀疏数据
GDA 低维连续数据、类别分离明显
贝叶斯网络 变量间有明确依赖关系
HMM 序列数据、时间序列
高斯过程 小样本回归、不确定性估计

医学应用举例

场景:疾病诊断与预后预测

为什么用贝叶斯方法?

  1. 朴素贝叶斯:根据症状快速筛查疾病
  2. 贝叶斯网络:建模疾病与风险因素之间的因果关系
  3. HMM:分析疾病进展过程,预测患者状态变化
  4. 高斯过程:预测药物剂量-效应曲线,提供置信区间

总结

本教程涵盖了机器学习的核心概念和方法:

第一章:机器学习基础概念,包括学习范式、问题类型、模型评估基础

第二章:数据预处理与特征工程,是机器学习成功的关键步骤

第三章:模型评估与验证,确保模型泛化能力

第四章:线性模型,机器学习的基础

第五章:树模型与集成学习,处理复杂非线性关系

第六章:支持向量机与核方法,高维数据的利器

第七章:最近邻与基于实例的方法,简单直观的非参数方法

第八章:贝叶斯与概率图模型,概率框架下的学习方法

学习建议

  1. 理论与实践结合:每学完一个方法,用真实数据实践
  2. 理解原理:不仅会用包,更要理解背后的数学原理
  3. 比较实验:同一问题尝试不同方法,比较优劣
  4. 关注医学应用:将机器学习方法与医学实际问题结合

进阶学习方向

  • 深度学习与神经网络
  • 强化学习
  • 因果推断
  • 联邦学习与隐私保护
  • 可解释机器学习