基于LASSO或ridge regression的Cox-Ph分析

介绍

glmnet提供了LASSO或ridge regression的Cox-PH分析模式,用于研究预测变量与生存时间的关系。具体做法先训练模型再使用最佳参数构建模型,最后评估预测精确度 C-index值判断预测的优劣。

加载数据

1
2
3
4
5
6
7
library(glmnet)
library(survival)
data(CoxExample)

phen <- y
rownames(phen) <- paste0("S", c(1:nrow(phen)))
head(phen) # 行样本名字,列是生存时间和状态
1
2
3
4
5
6
7
         time status
S1 1.76877757 1
S2 0.54528404 1
S3 0.04485918 0
S4 0.85032298 0
S5 0.61488426 1
S6 0.29860939 0
1
2
3
4
prof <- x
rownames(prof) <- rownames(phen)
colnames(prof) <- paste0("Feature", c(1:ncol(prof)))
head(prof)
1
2
3
4
5
6
7
    Feature1   Feature2    Feature3   Feature4    Feature5   Feature6    Feature7   Feature8    Feature9  Feature10
S1 -0.8767670 -0.6135224 -0.56757380 0.6621599 1.82218019 -1.0906678 -0.33186564 3.6754612 0.24580798 1.1382203
S2 -0.7463894 -1.7519457 0.28545898 1.1392105 0.80178007 1.8501985 0.30663005 -1.3729036 -0.03249051 0.7477848
S3 1.3759148 -0.2641132 0.88727408 0.3841870 0.05751801 -1.0917341 0.82119791 2.2960618 -0.44769567 -0.3046003
S4 0.2375820 0.7859162 -0.89670281 -0.8339338 -0.58237643 0.1874136 -0.58595131 0.4762090 -0.60580025 -1.2703322
S5 0.1086275 0.4665686 -0.57637261 1.7041314 0.32750715 -0.1211972 0.88537209 0.4505604 0.58878157 0.5504976
S6 1.2027213 -0.4187073 -0.05735193 0.5948491 0.44328682 -0.1191545 0.08097645 0.1645867 0.35648515 0.7186709

训练参数

设置alpha=0是ridge regression; alpha=1是LASSO;设置lambda,nlambda = 100或者lambda = 10^seq(3, -2, by = -.1);family选择不同数据分布情况,cox是cox-ph(其他选择”gaussian”, “binomial”, “poisson”, “multinomial”, “mgaussian”符合不同的数据)。

1
2
3
4
5
6
7
8
set.seed(123)
cv.fit <- cv.glmnet(prof, phen,
family = "cox",
type.measure = "C",
nfolds = 10,
alpha = 0,
nlambda = 100)
plot(cv.fit)

y轴坐标C-index原名是Harrell’concordance index,是用于评估模型的预测精度,常用于临床研究。x轴是lambda的log化结果,我们常选择最小的lambda值作为建模参数,也即是途中最大的C-index值。

C-index的计算方法是把所研究的资料中的所有研究对象随机地两两组成对子,以生存分析为例,两个病人如果生存时间较长的一位其预测生存时间长于另一位,或预测的生存概率高的一位的生存时间长于另一位,则称之为预测结果与实际结果相符,称之为一致。

计算C-index=K/M。

从上述计算方法可以看出C-index在0.5-1之间(任意配对随机情况下一致与不一致刚好是0.5的概率)。0.5为完全不一致,说明该模型没有预测作用,1为完全一致,说明该模型预测结果与实际完全一致。一般情况下C-index在0.50-0.70为准确度较低:在0.71-0.90之间为准确度中等;而高于0.90则为高准确度,跟相关系数有点类似。

构建模型

选择cv.fit最小lambda值

1
2
3
4
5
fit <- glmnet(prof, phen, 
family = "cox",
alpha = 0,
lambda = cv.fit$lambda.min)
summary(fit)

1
2
3
4
5
6
7
8
9
10
11
12
13
          Length Class     Mode   
a0 0 -none- NULL
beta 30 dgCMatrix S4
df 1 -none- numeric
dim 2 -none- numeric
lambda 1 -none- numeric
dev.ratio 1 -none- numeric
nulldev 1 -none- numeric
npasses 1 -none- numeric
jerr 1 -none- numeric
offset 1 -none- logical
call 6 -none- call
nobs 1 -none- numeric

  • 每个features的coefficient
    1
    coef(fit)
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    30 x 1 sparse Matrix of class "dgCMatrix"
    s0
    Feature1 0.5114403878
    Feature2 -0.1954830557
    Feature3 -0.2405974927
    Feature4 0.1957977902
    Feature5 -0.2074132864
    Feature6 -0.5056236002
    Feature7 0.3552934341
    Feature8 0.1057532384
    Feature9 0.4648827155
    Feature10 0.1375101610
    Feature11 -0.0194344395
    Feature12 0.0047078816
    Feature13 0.0410461245
    Feature14 0.0021407848
    Feature15 -0.0009016771
    Feature16 0.0001994629
    Feature17 -0.0372298071
    Feature18 -0.0100297637
    Feature19 0.0080887542
    Feature20 -0.0001315014
    Feature21 -0.0218432109
    Feature22 -0.0238062971
    Feature23 -0.0054209272
    Feature24 -0.0067311829
    Feature25 -0.0488808852
    Feature26 -0.0040028665
    Feature27 0.0286530927
    Feature28 -0.0136744813
    Feature29 0.0086380442
    Feature30 -0.0336064180
    可根据系数选择重要的features进行后续Cox-PH分析。

计算样本C-index

最佳lambda参数下构建模型的C-index值,反应模型的预测精度,越高越好

1
2
pred <- predict(fit, newx = prof)
apply(pred, 2, Cindex, y=phen)

1
2
       s0 
0.7344005

参考

  1. Regularized Cox Regression

  2. How and when: ridge regression with glmnet

  3. 一致性指数:Harrell’concordance index:C-index

  4. Calculate Concordance index

  5. Calculate Brier score on glmnet object

参考文章如引起任何侵权问题,可以与我联系,谢谢。

参考文章如引起任何侵权问题,可以与我联系,谢谢。


------------- The End Thanks for reading --------