Basic Statistic Method

0 Preface

在此记录下遇见的比较常用的或者比较有意思的统计方法。

1 数据处理方法

1 缺失值处理-MICE填补法

MICE(Multiple Imputation by Chained Equations)，数据集中指定的变量的缺失值使用其他变量作为预测变量进行估算。没有相关性的完全独立变量的数据集不会产生准确的插补。其基本步骤是：

• 初始插补：对每个缺失值进行初始插补（通常使用均值、回归或随机采样）。
• 迭代插补：逐变量插补：将缺失的变量视为被预测的目标变量，对于每个缺失的变量基于现有数据构建回归模型预测缺失值；
• 循环迭代：重复步骤（一般为5到10次），逐步更新插补结果，直到收敛（插补值不再发生显著变化）。

MICE可以使用线性回归、逻辑回归、预测均值匹配（PMM）的程序来选择要估算的值。

PMM（预测均值匹配Predictive Mean Matching），核心思想是通过匹配（注意，是匹配）预测值来插补缺失数据。首先构建回归模型，对缺失值通过回归模型计算其预测值，在原本无缺失值（观测数据中的真实值）的数据中找到与预测值最接近的值，这意味着插补的缺失值直接来源于已有数据中的观测值。

即通过其他变量回归进行预测，然后寻找真实数据中最接近的值。并且重复迭代插补，直到插补的值收敛。

2 数据预处理-RIN

Rank-Inverse Normal Transformation：逆秩正态转换，在处理数据分布非正态时，将原始数据的分布通过秩排序和正态变换转换为近似正态分布。

步骤：

• 秩排序：对数据集中每个观测值进行秩排序。假设有$n$个样本，根据观测值$x_i$为观测值本身分配一个秩$R_i$，最小的值对应最小的秩。
• 比例转换：将每个秩转换为一个百分比或累积分布的概率，根据样本的秩转换得到比例值：$p_i=\frac{R_i-0.5}{n}$。分子做了一个偏移，避免取到0或1。
• 逆正态转换：以比例值作为标准正态分布的分位数，得到最后的转换值：$y_i=\Phi^{-1}(p_i)$。

3 匹配对照方法-matchit()函数

matchit()函数属于MatchIt包，专门用于执行各种类型的倾向评分匹配（propensity score matching），允许基于一个或多个变量进行最近邻匹配、全匹配、卡尺匹配、倾向评分匹配。

matchit(
  formula,          # 通过公式的形式指定因变量（case）和自变量（匹配变量）
  data,             
  method = "nearest", # 匹配方法
  distance = "logit", # 估计倾向评分的模型
  caliper = NULL,     # 设定匹配宽限（仅用于连续变量，如年龄、血糖）
  ratio = 1,         # 每个病例匹配的对照数
  exact = NULL,      # 强制精确匹配的变量（如性别、种族）
  replace = FALSE,   # 是否允许有放回匹配
  verbose = TRUE ,   # 显示具体信息
  ...                # 其他参数
)

1. method：“nearest”：最近邻匹配，按举例最近的对照样本进行匹配；“exact”：精确匹配，仅匹配完全相同的值；“full”：全匹配，每个病例匹配所有对照，权重用于校正；，“subclass”：子分类匹配，通过倾向评分将样本分成多个子类；“optimal”：最优匹配，通过最小化总距离来匹配。
2. distance：“logit”：逻辑回归；“probit”；”linear”。
3. caliper：宽限值，限制某些变量的匹配范围；

在匹配成功之后，可以使用love plot（平衡图）来检查匹配的效应，主要观察的变量是各个指标的SMD（标准化均值差异）和VR（方差比）。也可以逐个变量使用显著性统计检验（连续性变量使用t检验，分类型变量使用卡方检验）。

也可以使用倾向得分评分来衡量：Distance：倾向得分：每个样本根据匹配变量计算出的概率，根据匹配公式中的自变量带入一个Logit模型或其他模型估计得出，表示每个样本被分配为“病例”的倾向评分。病例和对照的distance差异越小，越接近于配对。

4 相关系数热力图

对于感兴趣的变量，我们通常绘制相关系数图来表示其相关性，（要注意相关性只是初步的认知，对于低相关性的变量，我们可以认为其存在独立的可能性；对于高相关性的变量，我们需要在后续考虑其高相关性）在相关系数图中，可以使用渐变色、圆形、数字、相关性$p$值来表示。

接下来在此着重介绍绘制存在显著差异$p$值的相关系数热力图：

library(psych)
library(corrplot)

# df 中记载了我们感兴趣的变量的信息，通过columns_to_check向量索引到这些数据，并且提取记录在df_subset中
df_subset <- df[, ..columns_to_check]

# corr.test存储在psych库中，能够计算个体对之间的相关性以及相关显著检验。
corr_dftj <- corr.test(df_subset, use = "pairwise.complete.obs", method = "pearson")
# 提取相关性矩阵，用于后续绘制相关热力图
corr_dftj_matrix <- round(as.matrix(corr_dftj$r), 2)  # 获取相关系数矩阵
# 提取相关显著检验
P <- corr_dftj$p  # 获取 p 值矩阵

corrplot(corr_dftj_matrix, 
         p.mat = P, 
         method = "circle", # 显著星星时最好用圆
         type = "upper",
         col = colorRampPalette(c("#6495ED", "white", "#8B0000"))(200),  
         sig.level = c(.001, .01, .05), # 显著性阈值
         tl.pos = "lt",  # 标签位置
         tl.col = "black", # 标签颜色
         tl.cex = 2, # 标签字体
         tl.offset = 0.2, # 标签离图案间距
         cl.pos = "r", # 图例位置
         insig = "label_sig", # 标记显著性符号 
         pch.cex = 2, # 星星的大小
         pch.col = "red",  # 星星的颜色
         cl.cex = 1.5)    # 图例字体

# 在下三角区域显示数字相关系数
corrplot(corr_dftj_matrix , 
         type = "lower", 
         method = "number", # 相关系数值
         col = "coral4",  # 颜色
         tl.pos = "n",  # 变量
         cl.pos = "n",  # 图例
         number.cex = 1.4, 
         add = TRUE, 
         diag = FALSE)

psych 包中的 corr.test 函数返回的 p 值矩阵是针对相关系数的显著性检验的结果。这些 p 值是通过 t 检验 检验每个相关系数是否显著不同于零。

• 原假设$H_0$：变量之间没有线性关系；
• 备择假设$H_1$：变量之间存在线性关系。

检验统计量 t 统计量：

$t=\frac{r\sqrt{n-2}}{\sqrt{1-r^2}}$

其中 r 是相关系数，n 是样本大小。

• 推导过程：相关系数 r 衡量变量之间的线性相关程度，其通用计算公式为：

$r=\frac{\sum_{i=1}^n(X_i-\bar{X})(Y_i-\bar{Y})}{\sqrt{\sum_{i=1}^n(X_i-\bar{X})^2\sum_{i=1}^n(Y_i-\bar{Y})^2}}$

在原假设成立的情况下：r = 0，此时样本相关系数近似为一个 t 分布（利用到回归）

前提假设：变量X和Y都服从正态分布

2 差异性检验方法

2.1 数值型变量

Wald t-test

Wilcox 秩检验

Two-Sided Unpaired Welch’s t test

该方法应用于PMID: 37353004（Tiannan Guo: 《Proteome Landscapes of Human Hepatocellular Carcinoma and Intrahepatic Cholangiocarcinoma》）

Welch’s t检验是一种用于比较两组独立样本均值是否存在显著差异的统计方法。是独立样本t检验的改进方法，适用于两组方差不相等（违反方差齐性假设时）的情况。

Welch’s t统计量的计算如下：

$t=\frac{\bar{X}_1-\bar{X}_2}{\sqrt{\frac{s_1^2}{n_1}+\frac{s_2^2}{n_2}}}$

计算公式中的符号比较常见，再次就不过多赘述，分别表示样本均值、方差、样本量。该统计量服从的t分布的自由度近似计算如下（Welch-Satterthwaite方程）：

$df=\frac{\left(\frac{s_1^2}{n_1}+\frac{s_2^2}{n_2}\right)^2}{\frac{(\frac{s_1^2}{n_1})^2}{n_1-1}+\frac{(\frac{s_2^2}{n_2})^2}{n_2-1}}$

该自由度通常是非整数。

• 假设检验：

  • 原假设：两组均值相等。
  • 备择假设：两组均值不等。
  • 双侧检验：两组之间是否存在任意方向的差异。

• 适用情况：

  • 样本量不同，方差不齐
  • 两组数据独立，无配对关系
  • 数据服从正态分布

2.2 分类型变量

卡方检验

Fisher exact Test

Fisher精确检验适用于分析两个分类变量之间独立性的小样本数据的统计检验方法。适用于表格中某个变量具有较小的计数的情况。

• $H_0$：Group1和Group2的目标检测变量没有显著差异，即检测变量与组别无关。
• $H_1$：Group1和Group2的目标检测变量存在显著差异，即检测变量与组别相关。

Fisher精确检验依赖于计算每种可能的列联表的概率并累积其概率和。其底层逻辑基于超几何分布和列联表的排列组合，其核心思想是在列联表的所有可能配置中计算每种组合的概率，累加所有等于观察结果或者比观察结果还要极端的概率。因此确保了小样本情况下的结果精确性。

	Group 1	Group 2
1	a	c
0	b	d

考虑总共有$n=a+b+c+d$个个体，其中有$a+c$个个体是目标检测变量取值为1的，$b+d$个个体是取值为0的。给定的如上的$2\times 2$列联表，其概率由超几何分布给出：

$P(a|a+b,a+c,n)=\frac{\binom{a+b}a\binom{c+d}c}{\binom n{a+c}}$

Fisher精确检验的p值是基于计算观察到的列联表及所有比其更极端的表的概率。所谓“更极端”的表是指那些在原假设成立时比当前观察到的表更加偏离独立性的表。偏离独立性意味着更偏向于一种极端（如检出率极低或极高的情况）。因此，Fisher精确检验的p值为所有这些极端表的概率和：

