Cox 比例风险模型中HR和置信区间

Cox 回归是一种用于生存分析的统计模型，它可以用来估计某个因素对生存时间的影响。Cox 回归基于 Cox 比例风险模型，该模型假设风险比率是常数，即不随时间变化。在 Cox 回归中，我们使用最大似然方法来估计模型参数。

Cox 回归的计算过程：

1. 首先，我们需要准备数据集。数据集应包含每个观察值的生存时间、是否发生事件（例如死亡或治愈）以及其他可能影响生存时间的因素（例如年龄、性别、疾病状态等）。

2. 接下来，我们需要将数据集拆分为两个子集：训练集和测试集。训练集用于拟合 Cox 模型，测试集用于评估模型性能。

3. 对于每个观察值 $i$，我们定义一个风险函数 h i ( t ) h_i(t) hi​(t)，表示在时刻 $t$ 发生事件的概率。Cox 模型假设风险函数可以表示为以下形式：

   h i ( t ) = h 0 ( t ) ⋅ e ( β 1 x i 1 + β 2 x i 2 + . . . + β p x i p ) h_i(t) = h_0(t) \\cdot e^{(\\beta_1 x_{i1} + \\beta_2 x_{i2} + ... + \\beta_p x_{ip})} hi​(t)=h0​(t)⋅e(β1​xi1​+β2​xi2​+...+βp​xip​)

   其中 h 0 ( t ) h_0(t) h0​(t) 是基准风险函数，表示在所有协变量都为 0 时的风险。 β 1 , β 2 , . . . , β p \\beta_1, \\beta_2, ..., \\beta_p β1​,β2​,...,βp​ 是模型参数，表示每个协变量对风险的影响。 x i 1 , x i 2 , . . . , x i p x_{i1}, x_{i2}, ..., x_{ip} xi1​,xi2​,...,xip​ 是观察值 $i$ 的协变量值。

4. 我们使用最大似然方法来估计模型参数。最大似然方法的目标是找到一组参数，使得观察到的数据在这些参数下的概率最大。对于 Cox 回归，我们需要最大化以下似然函数：

   L ( β 1 , β 2 , . . . , β p ) = ∏ i = 1 n e ( β 1 x i 1 + β 2 x i 2 + . . . + β p x i p ) ∑ j ∈ R ( t i ) e ( β 1 x j 1 + β 2 x j 2 + . . . + β p x j p ) L(\\beta_1, \\beta_2, ..., \\beta_p) = \\prod_{i=1}^{n} \\frac{e^{(\\beta_1 x_{i1} + \\beta_2 x_{i2} + ... + \\beta_p x_{ip})}}{\\sum_{j\\in R(t_i)} e^{(\\beta_1 x_{j1} + \\beta_2 x_{j2} + ... + \\beta_p x_{jp})}} L(β1​,β2​,...,βp​)=i=1∏n​∑j∈R(ti​)​e(β1​xj1​+β2​xj2​+...+βp​xjp​)e(β1​xi1​+β2​xi2​+...+βp​xip​)​

   其中 $n$ 是观察值的数量， R ( t i ) R(t_i) R(ti​) 是在时刻 t i t_i ti​ 仍然存活的观察值的集合。

5. 我们可以使用梯度下降等优化算法来最大化似然函数，并找到最优参数 β ^ 1 , β ^ 2 , . . . , β ^ p \\hat{\\beta}_1, \\hat{\\beta}_2, ..., \\hat{\\beta}_p β^​1​,β^​2​,...,β^​p​。

6. 一旦我们找到了最优参数，我们就可以使用 Cox 模型来预测新观察值的生存时间。对于一个新观察值 x ∗ x^* x∗，我们可以计算其风险比率为：

   h ( x ∗ ) h 0 ( t ) = e ( β ^ 1 x 1 ∗ + β ^ 2 x 2 ∗ + . . . + β ^ p x p ∗ ) \\frac{h(x^*)}{h_0(t)} = e^{(\\hat{\\beta}_1 x^*_1 + \\hat{\\beta}_2 x^*_2 + ... + \\hat{\\beta}_p x^*_p)} h0​(t)h(x∗)​=e(β^​1​x1∗​+β^​2​x2∗​+...+β^​p​xp∗​)

   其中 x 1 ∗ , x 2 ∗ , . . . , x p ∗ x^*_1, x^*_2, ..., x^*_p x1∗​,x2∗​,...,xp∗​ 是新观察值的协变量值。

Cox 回归的计算过程比较复杂，需要对统计学和数学有一定的了解。在实践中，我们通常使用现成的软件包来拟合 Cox 模型，并进行预测和评估。常用的 Cox 回归软件包包括 R 中的 survival 包和 Python 中的 lifelines 包。

95CI计算方式

在 Cox 回归中，我们通常使用最大似然方法来估计模型参数。最大似然方法可以为每个参数提供一个点估计值，但是它并没有提供一个置信区间来描述这个点估计值的不确定性。因此，我们需要使用一些方法来计算参数的置信区间。

Cox 回归中常用的方法是基于 Wald 统计量的置信区间。Wald 统计量是一个用于检验模型参数是否显著的统计量，其定义为：

W = ( β ^ − β 0 ) S E ( β ^ ) W = \\frac{(\\hat{\\beta} - \\beta_0)}{SE(\\hat{\\beta})} W=SE(β^​)(β^​−β0​)​