$$p=\sum_\text{tables which are more extreme}P(a\mid a+b,a+c,n)$$。 ### CMH检验 CMH检验（Cochran-Mantel-Haenszel检验）是一种用于分析两个分类变量在不同层次（或分层）中的关系的统计检验方法，特别是在控制潜在混杂变量的情况下。这种检验常用于医学和社会科学研究中，以评估两个分类变量之间的关联是否在不同的层次中保持一致。 CMH检验的主要目的是比较在多个分层中，两个分类变量之间的关系。它通过计算一个综合的关联统计量，来判断在所有层次中是否存在显著的关系。 - 假设：原假设$H_0$：所有层次中，两个分类变量之间没有关联（条件独立）；备择假设$H_1$：至少有一个层次中，两个分类变量之间存在关联。 - 数据结构：CMH检验适用于具有分层结构的$2\times2$的列联表。通常情况下，数据会被组织成一个三维表格，其中每个层次对应一个$2\times2$列联表。 - 检测步骤： 1. 对于每个层次，构建一个$2\times2$列联表，显示两个分类变量的频数。 2. 计算每个层次的边际总和（列和和行和）。 3. 计算CMH统计量：$Q=\frac{(N\cdot(\mathrm{ad-bc})^2)}{\sum(a+b)(c+d)(a+c)(b+d)}$ 4. CMH统计量通常服从自由度为1的卡方分布。 只适用于两个分类变量之间的关系分析，数据是分层的。如果层次之间的样本量差异过大，会影响结果的稳定性。 ### Kruskal-Wallis H 检验 Kruskal-Wallis H 检验是一种非参数统计检验，用于**比较三个或多个独立组的中位数是否存在显著差异**。它不依赖于正态分布的假设，非常适合于样本量较小或数据不符合正态分布的情况。 - 原理：对所有观测值进行排名，比较这些排名在各组之间的分布情况。假设有$k$组独立样本，每组的样本量分别为$n_1,n_2,...,n_k$ - 假设：原假设$H_0$：所有组的中位数相等；备择假设$H_1$：至少有一组的中位数与其他组不同。 - 检测步骤： 1. 将所有观测值合并后进行排序，并为每个观测值分配一个排名。如果有相同值，则分配它们的平均排名。 2. 计算组的排名和：计算每个组的排名和$R_j,j=1,2,...,k$。 3. 计算检验统计量：$H=\frac{12}{N(N+1)}\sum_{j=1}^k\frac{R_j^2}{n_j}-3(N+1)$，$N$表示样本总数，$n_j$表示$j$组中的样本量。 4. 根据计算得到的统计量以及$k-1$自由度，根据卡方分布得到$p$值。 5. 做出决策。 ## 2.3 数值/分类型变量 ### 基于CLR的系数显著性检验 ## 2.4 趋势检验（P-trend） 趋势检验用于检测不同组间某种排序或分类的单调变化趋势。例如考察某种变量（比如剂量、时间或其他时间因素）对结果变量的单调影响。常见的方法包括Jonckheere-Terpstra检验和卡方检验。 ### Jonckheere-Terpstra检验（Jonckheere-Terpstra Test） 一种非参数趋势检验，适用于检测一个分组变量（有序分组）在数值或顺序响应变量上的单调趋势，这个检验尤其适用于连续或顺序数据。 适用于：1、自变量是有序的分类变量（低、中、高剂量；近、中、远发病时间）；2、因变量是连续变量或可以排序的变量；3、检测因变量在自变量水平的增加中是否出现单调趋势。 该检验基于秩和检验，假设组间存在一个特定的顺序方向。记目前有$k$组数据，组编号按有序变量的水平顺序排列，例如$Group_1,Group_2,...,Group_k$，每组数据量分别是$n_1,n_2,...,n_k$。对所有样本数据计算其秩（大小排序的顺位），对任意两组$G_i$和$G_j$（$i<j$）进行比较，统计从较低组别中取值较大的次数。定义统计量$J$为满足$x_{i,m}<x_{j,n}$的次数，$x_{i,m}$表示第$i$组中的第$m$个样本值。统计量的期望和方差可以通过分布来近似。根据样本量的大小，计算近似正态分布的统计量$Z$，然后计算对应的$p$值。 ### Chi-Square Test 卡方趋势检验常用于检测分组变量（如有序类别）对结果变量（通常是分类数据）的单调关系。卡方趋势检验对数据的期望是趋势单调的，但不必满足严格的线性关系。 - 自变量是有序的分类变量。 - 因变量是二分类或多分类变量。 - 目标是检测因变量的分类比例随自变量的水平增加或减少。 ### 2. 检验过程 卡方趋势检验可以通过构造趋势卡方统计量来评估趋势显著性。具体步骤如下： 1. **构建列联表**： - 将数据按自变量的水平和因变量的类别构建一个二维列联表。例如，行表示自变量的不同水平，列表示因变量的不同类别。 2. **计算每个单元格的期望频数**： - 对于每个单元格，根据行和列的边际总和以及总体样本量，计算期望频数$E_{ij}=\frac{行总和 \times 列综合}{总样本数}$ **计算卡方统计量**： - 卡方统计量可以表示为：$\chi^2=\sum\frac{(O_{ij}-E_{ij})^2}{E_{ij}}$，其中$O_{ij}$为观察到的频数，$E_{ij}$为期望的频数。 - 定义线性趋势卡方统计量：$\chi_{\mathrm{trend}}^2=\frac{\left(\sum_{i=1}^kx_i(O_i-E_i)\right)^2}{\sum_{i=1}^kx_i^2E_i}\sim\chi^2(1)$，其中$x_i$为自变量的有序编码，$O_i$为观察到的频数总和，$E_i$为期望的频数总和。 # 3 降维及聚类方法 ## 3.1 PCOA（主坐标分析） PCA（Principal component analysis）与PCoA（Principal coordianate analysis, 主坐标分析） PCA：利用正交变换把数据从线性相关变量转换为少数几个方差最大的线性无关变量，每个PC是对变量的线性组合，PCA没有要求变量之间线性无关，在进行PCA之前需要对变量进行标准化，统一量纲。 PCA要求的是最后得到的主成分之间是线性无关的，研究者希望在进行PCA之后得到的PC（主成分）能够解释较大的方差，这样的主成分是具有代表性的，是能够达到降维的目的的。 相比于PCA，PCoA并不是想尽可能多地保留方差，而是尽可能地保持数据之间的距离结构，并将其在欧氏空间中进行可视化。PCaA基于距离矩阵而不是基于原始数据矩阵来进行分析，适用于欧氏距离、布雷-柯蒂斯距离、Jaccard距离。 1. 计算距离矩阵 给定一组样本点$x_1,x_2,...,x_n$，计算样本对之间的距离，得到一个对称的$x\times x$距离矩阵$D$，其中$d_{ij}$表示样本$x_i$和$x_j$之间的距离。 2. 双中心化距离矩阵 将距离矩阵$D$进行双中心化，生成一个新的矩阵，用于将数据投射到欧氏空间中，对距离$d_{ij}$，转化为$a_{ij}=\frac{1}{2}d_{ij}^2$，在对矩阵$A$进行中心化，得到矩阵$\Delta_1=(\delta_{ij})$，其中$\delta_{ij}=a_{ij}-\bar{a}_i-\bar{a}_j+\bar{a}$，再对$\Delta_1$进行特征分解，计算器特征值，以及特征向量矩阵。 3. 可视化 选取前k个特征值以及对应的特征向量取值可以对对象进行降维，可以通过二维或三维图展示数据，样本间的距离关系可以得到好的保持。 PCOA能够对簇进行检验，检验两簇是否在统计意义上是相隔开的。 ### PCoA Plus: Permutation Test PERMANOVA是多元方差分析的非参数变体，被用来比较多组检测样本的统计指标值的异同。其利用**距离矩阵**对总方差进行分解，分析不同分组因素对样本差异的解释度，并且使用**置换检验**对每个变量解释的统计学意义进行显著性分析。有时也被称为Adonis检验。 其原假设为： - 每组样品在其检测指标构成的检测空间中的中心点和离散度没有差别。 拒绝原假设，即表示不同样本的分布显著不同。为此，该检验是根据距离矩阵完成的，即要预先计算样本在检测指标定义的多维空间的距离。 PERMANOVA对检测数据的分布没有限制，不受组间协方差的影响，其前提假设是： 1. 每个对象的数据可交换； 2. 可交换的对象相互独立； 3. 每个样本的数据具有一致的多变量分布。 检验统计值计算的是伪F值，计算方式是不同样本之间的距离平方和除以同一组样本之间距离的平方和。$F=\frac{SS_A/(a-1)}{SS_W/(N-a)}$。该统计量的显著性通过数据置换的方式计算，比较随机置换数据获得的F值，是否高于或者等于实际观测到的F值，如果多于5%随机置换得到的F值高于实际观测值，则记为p>0.05，即不同组的样本之间不存在显著差异。 ## 3.2 UMAP 应用场景：非线性、复杂、高维 Uniform Manifold Approximation and Projection是一种非线性降维技术， - Projection：通过投影点在平面、曲面或线上再现空间对象的过程或技术，即：将对象从高维空间映射到低维空间的映射。 - Approximation：假设只有一组有限的样本点，而不是构成流形的整个集合。 - Manifold：流形是一个拓扑空间，一维流形包括线和圆；二维流形包括平面、球体、环形等。 - Uniform：均匀性假设：数据样本均匀分布在流形上，但是在现实中很难实现，因此这个假设引出了在流形上距离是变化的概念：空间本身扭曲，空间根据数据显得更稀疏或更密集的位置进行拉伸收缩。 即UMAP：假设可用数据**均匀分布**在**拓扑空间**中，可以根据这些**有限数据近似**并映射到低维空间。 **步骤**： 1. 高维数据中的近邻关系建模，构建近邻图，通过**K-NN**算法查找每个数据点的局部邻域，使用模糊集合表示相邻点之间的关系权重，权重由**核函数**计算。例如高斯核函数： $w(x_i,x_j)=\exp\left(-\frac{\mathrm{dist}(x_i,x_j)-\rho_i}{\sigma_i}\right)$ 其中$\rho_i$表示点$x_i$的最近距离的阈值，以确保每个点都有一定的邻居；$\sigma_i$是一个尺度参数，用于调整权重分布，最终结果会生成一个加权图，节点表示数据点，边的权重表示邻域的相似度。 2. 定义低维空间中的目标分布，低维空间中同样使用加权图建模点之间的关系：$v(y_i,y_j)=\left(1+a\cdot\|y_i-y_j\|^{2b}\right)^{-1}$其中a和b是超参数，用于控制权重函数的性状。 3. 优化嵌入：通过最大化高维和低维空间的相似性，优化目标函数：$\mathcal{L}=\sum_{i,j}\left[w(x_i,x_j)\log(v(y_i,y_j))+(1-w(x_i,x_j))\log(1-v(y_i,y_j))\right]$目标是使高维空间中的点对权重与低维空间中的权重尽可能匹配，这一优化过程通常通过随机梯度下降完成。 主要参数： 1. **n_neighbors**（邻居数量）：决定近邻图的规模，影响局部结构的保留程度。较小值倾向于保留局部细节，较大值则关注全局结构。 2. **min_dist**（最小距离）：控制低维嵌入点之间的最小距离。较小值会导致点聚集，较大值则使点分布更加稀疏。 3. **metric**（距离度量）：用于高维空间中计算点间距离。支持多种距离度量，如欧几里得距离、余弦距离等。 4. **a, b**：控制低维空间中权重分布的形状，通常通过算法内部自动选择。 ## 3.3 WGCNA Weighted Gene Co-Expression Network Analysis是一种基于基因表达数据构建加权共表达网络的方法。 - 将基因作为网络节点，通过它们在不同样本中表达水平的相似度（一般用相关系数）来构建网络 - 通过设置软阈值，将基因之间的相似度转换为相应的权重，从而构建符合生物学规律的无尺度网咯 - 通过层次聚类或动态剪切方法，将表达模式相似的一组基因归纳成“模块”，每个模块通常聚集在一起，表达具有相同的功能、调控或生物学意义。随后，可以把网络模块与感兴趣的表型进行关联分析，发现与表现最为相关的核心模块，进一步筛选出关键基因。 该分析方法的目的是寻找协同表达的基因模块（module），探索基因网络与关注的表型之间的关联关系，以及网络中的核心基因。 #### 方法学 WGCNA分为表达量聚类分析和表型关联两部分，包括基因之间相关系数计算，基因模块的确定，共表达网络，模块与形状关联四个步骤。 1. 数据预处理：检查数据的质量，去除无意义的噪音数据。 2. 表达量聚类分析：任意计算两个基因之间的相关系数。为了衡量两个基因是否具有相似的表达模式，设置阈值后进行筛选，高于阈值则认为是相似的。WGCNA分析采用**相关系数加权值，对基因相关系数取N次幂**，使得网络中的基因之间的连接服从**无尺度网络分布（scale-free networks）**。这里的N被称为软阈值（soft-thresholding power），对于强相关的基因，经过N次幂计算之后会增强其相关性，对于弱相关的基因则会减少其贡献，使网络更符合无尺度网络的特性。 - 无尺度网络（Scale-Free Network）是一种网络拓扑结构，在许多自然系统和生物网络（如基因调控网络、蛋白质相互作用网络等）中广泛存在。它的主要特点是**节点的连接度分布遵循幂律分布（Power-Law Distribution）**，即大多数节点只有少数连接，而少数关键节点（**枢纽节点，Hub Nodes**）具有大量连接。 3. 计算拓扑重叠矩阵（TOM）：减少噪音影响，增强模块化结构。计算拓扑重叠度衡量基因间的相似性，不仅考虑基因之间的直接连接，还考虑他们之间是否有相同的邻居，然后计算基于拓扑重叠度的距离矩阵，使高度拓扑重叠的基因聚类到同一模块。 4. **识别基因模块（层次聚类）**：使用动态剪切树算法识别基因模块，如果两个模块的特征基因的相关性超过一定阈值则可以合并。至此则生成了一个基因共表达网络。 5. 计算模块特征基因（ME）：代表每个基因模块的表达模式，用于后续的表型关联分析。模块特征基因是模块中所有基因表达矩阵的第一主成分，代表了该模块的整体表达模式。计算基因与模块的相关性（Module Membership），相关性越大表示该基因对该模块的代表性越强。 6. 关联模块与表型：探索哪些基因模块与感兴趣的生物表型（如疾病状态、生存时间、临床指标）相关，通过**计算模块特征基因与表型的相关性系数**，得到关联的显著性模块。 7. 识别关键基因（Hub Genes）：筛选出可能在生物过程中起关键作用的基因，关键基因通常是有高模块成员关系（Module Membership，在模块内部起到核心作用）；高基因重要性（Gene Significance，与生物表型显著相关）；高连接度（kWithin，在模块中具有较强的连接性）。 ### 应用 具体参考R包。 # 4 统计模型 ## 4.1 条件逻辑回归 Breslow NE, Day NE. Statistical methods in cancer research. Volume I - The analysis of case-control studies. IARC Sci Publ. 1980;(32):5-338. PMID: 7216345. 似然函数能够描述前瞻性分析和病例对照实验。 当在研究中不引入条件（无条件回归模型）时：即使用无条件回归模型来估计相对风险时，会导致估计参数时，参数阶数（在估计时会引入每个分层的参数来描述不同分层之间的差异，这些是干扰参数，也是无关参数）与观察数接近，导致估计发生偏倚。 假设每个分层都由一个匹配的病例-对照构成，暴露记为$x=1$，未暴露记为$x=0$。每对的结果都可以表示为一个$2\times 2$表的形式，拟合的模型形式为： $$\mathrm{pr_i(y=1|x)=\frac{\exp(\alpha_i+\beta x)}{1+\exp(\alpha_i+\beta x)},}$$

而对于每个分层，可能存在四种情况（case暴露，control暴露；case暴露，control不暴露；case不暴露，control暴露；case不暴露，control不暴露）最终无条件模型与条件模型得到的估计会不同，无条件模型相对会低估结果。（书中有详细推导，在此不展开）

为此实际实验中通常设定一个case与一个或多个control进行匹配，此时条件似然函数被简化。

对于条件逻辑回归，一般来说一个病例匹配3~4个对照是比较合理的选择，虽然病例尽可能多地匹配对照，但是在增加到一定数目之后，边际效益会递减。

如下为个人观点：

为什么在设定1：4的匹配，与设定1：10的匹配的蛋白位点不同？

很直观的结果来看，当我们选择1：4匹配时，如果真正的阳性蛋白位点，其case普遍比control要高，那么其p值会很小。而当我们选择1：10匹配时，这一蛋白位点上，出现了较多的control其表达量反而比case要低，此时就降低了该蛋白位点的p值，使得该蛋白位点的阳性被“隐藏“。此时，在1：10匹配时，我们没有考虑完全所有的混杂因素，导致携带更多混杂的control进入了模型的系数估计的过程中。

但这并不是绝对的，我们如何保证1：4的匹配绝对控制混杂了呢？

• 在匹配过程中，要尽量匹配的变量需要满足：其对自变量X产生较大影响 或 对目标变量Y和X都产生影响的变量；
• 在建模过程中，要尽量校正的变量需要满足：对X影响较小（因为这会产生多重共线性），对Y影响较大的变量。

4.2 Cox比例风险模型

4.2.1 普通Cox比例风险模型

Cox模型不仅可以包含类别变量还可以包含数值变量。

1、生存函数：表示对象在时间点t仍然生存的概率，记作$S(t)$。

2、风险函数：表示对象的生存时间达到t后失败的概率，记作$h(t)$。

概率密度函数记作$f(t),h(t)=f(t)/S(t)$，累积分布函数满足$F(t)=1-S(t)$，因此累计风险函数$H(t)=-logS(t)$。

对于Cox模型，存在

$h(t,x)=h_0(t)exp(\beta X)$

即：$ln\frac{h(t,X)}{h_0(t)}=\beta X$

其中$X=(X_1,…,X_n)’$是独立于时间的可能影响生存时间的协变量。$h_0(t)$是一个非参数的基线风险函数，表示协变量全为0时的风险，对所有个体来说都是一样的。并且可知，Cox回归类似于线性回归，Cox回归的是个体风险与人群基线风险比值的对数。

前提条件：风险比值对数与协变量之间呈线性关系、风险比值对数与时间无关。

Tips：如果风险比值对数与时间有关，则需要使用依时协变量模型。

模型检验假设：1、比例风险假设检验：Schoenfeld残差检验，该残差应当与时间无关；若有证据证明某变量违反了比例风险假设，可以考虑将其设置为分层变量。2、线性关系假设：Martingale残差检验；自然样条检验。如违反线性假设，需要考虑对其进行非线性转换。

另：估计协变量的系数时，可以不涉及基线风险系数的估计。而如果涉及到后续预测，则需要估计基线风险，常用的方法是Breslow方法。

4.2.2 时依以及分层Cox比例风险模型

时依（协）变量（Time-dependent covariables）：该变量随时间的改变而发生变化。

当不满足等比例风险时，1、按照协变量分层分析，即采用分层的Cox回归；2、采用时依Cox回归。

4.2.3 加权Cox回归

2018, Daniela Dunkler, et al, Journal of Statistical Software.

在现实情形下，一般很难满足等比例风险假设，在模型构建时忽略等比例风险假设，最终将得到的估计值汇报为”平均效应“，这样就会导致在随访的某些情况下，估计值会高估真实的时间依赖效应，而在其他时段则会低估。当等比例假设不满足时，可以通过分层Cox回归、扩展Cox回归、加权Cox回归进行处理。

• 当只有一个分类协变量（类别较少）并且该协变量不是研究的主要变量时可以使用分层方法；
• 如果一个主要关注的协变量在短期和长期有不同的效应，可以通过将随访时间划分为多个阶段，并在每段中估计分段常数HR来表达其效应。但这种分析方法假设：每个时间段内HR恒定，两个时间段之间HR发生突变——这两种假设通常都不现实
• 扩展Cox回归模型通常可以包含多个协变量与时间的交互项，这种方法虽然灵活但需要较大的样本量。

对应的Schemper等人于2009年提出，在部分似然函数中，对每个事件时刻的贡献加权，加权因子为：”该时刻的生存函数值 × 该时刻随访累计概率的倒数“。因此，在每个事件时刻，其加权函数与”如果没有删失，理论上应该还在风险集中的人数“成正比。通过应用这些权重，加权Cox回归可以估计出平均HR，这个HR可以近似地表示”协调优势比“。

协调优势比：Odds of concordance, OC

$OC=\mathsf{P}(T_A<T_B)/\mathsf{P}(T_B<T_A)$

其中$T_A, T_B$表示从两个治疗组中随机选出的两个个体的生存时间。进而协调优势比可以转换为协调概率

$c=\mathsf{P}(T_A<T_B)=OC/(OC+1)$

这种方法能在不依赖比例风险假设的前提下回答哪个治疗组或因子水平更好。

具体来说，对对数似然函数加权如下：

$\frac{\partial\log L(\beta)}{\partial\beta_r}=\sum_{j=1}^mw(t_j)\left[x_{jr}-\frac{\sum_{l\in R_j}x_{lr}\exp(\hat{\beta}^\top x_l)}{\sum_{l\in R_j}\exp(\hat{\beta}^\top x_l)}\right]=0.$

其中N表示个体个数，$t_j$表示第j个时间点，一共有m个时间点，对于每个个体i都拥有一个行向量$x_i$其中包含k个元素代表着存在k个协变量，$R_j$表示$t_j$时刻时的风险集，$w(t_j)$表示加权函数，$\beta$表示回归系数。

• $w(t_j)=1$则退化为经典的Cox回归;
• $w(t_j)=G(t_j)^{-1}$表示平均回归效应；2000, Xu, O’Quigley, Biostatistics
• $w(t_j)=S(t_j)G(t_j)^{-1}$表示平均风险比；2009, Schemper, Statistics in Medicine

另外，在使用coxphw()进行加权Cox回归时，还能够设置使用Lin-Wei稳健估计方法用于显著性检验，或设置使用Jackknife方法估计协方差矩阵；并且在加权时能够设置对极端的权重进行截断以避免某些事件时刻对结果的过度影响。

4.2.4 解读Cox结果

1. 回归系数$\beta=log(HR)$：对数风险比，Cox模型中直接估计的参数，表示协变量每单位增加时，风险函数在对数尺度上的变化；
2. 风险比HR (Hazard Ratio), $HR=exp(\beta)$：表示暴露组与对照组的即时死亡（事件）风险之比，HR>1表示暴露风险增加，反之即反。
3. 置信区间CI：表示估计的不确定性，95%CI表示在95%的样本重复中，这个区间包含了真实的HR。

4.3 Meta分析

Meta-Analysis用于综合多个独立研究的结果，汇总不同研究的样本量和效应大小。Meta分析是一种定量研究方法，目的是提高统计效能、检测异质性（识别不同研究之间的差异，分析差异来源）、探寻亚组效应（探究变量如何影响总体效应）。

具体步骤如下：

• 明确研究问题与定义效应指标：在明确研究问题时通常涉及某种干预措施、暴露因素或生物标志物对某个结局的影响；效应指标通常为相对风险（RR）、比值比（OR）、均数差（MD）或标准化均数差（SMD）。
• 文献检索与研究选择：系统检索现有的文献；设定明确的纳入标准和排除标准，以确保选取的研究与研究问题相关。通常会考虑实验设计、暴露因素、结局指标等方面；最终筛选文献是符合标准的高质量研究。
• 数据提取与质量评估：从纳入的研究中提取数据，如效应大小、标准误、样本量、随访时间等；再使用适当的工具评估研究的质量，包括随机化、盲法、失访率等，评估偏倚风险。
• 选择适当的Meta分析模型：固定效应模型or随机效应模型。
• 效应合并与异质性分析：通过统计学方法将不同研究的效应大小进行加权平均，得到总体效应值，加权通常基于每个研究的样本量或标准误；异质性分析用来分析各个研究之间的差异是否显著，包括Q检验、$I^2$统计量。（异质性显著时建议采用随即效用模型，否则可以选择固定效应模型）
• 亚组分析与敏感性分析：将数据按照某些特征进行分组，探索不同亚组间的效应差异；敏感性分析用于评估某些特定研究对总体效应估计的影响。
• 发表偏倚评估：一般来说只有显著结果的研究更容易发表，没有显著结果的研究可能会被遗漏，通过漏斗图或Egger回归检验来评估是否存在发表偏倚。
• 报告与解释结果：最终结果通常以森林图展示，展示所有研究的效应大小和置信区间。

4.3.1 固定效应Meta分析

• 基本假设—唯一真实效应：假设所有纳入的研究都在估计相同的真实效应值，研究间的差异完全由样本内的随机误差引起，意味着每个研究实际上都在估计同一效应值。

• 权重分配：通常样本量大的研究会赋予更大的权重。加权方式是通过逆方差加权法来分配的，加权公式是：

  $w_i=\frac{1}{SE_i^2}$，其中$SE_i^2$是该研究效应估计的标准误。

• 效应合并：将个研究的效应估计值按其权重加权平均，得出合并效应值，

  $\hat\theta=\frac{\sum{w_{i}\theta_i}}{\sum{w_{i}}}$

适用于异质性较小的情况，假设所有研究估计同一个效应。