其中 β ^ \\hat{\\beta} β^​ 是模型参数的点估计值， β 0 \\beta_0 β0​ 是要检验的假设值（通常为 0）， S E ( β ^ ) SE(\\hat{\\beta}) SE(β^​) 是 β ^ \\hat{\\beta} β^​ 的标准误差。

根据正态分布理论，当样本容量足够大时，Wald 统计量近似服从标准正态分布。因此，我们可以使用标准正态分布表来查找相应的临界值，并将其乘以标准误差来得到置信区间。

具体地说，在 Cox 回归中，我们可以使用以下公式来计算 95% 置信区间：

β ^ ± 1.96 ⋅ S E ( β ^ ) \\hat{\\beta} \\pm 1.96 \\cdot SE(\\hat{\\beta}) β^​±1.96⋅SE(β^​)

其中 1.96 是标准正态分布在 95% 置信水平下的临界值。

需要注意的是，Wald 置信区间有一些局限性。当样本容量较小时，Wald 置信区间可能过于宽泛或过于狭窄，导致估计结果不准确。此外，当模型参数存在较大偏差时，Wald 置信区间也可能不准确。因此，在实践中，我们通常会使用其他方法来计算参数的置信区间，例如基于 Bootstrap 的置信区间和基于贝叶斯统计的置信区间。

Bootstrap 置信区间是一种非参数方法，它通过从原始数据集中重复抽样来估计参数的分布，并使用这个分布来计算置信区间。具体地说，在 Cox 回归中，我们可以使用以下步骤来计算 Bootstrap 置信区间：

1. 从原始数据集中随机抽取一个样本，将其作为新的数据集。

2. 使用 Cox 回归模型拟合新的数据集，并记录模型参数。

3. 重复步骤 1 和步骤 2 多次（例如 1000 次），得到模型参数的分布。

4. 根据模型参数的分布计算置信区间。通常情况下，我们可以使用第 2.5% 和第 97.5% 分位数作为 95% 置信区间的下限和上限。

Bootstrap 置信区间不依赖于任何假设，因此它对样本容量和模型偏差不敏感。但是，它需要进行大量重复抽样，并且可能需要较长时间才能得到结果。

贝叶斯置信区间是一种基于贝叶斯统计理论的方法，它使用先验分布和后验分布来描述参数的不确定性，并计算置信区间。贝叶斯置信区间可以提供更准确的估计结果，但需要选择合适的先验分布，并进行复杂的计算。

在实践中，我们通常会使用多种方法来计算参数的置信区间，并比较它们的结果。如果不同方法得到的置信区间相似，则可以认为估计结果比较可靠。如果不同方法得到的置信区间差异较大，则需要进一步检查数据和模型，以确定哪种方法更适合当前的问题。

需要注意的是，置信区间只是对参数估计结果的不确定性进行了量化，并不能保证参数估计结果的准确性。在解释参数估计结果时，我们应该考虑到置信区间以外的可能性，并结合实际情况进行分析。此外，我们还应该注意避免过度解读统计显著性和置信区间，而应该将它们作为辅助工具来帮助我们理解数据和模型。

python中lifelines 库的CoxPHFitter类HR值计算方式

在 Python 的 lifelines 库中，CoxPHFitter 类可以用来拟合 Cox 比例风险模型，并计算相应的风险比（hazard ratio，HR）和置信区间。具体地说，在 CoxPHFitter 中，HR 的计算方法如下：

1. 对于每个自变量 X i X_i Xi​，计算其对应的系数 β i \\beta_i βi​。

2. 计算每个自变量 X i X_i Xi​ 的指数函数 e β i e^{\\beta_i} eβi​。

3. 计算两组个体在某个时间点 $t$ 的风险比（hazard ratio）：

   H R ( t ) = h 1 ( t ) h 0 ( t ) = ∏ i = 1 k e β i X i , 1 ∏ i = 1 k e β i X i , 0 = e ∑ i = 1 k β i ( X i , 1 − X i , 0 ) HR(t) = \\frac{h_1(t)}{h_0(t)} = \\frac{\\prod_{i=1}^{k} e^{\\beta_i X_{i,1}}}{\\prod_{i=1}^{k} e^{\\beta_i X_{i,0}}} = e^{\\sum_{i=1}^{k}\\beta_i(X_{i,1}-X_{i,0})} HR(t)=h0​(t)h1​(t)​=∏i=1k​eβi​Xi,0​∏i=1k​eβi​Xi,1​​=e∑i=1k​βi​(Xi,1​−Xi,0​)

   其中 $k$ 是自变量的数量， X i , 1 X_{i,1} Xi,1​ 和 X i , 0 X_{i,0} Xi,0​ 分别是第一组和第二组个体在自变量 X i X_i Xi​ 上的取值。

4. 计算 HR 的置信区间。通常情况下，我们可以使用 Wald 置信区间或 Bootstrap 置信区间来计算 HR 的置信区间。

需要注意的是，在 CoxPHFitter 中，HR 的计算依赖于 Cox 比例风险模型的假设，即风险比在时间上是恒定的。如果这个假设不成立，例如存在时间依赖性或非比例风险等情况，那么 CoxPHFitter 得到的 HR 可能不准确。在实践中，我们应该对数据和模型进行充分的检查和验证，并选择合适的方法来处理非比例风险或时间依赖性等问题。