4.3.2 随机效应Meta分析

• 基本假设—多重真实效应：假设每个研究可能估计不同的真实效应值。意味着研究之间存在一定的变异，且这种变异不仅仅由随机误差引起，可能因为研究设计、对象、干预措施等方面的不同；研究之间的差异来自每个研究内的随机误差或研究间的真实效应值差异。

• 权重分配：随机效应模型还考虑了研究间的变异，

  $w_i=\frac{1}{SE^2_i+\tau^2}$

  其中$\tau^2$是研究间的变异，当研究之间差异较大时，会降低各研究的权重。

• 效应合并，同样根据加权合并。

• 异质性估计：DerSimonian-Laird法衡量研究间差异，$I^2$统计量。

适用于异质性较大的情况。

4.3.3 Meta分析的原理

Meta分析的理论基础主要建立在以下几个统计学原理上：

1. 效应大小的汇总
   Meta分析的一个关键思想是将多个独立研究的效应大小进行加权平均。每个研究的权重通常与其样本量成正比（样本量越大，研究结果越可靠，其权重越大）。具体的加权方式取决于是否使用固定效应模型或随机效应模型：
2. 异质性与一致性
   Meta分析假设各个研究所估计的效应大小在某种程度上是可比的，但实际上不同研究可能因为方法、样本、时间、地域等原因存在差异，这些差异就表现为异质性。异质性的存在是对结果一致性的挑战，因此在Meta分析中要首先对异质性进行量化和解释。
3. 统计功效与精度
   Meta分析通过汇总多个研究的数据来提高统计功效（statistical power），即提高发现显著效应的能力。这也是Meta分析的重要优势之一，因为单一研究的样本量有限，可能导致统计功效不足，而Meta分析能够综合多项研究，从而提高检测效应的精度和可靠性。
4. 偏倚与校正
   在进行Meta分析时，研究人员需要警惕各种偏倚。发表偏倚是其中最重要的一种。通过多种技术（如漏斗图、对未发表研究进行估计等），Meta分析可以对偏倚进行检验并作出相应调整。

4.3.4 数学原理

1 效应大小

在Meta分析中，每个独立研究通常都会给出一个效应大小(effect size)。这些效应大小可以是均值差 (Cohen’s d)、标准化均差、相关系数、优势比 (Odds Ratio, OR) 等。假设我们有$k$个独立研究，每个研究给出了一个效应大小$\hat{\theta}_i$和相应的标准误差$\hat{\sigma}_i$。

Meta分析的主要目标是对这些效应大小进行加权平均，从而获得总体效应大小的估计值$\hat{\theta}$。在加权平
均时，每个研究的权重$w_i$通常设为其效应估计精度的倒数平方，即：

$$w_i=\frac1{\hat{\sigma}_i^2}$$

总的效应大小$\hat{\theta}$ 通过加权平均计算：

$$\hat{\theta}=\frac{\sum_{i=1}^kw_i\hat{\theta}_i}{\sum_{i=1}^kw_i}$$

加权平均后，总效应大小的标准误差$\hat{\sigma}_{\hat{\theta}}$则为：

$$\hat{\sigma}_{\hat{\theta}}=\sqrt{\frac1{\sum_{i=1}^kw_i}}$$

2 异质性检验

Meta分析的一个关键步骤是检验各研究之间是否存在显著异质性（heterogeneity）。即，不同研究的效应大小是否相似。如果存在显著异质性，则可能表明各研究的效应大小来源于不同的总体。

1. Cochran’s Q检验

   Cochran’s Q检验用于判断是否存在异质性，其公式如下：

   $$Q=\sum_{i=1}^kw_i(\hat{\theta}_i-\hat{\theta})^2$$

   该统计量服从$\chi^2$分布，自由度为$k-1$。如果$Q$值显著大于$\chi^{2}$分布的临界值，则表明存在异质性。

2. $I^{2}$统计量

$$I^2=\frac{Q-(k-1)}{Q}\times100\%$$

其中，$k$是研究的数量。如果$I^{2}$较大 (通常超过 50%),则表明异质性显著。

3 效应模型

略。

4 偏倚检验

Meta分析需要考虑偏倚，特别是发表偏倚（publication bias），即倾向于发表有显著结果的研究。偏倚会导致Meta分析结果失真。

1. 漏斗图(Funnel Plot)

漏斗图用于可视化偏倚。如果研究结果不存在偏倚，效应大小与标准误的散点图应呈对称的倒漏斗形状。如果存在发表偏倚，图形会呈现不对称性。

1. Egger回归检验

Egger回归检验是一种定量检验方法。通过线性回归分析漏斗图中的不对称性，检验是否存在发表偏倚。其回归模型为：

$$Y_i=\beta_0+\beta_1X_i+\epsilon_i$$

其中，$Y_i$是效应大小标准误的倒数，$X_i$是效应大小。显著的$\beta_1$表明存在偏倚。

4.4 LASSO回归

LASSO回归是加入L1正则项以起到压缩变量前系数（自动化变量选择）的目的。LASSO回归的目标是最小化损失函数：

$\mathrm{Loss}=\mathrm{RSS}+\lambda\sum_{j=1}^p|\beta_j|$

其中RSS表示残差平方和：$\mathrm{RSS}=\sum_{i=1}^n(y_i-\hat{y_i})^2$。

1、使用交叉验证的方法得到最好的Lambda，lambda的遍历是对数间隔（在对数尺度上是均匀分布的）具体来说$\lambda_i=\lambda_{\max}\times\left(\frac{\lambda_{\min}}{\lambda_{\max}}\right)^{\frac{i}{100-1}}\quad\text{where}\quad i=0,1,2,…,99$ 而其中lambda的最大值是根据标准化后的数据自动计算的，该最大值能够将所有的系数压缩为0，表示最严格的惩罚项。（正则化效应的变化在数值范围上通常是非线性的，对数尺度能够更好地表示正则化参数的影响。对于较小的 lambda 值，模型的变化更加显著，而对于较大的 lambda 值，模型的变化相对较小。因此，在对数尺度上分布的 lambda 值能够更均匀地覆盖模型的行为）

2、衡量标准可以是均方误差MSE，也可以是AUC，也可以是偏差Deviance

3、交叉验证的方法会输出lambda_min和lambda_1se：lambda_1se是lambda_min的一个更保守版本，使得验证误差在最优lambda_min 的 标准误差范围内。使得模型更加简单，避免了过拟合，可能会有较大的偏差，但通常具有较好的泛化能力。

4.5 ANCOVA

ANOCOVA, Analysis of Covariance, 协方差分析：结合了方差分析和回归分析的方法，主要用于比较两个或者多个组别的均值，同时控制一个或者多个协变量的影响。

• 核心思想：通过回归调整消除混杂因素的影响，使得不同组之间的均值比较更准确。
• 数学模型：$Y_{ij}=\mu+\tau_i+\beta X_{ij}+\epsilon_{ij}$，i表示组别，j表示相应的个体。$\mu$表示总体的均值，$\tau_i$表示第i组的固定效应，$X_{ij}$表示协变量，假设误差项独立同分布于$N(0,\sigma^2)$。根据这个模型可知，因变量受到分类变量以及协变量的影响。

步骤：

1. 检验协变量与因变量的关系，在进行协方差分析之前，需要检查协变量X是否与因变量Y相关，这意味着协方差分析中的混杂因素应当是与因变量具有相关性的。方法：回归分析、相关系数等。
2. 检验协变量与因变量之间的同质性：协方差分析假设不同组别之间的协变量的回归系数是相同的（平行性假设）。这意味着协变量X对因变量Y的影响在不同组别之间是一致的，换句话说，协变量对因变量的影响不受组别的划分发生变化。方法：引入协变量和组别的交互项进行检验：$Y=\mu+\tau+\beta X + \gamma(X \times Group)+\epsilon$，如果交互项前的系数显著，则说明不满足平行性假设。
3. 回归调整：利用原回归模型估计因变量和协变量之间的关系，并且据此调整因变量（以此减弱协变量对对最终均值比较的影响，这也是这个方法的核心所在）。
4. 方差分析：在调整后的因变量上执行传统的单因素方差分析，检验不同组别之间的均值是否有显著差异。

一些核心假设：

1. 因变量与协变量之间线性相关；
2. 协变量在不同组别之间的均值相近（t检验）；
3. 协变量的回归斜率在各组别之间相同（交互项斜率显著性检验）；
4. 误差项符合正态分布（Q-Q图胡总Shapiro-Wilk检验）；
5. 误差项方差齐性（Levene检验）；
6. 观测值独立（研究设计）。

4.6 限制性立方样条

Restricted Cubic Spline, RCS：是一种非参数回归方法，用于分析连续变量与结局之间的非线性关系。

(To be continued)

4.7 Kaplan-Meier生存曲线估计

乘积极限法，适用于小样本数据。

• 生存函数：Survival Function: $S(t)=P(T>t)$，表示个体存活时间超过时刻t的概率。

其思想是将整个观察时间轴按照事件发生的时间点划分为若干个时间段；在每个事件发生的时间点计算存活概率；将这些存活概率连乘得到在该时间点的累计生存概率。因此KM方法估计得到的生存函数为：

$\hat{S}(t)=\prod_{t_i\leq t}\left(1-\frac{d_i}{n_i}\right)$

其中，$t_i$表示第i个事件发生的时刻，$d_i$表示在时间$t_i$时有多少事件发生，$n_i$表示在该时刻有多少个个体仍然处于风险中（未发生事件也未被删失）。

KM方法能够处理右删失数据，对于删失的个体只在其删失时间之前作为风险集的一部分，在删失之后则不再计入为$n_i$。这样的处理依赖于KM对删失的假设，即非信息性删失，假设删失与生存时间是独立的。

• 对于$\hat{S}(t)$的估计的精确程度可以估计置信区间，使用Greenwood公式：

$\mathrm{Var}[\hat{S}(t)]=\hat{S}(t)^2\sum_{t_i\leq t}\frac{d_i}{n_i(n_i-d_i)}$

• KM本身用于估计，若要比较不同组的生存曲线是否显著不同，可以使用Log-Rank检验，也可以使用Wilcoxon检验等加权方法。
• 缺点：不能调整协变量、无法识别复杂的生存模式、难以处理删失比例高的情况。

5 模型统计指标

5.1 ROC曲线及曲线下的面积AUC

ROC曲线（接收者操作特征曲线）通过绘制假阳性率（FPR）和真阳性率（TPR）的变化来描述分类模型的性能。其基本思想是随着阈值的变化，模型的分类的结果也会发生变化。ROC曲线则描绘了这种变化对分类模型的性能的影响。

• X轴：假阳性率（False Positive Rate, FPR）：与特异度（Specificity）互补（和为1），被错误地预测为阳性的负类样本的比例：

  $FPR=\frac{False \ Positives}{False \ Positives + True \ Negetives}$

  分母则为所有的真实为负类样本的个数

• Y轴：真阳性率（True Positive Rate, TPR）：也称为召回率（Recall）或灵敏度（Sensitivity），表示正确地预测为阳性的正类样本的比例：

  $TPR=\frac{True \ Positives}{True \ Positives + False \ Negetives}$

  分母则为所有的真实为阳性样本的个数

绘制过程：分类模型会输出每一个样本属于某个类别的概率，对于预测给出的概率值，可以设定不同的阈值决定哪个概率值大于阈值时为阳性。随着阈值的改变，FPR和TPR会发生相应的变化。对于每一个可能的阈值计算对应的FPR和TPR在二维坐标系中标出。

ROC曲线：曲线左下角(1,0)表示分类器完全不能区分正负样本，FPR为1，TPR为0；曲线右上角(0,1)意味着分类器完全区分正负类，所有样本都被正确分类。从左至右，代表FPR越来越小，即阈值越来越严格，越来越少的样本被错误预测为阳性，沿着ROC曲线，越来越多的样本被正确分类为真阳性。

5.1.1 ROC-AUC

ROC曲线下的面积被称为AUC，表示模型区分正负样本的能力，从ROC曲线的定义出发可以很好的理解，ROC曲线下则代表了在所有阈值下TPR大于FPR的概率，AUC越大，模型的分类性能越好。一般与AUC=0.5作比较，AUC=0.5代表随机猜测。事实上，没有模型会出现AUC小于0.5的情况，这说明模型的分类完全反向，此时应当是模型的设置出现了问题。

5.1.2 Delong’s Test

用于比较两个模型的AUC的非参数统计检验方法，以判断是否存在显著差异。

原理：Delong检验基于非参数方法计算AUC的方差，并使用AUC差异的方差来评估两种模型之间的差异是否显著。Delong检验首先估计两个模型的AUC，然后计算它们的差值，再对这个差值进行统计显著性检验，以确定差值是否为零。

步骤：

• 假设：原假设$H_0$：两个模型的AUC无显著差异；备择假设$H_1$：两个模型的AUC有显著差异。

• 步骤：

  计算每个模型的AUC：

  • 使用ROC曲线计算两个模型在同一数据集上的AUC值，分别记为 $AUC_1, AUC_2$。

  计算每个AUC的得分：

  • 对于每个样本，计算模型预测的分数（例如预测概率或得分）。通过所有正负样本对的比较，计算每个样本对AUC贡献的“得分”。

  计算AUC差异的方差：

  • 使用U统计量的方法，将每个模型的AUC得分映射到样本上。
  • 对于两个模型，计算每个样本的得分差异，并基于这些得分差异估计出AUC差异的方差。

  计算检验统计量：

  • 计算两个AUC差异的标准误：

    $SE(\Delta AUC)=\sqrt{\mathrm{Var}(AUC_1)+\mathrm{Var}(AUC_2)-2\cdot\mathrm{Cov}(AUC_1,AUC_2)}$。

  计算z值：

  • 使用以下公式计算z统计量：$z=\frac{\Delta AUC}{SE(\Delta AUC)}$。

5.2 P-R曲线

相比于ROC曲线，PR曲线（Precision-Recall Curve）更适合用于在数据不平衡的情况下衡量分类器的性能。ROC曲线更用于评估分类器的整体性能，PR曲线则更加聚焦于正类分类的情况。

• Why？在不平衡数据集中，（一般是负类样本要远多于正类样本）FPR的值受负类样本数量的影响较大。例如极端情况下，所有样本为负类，FPR仍然会很低，这种情况下ROC曲线会高估模型的性能。

• X轴：精确率（Precision）：所有样本中被预测为阳性的样本中，本身就是阳性的比例，

  $Precision=\frac{True \ Positive}{True \ Positive + False \ Positive}$

• Y轴：真阳性率

PR曲线的形状可以反映分类器在不同阈值下的 精确率 和 召回率 的权衡：

• 精确率和召回率的权衡：通常精确率和召回率是存在权衡，提高一个指标（如召回率）可能会导致另一个指标（如精确率）的下降。PR曲线可以帮助我们直观地看到这一点。
• 理想的PR曲线：PR曲线理想情况下应该是 尽可能接近右上角。在这种情况下，精确率 和 召回率 都比较高，意味着分类器不仅能准确地预测阳性样本，而且能找到大多数阳性样本。
• 曲线靠近对角线：如果PR曲线靠近对角线（即斜对角线），意味着模型的表现较差，精确率和召回率较低。

在二分类问题中，随机分类器的精确率等于样本正类比例，随机分类器在预测正负类时没有利用任何特征信息，仅依赖样本的总体分布。因此，在随机分类的情况下，预测为正的样本中实际为正类的比例会接近样本的正类比例。在P-R曲线中可以对应标记出。

5.3 曲线绘制逻辑

不管是ROC曲线还是P-R曲线时，其绘制的逻辑都是：当判断为阳性的阈值变化时混淆矩阵的值都会发生相应的变化。

混淆矩阵：

当模型将样本判断为阳性的阈值降低时，模型判断为阳性的标准越来越宽松，此时越来越多的样本会被轻松的判断为阳性，但很显然这并不是一个线性的过程（这也是ROC和P-R曲线会出现平台的原因）。当预测为阳性的样本的数量升高时，TP和FP的数量也会升高，相应的FN和TN的数目会降低。

至此，我们可以得到混淆矩阵的变化方向，再可以落实到Precision, Recall, Sensitivity以及Specificity的变化情况。

• $Recall=Sensitivity=TRP=\frac{True \ Positives}{True \ Positives + False \ Negetives}$,可知其分母为所有的真实为正类的样本，所以阈值降低时，TPR是会升高的。

• $FPR=\frac{False \ Positives}{False \ Positives + True \ Negetives},Specificity=1-FPR$，当阈值降低时FP升高，TN降低，此时FPR整体变大，因此Specificity整体会降低。

• $Precision=\frac{True \ Positive}{True \ Positive + False \ Positive}$，当阈值降低时TP,FP都升高，此时就涉及到了TP和FP的变化速度的问题，但是我们能够得到一个整体的变化趋势：当阈值极端为1时，模型会将所有的样本划分为负类，此时FP和TP都为0，Precision趋于1；当阈值极端为0时，模型会将所有样本分为正类，此时会有大量的负样本涌入，此时FP较大，TP在阈值的某一个点后不再会发生变化，因此Precision也会靠近于0（不一定趋于0）。（上述思考是在阳性样本少于阴性样本的情况下讨论的，这也是符合我们日常需求的，如果真实情况相反，只需要相反标记即可）

5.4 IDI（Integrated Discrimination Improvement）综合判别能力改进

定义：IDI用于评估新模型（如添加新的预测因素后）在区分事件（如疾病发生）和非事件（如疾病未发生）个体上的综合判别能力是否改进。IDI计算了新模型与原模型之间预测概率的平均差异。

计算原理：在新模型中，事件（如患病）和非事件（如未患病）个体的预测概率应该更为分离，即事件个体的预测概率应该更高，非事件个体的预测概率应该更低。IDI通过计算事件和非事件个体预测概率的变化来量化这一改善。

事件组的预测概率改进：$\Delta p_{\mathrm{event}}=\bar{p}_{\mathrm{new,~event}}-\bar{p}_{\mathrm{old,~event}}$

非事件组的预测概率改进：$\Delta p_{\mathrm{nonevent}}=\bar{p}_{\mathrm{new,~nonevent}}-\bar{p}_{\mathrm{old,~nonevent}}$

$\mathrm{IDI}=\Delta p_\text{event}-\Delta p_\text{nonevent}$

IDI的值越大，说明新模型相比旧模型，区分事件和非事件个体的能力越强。如果IDI接近0，表示新模型在判别能力上并无明显改进。

应用：在评估新的生物标志物或风险评分（如MRS，PRS）是否提升模型的预测性能时，IDI可以作为判断模型改进的一项指标。

5.5 NRI（Net Reclassification Improvement）净重新分类改进

定义：NRI是一种评估新模型相对于原模型在重新分类方面的改进程度的指标。它通过考察新模型是否能够将事件和非事件个体更准确地重新分类到风险类别中来量化改进。

计算原理：NRI将个体分为“事件组”和“非事件组”，并在新模型和原模型之间进行分类比较。根据新模型是否将事件个体划分到更高风险类别，或将非事件个体划分到更低风险类别，NRI分别计算事件和非事件的分类改进率，然后取两者的和，得到净重新分类改善值。

计算步骤：

1. 定义风险类别：

   • 根据实际需要，将预测概率划分为低风险、中风险和高风险等不同类别。例如：低风险（<10%），中风险（10%-20%），高风险（>20%）。

2. 计算事件组的重新分类改进：

   • 对于实际发生事件的个体，如果新模型将个体重新分类到更高风险类别，则视为分类改进。

   • 如果新模型将个体重新分类到更低风险类别，则视为分类退步。

   • 事件组NRI为重新分类到更高类别的人数比例减去重新分类到更低类别的人数比例：

     $NRI_{event}=\frac{重新分类到更高类别的人数-重新分类到更低类别的人数}{事件组人数}$

3. 计算非事件组的重新分类改进：

   • 对于实际未发生事件的个体，如果新模型将个体重新分类到更低风险类别，则视为分类改进。

   • 如果新模型将个体重新分类到更高风险类别，则视为分类退步。

   • 非事件组NRI为重新分类到更低类别的人数比例减去重新分类到更高类别的人数比例：

     $NRI_{nonevent}=\frac{重新分类到更低类别的人数-重新分类到更高类别的人数}{事件组总人数}$

4. 计算总NRI：

   • NRI为事件组和非事件组的分类改进和： $\mathrm{NRI=NRI_{event}+NRI_{nonevent}}$

如果NRI > 0，则说明新模型在重新分类方面相较于旧模型有所改善，特别是在将事件个体划归更高风险类别或将非事件个体划归更低风险类别上表现更好。如果NRI接近0或为负，表明新模型在重新分类方面没有优势或有退步。

应用：NRI常用于临床预测模型的评估。例如，在ACS风险预测中，通过NRI可以衡量MRS模型是否更好地将高风险和低风险个体重新归类。

5.6 C-Index（Concordance Index）（分类情况）

又称Harrell’s C-index，C statistics。

C统计量用于评价模型的预测准确性，针对二分类逻辑回归模型，C-Index退化为“正样本被预测为阳性的概率”大于“负样本预测为阳性的概率”的可能性。经过证明，二分类逻辑回归的C-Index等价于AUC的值。

可以简单的理解为：C-Index是AUC的一般情况，AUC是C-Index在二分类模型的特殊情况。

C指数的一般定义：把所有样本中的对象随机两两组队，共有$C_n^2$种，计算预测结果一致的组队的占比。C-Index本质上估计了预测结果与实际观察到的结果相一致的概率。在分类问题种则描述了正类样本的预测概率大于负类样本的预测概率的可能性。

C-Index的一般形式一般用于生存分析中，在此不展开。

5.7 DCA（Decision Curve Analysis）曲线

决策曲线分析是一种用于评估预测模型临床实用性的方法，特别适用于二分类问题。DCA的核心思想是基于净收益（Net Benefit，NB）来衡量不同模型的优劣，并且通过决策曲线展示模型在不同风险阈值下的临床价值。

相比于传统的评估方法（ROC曲线、AUC、校准曲线）：

• AUC仅仅衡量模型的区分能力，但是无法反应模型的临床价值；
• 校准曲线评估了预测概率的准确性，但是没有考虑实际决策中的收益；
• 传统决策分析方法需要明确的成本和收益，而DCA基于相对收益的概念，不需要明确的成本数据。

原理：$NB=\frac{TP}{N}-\frac{FP}{N}\times\frac{p_t}{1-p_t}$

其中第一项表示正确筛查出的患者的比例，即模型带来的真实收益；第二项代表错误筛查出的或者带来的相对损失，其中$p_i$表示决策阈值（医生或决策者愿意采取干预措施的概率阈值），赋予的系数表示错判一个阴性个体为阳性个体相对于正确判断一个阳性的收益，因为假阳性会导致不必要的干预。当NB大于0时说明该模型在该决策阈值下有实际价值。

绘制：在绘制时除了绘制DCA曲线外，还会绘制无模型基准线（假设不筛查任何人，NB=0）以及全员治疗基准线（假设筛查所有人）。

6 因果推断方法

6.1 MR分析

孟德尔随机化（Mendelian Randomization, MR）是一种基于遗传学的因果推断方法，用于评估暴露因素与结局之间的因果关系。该方法利用了基因变异的随机分布，这类似于随机对照试验中的随机化。具体来说，孟德尔随机化将遗传变异视为工具变量，用以检测暴露变量（如某种生活方式、代谢指标）与感兴趣的结局变量（如疾病）的因果关系。

一、方法原理

孟德尔随机化基于孟德尔遗传定律，即父母的基因随机地传递给后代，这在种群水平上相当于一个自然的随机对照试验。遗传变异在个体之间的分布是随机的，因此这些基因变异可以作为工具变量（Instrumental Variable, IV），用来控制潜在的混杂因素，从而推断暴露与结局之间的因果关系。

一个典型的MR研究包括以下三个主要假设：

1. 相关性假设（Relevance Assumption）：工具变量与暴露变量密切相关，基因变异显著影响暴露。
2. 独立性假设（Independence Assumption）：工具变量与任何潜在的混杂因素无关，换句话说，基因变异与暴露之外的其它因素无相关性。
3. 排除限制假设（Exclusion Restriction Assumption）：工具变量通过暴露影响结局，不能直接影响结局，即基因变异对结局的影响完全通过暴露介导。

二、步骤

选择工具变量：基因变异（通常是单核苷酸多态性，SNPs）被选作工具变量。要选择合适的SNP，这些基因位点与感兴趣的暴露变量（如胆固醇水平、体重指数）有显著关联。

确定暴露与结局的关联：收集关于暴露和结局的遗传相关数据。这些数据通常来源于全基因组关联研究（GWAS），通过对群体中SNP与暴露变量、结局变量的关联进行分析。

进行因果推断：使用工具变量估计暴露与结局之间的因果效应。MR分析常用的统计方法包括：

• 双阶段最小二乘法（2SLS）：首先估计工具变量与暴露的关系，再估计暴露与结局之间的因果关系。
• Wald比值法（Wald Ratio）：对于单个SNP的分析，该方法通过暴露对SNP的效应与结局对SNP的效应比值得出因果估计。
• 加权中值估计：当多个SNP用作工具变量时，这种方法通过对各SNP的估计值进行加权平均以提供更稳健的因果效应估计。

敏感性分析：由于基因变异可能并非完全符合MR假设，敏感性分析如MR-Egger回归等方法用于检测并纠正潜在的基因水平效应（如基因多效性，pleiotropy）。

三、优点

控制混杂因素：由于基因的随机分配，孟德尔随机化在控制混杂因素方面具有天然优势，尤其是在处理难以控制的环境因素时。

减少反向因果关系：因为遗传变异在出生时已确定，因此它们不受结局的反向影响，降低了反向因果问题。

无偏效应估计：在严格满足工具变量假设的前提下，MR提供无偏的因果效应估计。

四、缺点

弱工具变量问题：如果工具变量与暴露变量的关联较弱，MR分析的统计效能可能不足，导致因果估计不可靠。

基因多效性：如果某个SNP同时影响暴露变量和结局变量，可能违反独立性假设，导致因果推断偏差。这种情况可以通过敏感性分析部分校正。

外部效应：工具变量的效应可能受到外部因素（如环境因素）的影响，导致因果关系估计不准确。

五、R包

install.packages("remotes")
remotes::install_github("MRCIEU/TwoSampleMR")

使用TwoSampleMR包进行MR分析，需要准备两个主要的数据集：暴露数据（exposure data）和结局数据（outcome data）：

暴露数据包含：SNP的rsid号，effect_allele效应等位基因，other_allele非效应等位基因，beta（SNP与暴露的关联效应值），se（beta的标准差），pval（SNP与暴露关联的p值），samplesize（非必要）

结局数据包含：SNP的rsid号，effect_allele效应等位基因，other_allele非效应等位基因，beta（SNP与结局的关联效应值），se（beta的标准差），pval（SNP与结局关联的p值），samplesize（非必要）

MR的基本原理是利用遗传变异（SNPs）作为工具变量，这些SNPs应该是随机分配的（孟德尔遗传定律），通过这种随机化来模拟平常不易开展的随机对照试验（RCT）。

其中包括三个核心假设：相关性假设：SNPs与暴露因素显著相关（因此要提供SNP-exposure的beta和p值）；独立性假设：SNPs不通过其他途径影响结局（因此要选择独立的SNPs）；排他性假设：SNPs只通过暴露因素影响结局。

具体路径：SNPs 👉 暴露因素 👉 结局：通过比较SNPs对结局的效应和SNPs对暴露因素的效应来估计暴露因素对结局的因果效应。此时：effect_allele和other_allele是为了确保效应方向的一致性。

6.2 基于逆方差加权法IVW的MR因果推断方法

在 MR 分析中，逆方差加权法（Inverse-Variance Weighted, IVW） 是一种常用的方法，用来整合所有 SNP 的效应比值，进而估计暴露（如蛋白质水平）对结局（如疾病）的总体因果效应。IVW 的底层逻辑可以分为以下几个关键步骤：

1. 计算每个 SNP 的效应比值

假设有多个 SNP，每个 SNP 的 暴露效应（对蛋白质水平的影响） 和 结局效应（对疾病的影响） 分别用以下符号表示：

对暴露的效应：$\beta_{\text{exposure,}i}$

对结局的效应：$\beta_{\text{outcome,}i}$

对结局效应的标准误：$SE_{\text{outcome,}i}$

对于每个SNP，其效应比值可以计算为：$\mathrm{Ratio}_i=\frac{\beta_{\text{outcome,}i}}{\beta_{\text{exposure,}i}}$，表示 SNP 通过暴露对结局的影响程度，或暴露每增加一个单位对结局产生的效应。

2. 加权因果效应估计

由于每个 SNP 的效应比值存在不确定性（即标准误），希望对这些效应比值进行加权平均。IVW 的核心思路是对每个效应比值使用其逆方差作为权重，因为方差越小的 SNP（即不确定性越低）应该在总体估计中占更大比重。对于第$i$个SNP，权重$w_i$由该SNP在结局效应上的逆方差决定：$w_i=\frac1{\text{SE}_{\text{outcome,}i}^2}$

然后，使用这些权重对每个SNP的效应比值$Ratio_i$求加权平均，从而得到暴露对结局的总体因果效应估计值：$\hat{\beta}_{\mathbf{IVW}}=\frac{\sum_iw_i\cdot\mathrm{Ratio}_i}{\sum_iw_i}$。

3. 估计总体因果效应的标准误

加权平均的总体因果效应估计值$\hat{\beta}_{\mathbf{IVW}}$本身存在不确定性：$\mathrm{SE}_{\mathrm{IVW}}=\sqrt{\frac1{\sum_iw_i}}$，反应了总体因果效应估计值的不确定性，当权重越大（即各个SNP的方差越小时）标准误更小，结果更稳健。$\mathrm{Var}(\hat{\beta}_\mathrm{IVW})=\mathrm{Var}\left(\frac{\sum_iw_i\cdot\mathrm{Ratio}_i}{\sum_iw_i}\right)=\frac{\sum_iw_i^2\cdot\mathrm{Var}(\mathrm{Ratio}_i)}{\left(\sum_iw_i\right)^2}=\frac{\sum_iw_i^2\cdot\frac1{w_i}}{\left(\sum_iw_i\right)^2}=\frac1{\sum_iw_i}$ （独立性假设）

4. 计算p值，检验因果效应的显著性

通过因果效应估计值以及标准误，能够计算其p值，检验暴露对结局的因果效应是否显著：$\text{Z-score}=\frac{\hat{\beta}_\mathrm{IVW}}{\mathrm{SE}_\mathrm{IVW}}$。

6.3 基于Weighted Median的MR因果推断方法

加权中位数估计法是一种基于中位数的MR因果推断方法，通过对多个遗传工具变量的效应进行加权中位数计算，从而估计暴露变量与结局变量之间的因果效应。

核心思想：只要超过50%的遗传工具变量有效，就可以得到稳健的因果效应估计。

步骤：1.将每个工具变量对暴露变量的效应大小排序；2.对每个工具变量的估计因果效应赋予权重，权重通常是工具变量效应的标准误的倒数。3.通过加权中位数值，得到暴露变量对结局变量的因果效应估计。

优势：1.稳健性：在存在一定比例的无效工具变量时仍然能够提供稳健的因果效应估计；2.减少多效性偏倚；3.对异常值的敏感性较低。

局限性：相比较于IVW，需要较大的样本量；无法完全消除多效性的偏倚；方法依赖于有效工具变量超过50%的假设。

6.4 基于MR Egger的MR因果推断方法

MR-Egger是一种用于评估和校正水平多效性（pleiotropy）的工具，水平多效性指的是工具变量通过多种不同的生物途径影响因变量，可能会引入偏倚影响因果推断的准确性。

MR-Egger的截距项：截距项如果显著偏离0，则表明存在着系统的水平多效性，代表了工具变量不仅通过暴露变量影响结局变量，还通过其他途径产生影响。

MR-Egger的斜率项：用于估计暴露变量对结局变量的因果效应，表示在调整了水平多效性之后的暴露变量与结局变量之间的关联。

作用：MR-Egger在因果估计中可以调整多效性带来的影响，该方法假设遗传工具变量（比如SNP位点）与多效性效应的关联是随机的（无效工具强度独立于直接效应InSIDE），在满足该假设下，MR-Egger能够给出相对稳健的因果效应的估计。

• InSIDE假设：遗传工具变量的效应强度（工具与暴露变量的关联强度）与工具变量通过其他途径直接影响结局变量的效应独立。

局限性：相比于IVW需要更大的样本量；过度校正多效性带来的影响；InSIDE假设的限制

7 共定位分析

MR分析可能存在多效性，同一个位点可能通过不同机制影响暴露和结局，也有可能是LD的SNP而非是真正的因果变异。共定位分析用于检验暴露和结局的遗传关联是否来自相同的因果变异，帮助确认是否是同一个因果变异驱动了两种表型，从而降低假阳性的风险。

共定位分析用于检验暴露和结局的GWAS信号是否共享相同的因果变异，使用贝叶斯方法估计：

H0：该区域没有任何变异与两个形状相关的关联信号

H1：仅存在与暴露相关的变异

H2：仅存在与结局相关的变异

H3：暴露和结局各自有独立的因果变异

H4：暴露和结局共享同一个因果变异

PP4>0.8，强烈支持共定位；PP3>0.8，支持不同因果变异。

因此：P4>0.8且MR显著：最强的因果关系证据

同样的，共定位分析需要用到两个GWAS的summary statistic的数据。每个数据集需要包括：SNP的rsid号，染色体位置，效应等位基因、其他等位基因、效应大小、标准误、p值、等位基因频率、样本量等。

8 特征工程-变量选择

最基础的有向前选择法，向后选择法，双向选择法，以及基于Lasso的自动选择变量的方法。

前置概念：

• 包装器（Wrapper）：代指特征选择中的一种方法论类别，其工作原理是将特征选择过程”包装“在学习算法周围，使用目标学习算法的性能作为特征子集评价的标准，直接针对特定的学习算法优化特征选择。其基本步骤是：生成候选特征子集，使用学习算法训练模型，评估模型性能，根据性能指标决定是否保留该特征子集。例如前向选择，Bruta算法

• 过滤器（Filter）：独立于学习算法，使用统计指标来选择特征；例如相关系数的方法

• 嵌入式（Embedded）：特征选择在模型训练过程中完成。例如Lasso回归的方法

8.1 基于随机森林的Boruta算法

• 阴影特征（Shadow Features）：通过随机打乱原始特征的值所人工创建的特征，保持了与原始特征相同的统计分布，但是破坏了与目标变量的关系。创建方法：通过对原始特征的值进行随机排列，为每个原始特征创建对应的阴影版本，保持数据分布的特征，只改变顺序。如果原始特征确实重要，其重要性分数应显著高于阴影特征。提供了一个数据驱动的特征选择标准

一种基于随机森林的分类包装器算法（包装器算法指的是排除不太重要的特征）。

其目的是选择出所有与响应变量相关的特征，而不是选择使损失函数最小的特征集合。其核心思想是：如果一个特征是相关的重要变量的话，那么它的原始特征重要性要优于其随机版本（阴影特征）。

算法原理：

1. 特征重要性评估基础

• 利用随机森林的平均精确度下降(Mean Decrease Accuracy)或平均Gini系数下降(Mean Decrease Gini)作为特征重要性度量
• 通过随机排列特征值创建阴影特征，以模拟随机性背景

1. 具体执行步骤：

• 对原始数据集中的每个特征，创建对应的阴影特征(通过随机打乱原特征值)
• 使用扩展后的数据集(包含原始特征和阴影特征)训练随机森林模型
• 计算所有特征(原始和阴影)的重要性分数
• 找出阴影特征中重要性最高的特征(maxShadow)，将其作为阈值
• 比较原始特征与maxShadow：
  • 显著高于maxShadow的特征被标记为”重要”
  • 显著低于maxShadow的特征被标记为”不重要”
  • 其他特征保持”暂定”状态
• 移除被标记为”不重要”的特征
• 重复以上步骤，直到：
  • 所有特征都被标记为”重要”或”不重要”
  • 达到预设的最大迭代次数

1. 统计显著性判断：

• 使用双边Z检验比较特征重要性与maxShadow
• 通常采用95%或99%的置信水平
• 多次迭代中累积的决策用于最终确定特征的重要性

8.2 mRMR（最小冗余最大相关）方法

mRMR, Minimum Redundancy Maximum Relevance：核心思想是在选取特征时，既要保证选出的特征与目标变量高度相关，又要保证特征之间的冗余度最小。

假设有一个数据集$\mathcal{D}=\{X_1,X_2,...,X_n,Y\}$，其中X代表特征，Y代表目标变量，mRMR方法通过互信息（Mutual Information, MI）衡量变量间的相关性和冗余性。

• 最大相关性：相关性用互信息I(X; Y)来衡量：$\max S=\frac{1}{|S|}\sum_{X_i\in S}I(X_i;Y)$，其中S表示已选的特征子集；
• 最小冗余性：$\min R=\frac{1}{|S|^2}\sum_{X_i,X_j\in S}I(X_i;X_j)$，衡量了特征之间的互信息。
• 目标函数：$\max(S-R)=\max\left(\sum_{X_i\in S}I(X_i;Y)-\sum_{X_i,X_j\in S}I(X_i;X_j)\right)$，该形式也被称为互信息差值（MID）。MID在特征维数较高时比较稳定，相对应的有互信息商（MIQ）：$\max\frac{S}{R}=\max\frac{\sum_{X_i\in S}I(X_i;Y)}{\sum_{X_i,X_j\in S}I(X_i;X_j)}$，该方案在高相关特征数据集上效果更好。

算法实现策略，贪心搜索算法，通过迭代方法选择特征：

1、计算与因变量最相关的特征作为第一个特征；

2、对于剩余的变量，计算特征与目标变量的互信息，计算与去除改变了之后的剩余其他特征的平均互信息，计算mRMR评分。

3、选择具有最高mRMR评分的特征加入特征子集。

4、重复步骤2和3，直至满足终止特征（比如特征贡献度下降）。

8.3 遗传算法（Genetic Algorithm, GA）

遗传算法是一类基于自然选择和遗传机制的搜索和优化算法。主要模拟生物进化过程，包括选择、交叉、突变等。在特征选择中，GA算法将特征子集表示为染色体，通过遗传操作逐步优化出最优的特征组合。

步骤：

1. 编码（Chromosome Representation）：在GA中，每个个体代表一个特征子集，通常采用二进制编码，长度代表特征总数，每个位（基因）表示某个特征是否被选中。例如，对如一个3个特征的特征集，一个染色体可能是{1, 0, 1}表示选择了第1、3个特征。
2. 适应度函数（Fitness Function）：衡量染色体的质量，常见的适应度函数包括：分类准确率、AUC、均方误差等。为了防止选择过多的特征，对特征数目采取正则化的措施。

3. 选择：选择操作决定哪些个体能够进入下一代，通常采用：

   • 轮盘赌选择：根据适应度值的大小以概率方式选择个体；
   • 锦标赛选择：随机挑选部分个体，选择适应度最高的；
   • 精英保留：直接保留适应度最高的个体。

4. 交叉（Crossover）：交换染色体部分基因生成新的个体，通常采用：

   • 单点交叉：在随机位置切割并交换部分基因；
   • 两点交叉：选择两个点进行基因交换；
   • 均匀交叉：每个基因位点以某个概率进行交换。
5. 变异（Mutation）：对某些位进行随机翻转，增加种群多样性，防止陷入局部最优解，一般采用低概率进行变异。
6. 终止条件：迭代达到最大代数；适应度值收敛。

缺点：计算复杂度高，计算成本高，难以调参，不保证寻找到全局最优解。

优点：比贪心算法更容易跳出局部最优解；适用于高维数据；不依赖任何分类器。

8.4 模拟退火算法（Simulated Annealing，SA）

模拟退火算法（Simulated Annealing, SA）是一种基于物理退火过程的全局优化算法。其核心思想源自金属材料退火的物理过程：金属加热到高温后，逐步冷却，使其内部粒子逐步趋于低能量稳定状态。据此，在优化问题中，SA 通过引入“温度”参数，在搜索初期允许接受次优解以避免陷入局部最优解，并在温度逐渐降低的过程中收敛到最优解。

步骤：

• 状态表示，一个状态对应一个特征子集，常见的表示方法有：二进制编码（参考遗传算法）；特征索引列表：直接存储当前选择的特征索引集合。
• 目标函数：用于评估当前特征子集的质量，通常基于模型的预测性能和特征子集的复杂度。
• 初始状态：随机选择一组特征；或者全选；或者全不选；基于某些启发式方法得到的初始特征子集。
• 邻域搜索策略：随机添加或删除一个特征；随机替换一个特征；交换两个特征。
• 迁移概率（Metropolis准则）：如果新特征子集的目标函数值优于当前特征子集，则接受；否则以一定概率P接受较差的解：$P=\exp\left(\frac{f(S^{\prime})-f(S)}{T}\right)$，T表示当前温度，会随着迭代逐渐降低。这一策略允许在初始阶段接受较差解，以跳出局部最优解。

• 降温策略：温度T的更新方式决定了算法的搜索效率，常见的降温方案包括：

  • 几何降温：$T_{k+1}=\alpha T_k$，$\alpha$通常取0.8~0.99；
  • 对数降温：$T_k=\frac{T_0}{1+\beta\log(1+k)}$
  • 线性降温：$T_k=T_0-k\Delta T$

• 终止条件：温度降到某个阈值；目标函数多次迭代未改善；达到最大迭代次数。

优点：跳出局部最优解；适用于非凸优化问题；可以灵活调整目标函数；

缺点：成本高；调参复杂；结果稳定性依赖于初始温度和冷却策略。

8.5 RFE递归特征消除

Recursive Feature Elimination, RFE主要用于选择对模型最重要的特征。其核心思想是训练模型（评估所有特征的重要性），消除最不重要的特征，重复步骤。RFE主要依赖于一个基础学习器来评估每个特征的重要性，并且采用递归的方式减少特征。

QTL

QTL（Quantitative Trait Locus）分析，即“数量性状基因座分析”，是一种用于定位与数量性状相关的遗传区域（基因座）的遗传学方法。数量性状（如身高、血压、产量、体重等）往往由多个基因和环境因素共同作用，因此无法简单地归结于某一单一基因。QTL分析通过遗传标记与性状数据的关联分析，来寻找影响这些复杂性状的基因区域，是作物育种、动物育种以及人类疾病遗传学研究的重要工具。

pQTL

pQTL（Protein Quantitative Trait Loci） 是一种用于探究基因变异对蛋白质水平的影响的研究方法，具体来说，它是通过统计学分析来识别和定位与特定蛋白质表达水平相关的基因组位点（即 SNPs）。它类似于 GWAS（基因组全关联研究），但重点不在于性状（如疾病），而是蛋白质的浓度或表达水平。pQTL 的研究可以帮助我们理解哪些基因变异会影响特定蛋白质的表达，从而揭示蛋白质调控背后的遗传基础。

理解：可以把 pQTL 想象成在基因组中寻找影响“蛋白质开关”的基因位点。某些基因变异（SNPs）就像是“开关”，它们的不同状态（等位基因）会影响蛋白质的表达量。例如，一个 SNP 的特定等位基因可能导致特定蛋白质的表达量升高或降低，从而影响生物体的功能。

核心目的：pQTL最后生成的文件是为了找到与某蛋白质表达水平显著相关的遗传位点。

目的：生成 pQTL 结果文件的主要目标是鉴定影响蛋白质水平的遗传因素。这项研究的核心在于：

• 揭示基因调控网络：找出影响蛋白质表达的关键遗传变异，理解基因如何通过蛋白质影响生物过程。
• 推动药物研发：pQTL 研究可以帮助识别与疾病相关的蛋白质，并进一步指向潜在的药物靶点。例如，如果某个 SNP 显著影响与心脏病相关的蛋白质水平，这个 SNP 及其关联的基因可能是潜在的治疗靶点。
• 疾病因果研究：通过 pQTL，研究者可以在孟德尔随机化分析中评估这些蛋白质是否具有因果效应。鉴于基因是天然的“随机对照”变量，pQTL 提供了一个接近因果推断的机会。

步骤

1. 样本和数据准备

pQTL 分析的第一步是选择样本和收集遗传信息及蛋白质表达信息。通常需要以下数据：

• 蛋白质表达水平：在血液、血浆或组织样本中测量蛋白质浓度或表达水平。蛋白质浓度可以通过高通量技术，如 Olink、SOMAscan 或 质谱分析 等来获得。
• 基因型数据：通过基因分型芯片或全基因组测序获取样本的遗传变异信息（如 SNP 数据）。

样本数量

• 进行 pQTL 分析时，需要足够大的样本量以确保统计效力。通常 pQTL 分析要求几千到数万个样本，这样才能更有可能检测到弱效的基因变异对蛋白质表达的影响。

2. 蛋白质测量和遗传变异数据收集

pQTL 分析的核心在于蛋白质水平和遗传变异的关联。下面是这两类数据的采集方法：

• 蛋白质浓度测量：使用高通量的蛋白质分析平台测量大量蛋白质的浓度。例如：
  • Olink：通过抗体和 PCR 技术的结合来检测蛋白质浓度，灵敏度高，适合血浆和血清中的蛋白质。
  • SOMAscan：使用适配子技术检测蛋白质，覆盖数千种蛋白。
  • 质谱分析：适合检测多种类型的蛋白质，尤其是针对特定组织样本。
• 基因分型：通常使用基因分型芯片或者全基因组测序技术获得基因型数据，捕捉样本中可能存在的遗传变异信息。

数据对齐：

蛋白质测量和遗传变异数据需要在同一批样本中测得，以便可以将蛋白质表达数据与基因型数据相匹配。

3. 数据预处理

在数据分析之前，需要对蛋白质和基因型数据进行严格的预处理：

• 蛋白质数据预处理：
  • 去除异常值：检测并去除蛋白质浓度中的异常值，避免极端值对分析的影响。
  • 标准化：将蛋白质浓度数据进行标准化，通常使用 log 转换或 z-score 标准化，使得蛋白质水平具有更接近的分布，减少测量误差的影响。
• 基因型数据预处理：
  • 质控：检查基因型数据的质量，例如 SNP 的最小等位基因频率（MAF）过滤、遗传距离过滤等，以确保数据的准确性。
  • 填补缺失数据：使用统计方法（如最大似然法）填补基因型中的缺失值，保证每个样本的基因型数据完整。

经过预处理的蛋白质数据和基因型数据才能进入下一步的关联分析。

4. pQTL 关联分析

在关联分析中，pQTL 研究的核心是检测遗传变异（SNP）与蛋白质表达水平之间的关联。常用的分析步骤如下：

• 单位点关联分析：
  • 对于每个蛋白质，针对每个 SNP 进行回归分析。例如，假设蛋白质浓度为因变量、SNP 基因型为自变量，使用加性模型来评估 SNP 与蛋白质表达的关系。
  • 回归模型中通常会调整潜在混杂变量（如年龄、性别、测量批次等）以减少偏差。
  • 该分析会生成一个回归系数（效应量）、标准误差和显著性 p 值。
• 顺式 pQTL 与反式 pQTL：
  • 顺式 pQTL（cis-pQTL）：SNP 位于编码该蛋白质的基因附近（通常定义为基因附近 1 Mb 范围内）。顺式 pQTL 可能直接影响基因表达，具有更直接的生物学意义。
  • 反式 pQTL（trans-pQTL）：SNP 位于基因组的其他位置，可能通过间接机制（如信号通路）影响蛋白质水平。反式 pQTL 通常比顺式 pQTL 更难解释，但可能揭示远程调控机制。

5. 多重检验和显著性阈值

在进行数百万个 SNP 和数千种蛋白质的关联分析时，多重检验是一个重要的问题。

• 多重检验校正：由于分析了大量的 SNP，需要使用多重检验校正（如 Bonferroni 校正或 FDR 校正）来减少假阳性结果。
• 显著性阈值：通常设定 p 值阈值为 5e-8（GWAS 标准）或通过校正后的显著性阈值来判断 SNP 是否显著关联。
• 统计功效：较大的样本量可以提高检测到显著性关联的功效，尤其是对效应较小的 SNP。

6. 结果解读和生物学验证

在获得 pQTL 的显著结果后，下一步是对这些结果进行生物学解释和验证。

• 注释和功能分析：将显著 SNP 注释到基因组上，查看它们是否位于基因内部、附近调控区或远程调控区。
• 功能验证实验：使用分子生物学实验（如 CRISPR 基因编辑、RNAi）验证 SNP 对蛋白质表达的影响。
• 网络分析：对显著 SNP 关联的蛋白质进行通路和网络分析，识别潜在的调控网络和关键节点。

7. 结果应用

pQTL 分析的结果应用非常广泛，包括以下几方面：

• 孟德尔随机化：pQTL 结果可以用于孟德尔随机化分析，将显著的 pQTL SNP 用作工具变量，评估蛋白质是否对特定疾病有因果影响。
• 药物靶点识别：与蛋白质相关的显著 pQTL SNP 可以指向潜在的药物靶点，特别是那些与疾病密切相关的蛋白质。
• 生物标志物发现：显著 pQTL 可以帮助发现与特定表型（如疾病风险）相关的蛋白质生物标志物，用于疾病诊断或预后。
• 个性化医疗：了解遗传变异如何影响蛋白质表达，能够为个性化治疗和药物剂量优化提供信息。

eQTL

表达定量性状基因座（eQTL）是定量性状基因座（QTL）的一种，是与特定可量化性状的表型变异相关的基因组基因座（DNA 区域）。 虽然 QTL 一词可以指各种表型性状，但更具体的 eQTL 是指通过基因表达（如 mRNA 水平）测量的性状。虽然被命名为 “表达 QTL”，但并非所有的基因表达测量都可用于 eQTL。 例如，通过蛋白质水平量化的性状被称为蛋白质 QTLs（pQTLs）。

eQTL：表达定量性状位点是指影响基因表达水平的特定遗传位点，这些位点通常是指染色体上的特定位置，通常通过SNP这种遗传变异形式存在，即eQTL表示的是遗传变异与基因表达水平之间的关联。

• 定量形状位点（QTL）：指染色体上某个影响表型的遗传位点，定量形状是指表现为连续变化的性状。
• 基因表达定量性状（eQT）：基因表达水平通过RNA-seq等技术测量，表达量是一个定量的性状，由基因型和环境共同影响的。
• eQTL即负责调控基因表达水平的遗传位点，遗传变异通过影响这些位点从而改变某个基因的表达水平。

eQTL可以分为：1、cis-eQTL：当eQTL位点靠近它所调控的基因时，被称为cis-eQTL，通常意味着这个变异可能直接作用于基因的启动子区域、增强子区域或其他调控元件，直接影响该基因的表达；2、trans-eQTL：如果eQTL位点远离它调控的基因，甚至在不同染色体上，则被称为trans-eQTL，这种情况下，遗传变异可能通过更复杂的机制间接影响基因表达。

mQTL

mQTL（methylation Quantitative Trait Loci）分析是一种用来研究DNA甲基化与基因组中的遗传变异之间关联的分析方法。其核心目的是通过寻找基因组中影响DNA甲基化水平的遗传变异位点，揭示遗传因素对表观遗传修饰的调控机制。

在mQTL分析中，研究人员通常会检测个体的全基因组甲基化水平（比如CpG位点的甲基化程度），并与其全基因组的单核苷酸多态性（SNPs）数据进行关联分析。这种关联分析能够帮助识别哪些SNP位点（称为mQTL）与特定基因区域的甲基化水平显著相关。

mQTL分析的应用

1. 遗传调控的理解：mQTL分析揭示了遗传变异如何通过调控DNA甲基化影响基因表达，从而帮助我们理解基因调控网络和遗传调控机制。
2. 疾病研究：许多复杂疾病（如癌症、糖尿病等）都涉及遗传和表观遗传因素。mQTL分析有助于找到与疾病相关的甲基化调控位点，从而为疾病的发病机制提供线索。
3. 药物开发：通过识别mQTL，研究人员可以找到潜在的药物靶点，或设计靶向特定表观遗传修饰的药物。
4. 环境与遗传交互作用研究：mQTL研究也可以用于探讨环境因素（如吸烟、饮食）与遗传变异之间对DNA甲基化的影响，帮助理解环境和遗传交互对表观遗传状态的影响。

mQTL分析的主要步骤

1. 数据准备：收集个体的基因型数据（如SNP数据）和甲基化水平数据（常见的是CpG位点的甲基化水平）。
2. 关联分析：通过统计学方法（如线性回归或混合效应模型）对每个SNP-甲基化位点对进行关联分析，评估SNP是否对甲基化水平有显著影响。
3. 多重检验校正：由于分析过程中涉及大量的SNP和甲基化位点，需要进行多重检验校正（如Bonferroni校正或FDR方法），以减少假阳性结果。
4. 结果解释：对于显著关联的SNP-甲基化位点对，进一步分析其在生物学和疾病机制中的意义。

QTM

QTM分析，即“定量性状基因座作图”（Quantitative Trait Mapping，QTM），是一种用于定位控制复杂性状的基因区域（基因座）的分析方法。复杂性状（如身高、血压、产量等）往往由多个基因和环境因素共同影响，因此难以直接观察和分析。QTM分析通过基因组数据和性状表现数据之间的统计关联，寻找与这些性状相关的基因位点，通常用于农作物、牲畜、以及人类疾病的遗传学研究。

eQTM（表达定量性状基因座作图）

eQTM，即表达定量性状基因座作图，通过分析基因表达水平与基因组标记之间的关联来寻找控制基因表达的遗传变异。eQTM的目的是识别基因表达水平的遗传调控因素，帮助理解基因表达调控的遗传机制。

• 研究目标：识别基因表达水平（mRNA水平）的变异来源，确定调控基因表达的遗传位点（eQTLs）。
• 方法：在特定组织或细胞中测量大量基因的表达水平，将这些表达水平与基因组标记数据进行关联分析，常用的方法包括eQTL分析。
• 应用：eQTM广泛应用于疾病遗传学研究，通过识别与基因表达相关的基因变异，揭示复杂性状或疾病的潜在调控机制。它可以为理解表型变异的分子基础提供线索，尤其是当基因表达的变化与疾病表型密切相关时。

pQTM

pQTM，即蛋白质定量性状基因座作图，主要用于研究蛋白质水平上的复杂性状。pQTM的目标是定位调控蛋白质丰度或活性的遗传位点，从而揭示蛋白质层面的遗传调控。

• 研究目标：确定控制蛋白质丰度、结构或功能的遗传变异，识别调控蛋白质丰度或功能的基因位点（pQTLs）。
• 方法：测量蛋白质丰度或功能性相关参数（如酶活性等），将其与基因组标记数据进行关联。数据常来源于质谱技术（如LC-MS/MS），并结合pQTL分析模型。
• 应用：pQTM可用于阐明蛋白质调控网络，揭示蛋白质丰度如何影响复杂性状，尤其在药物靶标识别、疾病分子机制解析中有重要意义。通过pQTM分析，可以探索特定蛋白质的调控途径，找到潜在的治疗靶点。

EWAS

EWAS（Epigenome-Wide Association Study，全基因组范围表观遗传关联研究）是一种用于探究特定表观遗传修饰（如DNA甲基化）与特定表型（例如疾病、环境因素或行为特征）之间关联的研究方法。其目的是在全基因组范围内筛查出与特定表型显著相关的表观遗传标记，类似于GWAS（全基因组关联研究）在基因变异与表型关联方面的作用。

EWAS分析的关键要素

1. 研究对象和表型选择：通常，EWAS分析会涉及两个或多个表型组（如健康对照组和疾病组），以便比较表观遗传修饰水平差异。
2. 表观遗传修饰数据：主要是DNA甲基化数据，通常采用高通量技术（如Illumina Infinium HumanMethylation450K或850K芯片）来测量全基因组范围内的CpG位点甲基化水平。
3. 数据分析：在统计分析中，每个CpG位点的甲基化水平都会与表型变量进行关联分析。线性回归模型和混合效应模型是常用的方法，以控制年龄、性别、技术批次等混杂因素的影响。
4. 多重检验校正：由于EWAS分析涉及成千上万的CpG位点，所以需要进行多重检验校正（如Bonferroni校正或FDR校正）来降低假阳性风险。
5. 结果解释与生物学验证：对显著关联的CpG位点进一步解读，通常还需结合功能性注释分析，或在独立人群中验证结果的重复性和可靠性。

EWAS的应用

EWAS分析广泛应用于多种复杂性疾病和环境研究领域，包括但不限于以下几个方面：

1. 疾病机制研究：通过EWAS可以识别特定疾病（如癌症、糖尿病、阿尔茨海默病等）相关的甲基化改变，从而揭示疾病的表观遗传机制。
2. 环境因素研究：EWAS用于研究环境暴露（如吸烟、污染、饮食）对DNA甲基化的影响，探索环境因素如何通过表观遗传机制影响健康。
3. 老化和发育：EWAS也应用于探索表观遗传在老化和发育过程中的角色，帮助理解表观遗传随着年龄的动态变化及其对健康的影响。
4. 个体化医学：通过找到与特定疾病或表型相关的甲基化位点，可以进一步开发表观遗传标记用于疾病预测、早期筛查和诊断。

PheWAS

PheWAS从基因变异出发，寻找与之相关的多种表型（疾病或性状），首先选择特定的基因变异，然后系统地检查这些变异与数千种可能的表现之间的关系。适合用于发现基因的多效性，即一个基因变异可能影响多个不同的表型。

一句话来说，PheWAS关注的是特定基因变异与多种表型之间的关系。

GWAS

全基因组关联研究研究整个基因组范围内的遗传变异（包括编码和非编码区域），主要研究常见变异（MAF>1%），GWAS可以发现包括调控区域在内的更广泛的变异，需要对整个基因组进行测序或基因分型，成本较高。

一句话来说，GWAS关注的是整个基因组范畴的遗传变异与多种表型之间的关系。

ExWAS

全外显子关联研究专注于研究人类基因组中编码蛋白质的区域（即外显子）中的变异与特定表型之间的关联。它只关注基因组中的编码区域（可以研究罕见变异和常见变异，特别适合研究罕见致病变异），这些区域直接参与蛋白质的合成。相比于全基因组关联研究，ExWAS能够更有效地识别可能直接影响蛋白质功能的变异。

一句话来说，ExWAS关注的是外显子区域的所有变异与特定表型之间的关系