数据分析：基于glmnet的Cox

glmnet提供了LASSO或ridge regression的Cox-PH分析模式，用于研究预测变量与生存时间的关系。

加载数据

library(glmnet)library(survival)data(CoxExample)phen <- y rownames(phen) <- paste0("S", c(1:nrow(phen)))head(phen)  # 行样本名字，列是生存时间和状态
time statusS1 1.76877757      1S2 0.54528404      1S3 0.04485918      0S4 0.85032298      0S5 0.61488426      1S6 0.29860939      0
prof <- xrownames(prof) <- rownames(phen)colnames(prof) <- paste0("Feature", c(1:ncol(prof)))head(prof)
Feature1   Feature2    Feature3   Feature4    Feature5   Feature6    Feature7   Feature8    Feature9  Feature10S1 -0.8767670 -0.6135224 -0.56757380  0.6621599  1.82218019 -1.0906678 -0.33186564  3.6754612  0.24580798  1.1382203S2 -0.7463894 -1.7519457  0.28545898  1.1392105  0.80178007  1.8501985  0.30663005 -1.3729036 -0.03249051  0.7477848S3  1.3759148 -0.2641132  0.88727408  0.3841870  0.05751801 -1.0917341  0.82119791  2.2960618 -0.44769567 -0.3046003S4  0.2375820  0.7859162 -0.89670281 -0.8339338 -0.58237643  0.1874136 -0.58595131  0.4762090 -0.60580025 -1.2703322S5  0.1086275  0.4665686 -0.57637261  1.7041314  0.32750715 -0.1211972  0.88537209  0.4505604  0.58878157  0.5504976S6  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"符合不同的数据）。

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值

fit <- glmnet(prof, phen,               family = "cox",               alpha = 0,              lambda = cv.fit$lambda.min)summary(fit)
Length Class     Mode   a0         0     -none-    NULL   beta      30     dgCMatrix S4     df         1     -none-    numericdim        2     -none-    numericlambda     1     -none-    numericdev.ratio  1     -none-    numericnulldev    1     -none-    numericnpasses    1     -none-    numericjerr       1     -none-    numericoffset     1     -none-    logicalcall       6     -none-    call   nobs       1     -none-    numeric

  • 每个features的coefficient

coef(fit)
30 x 1 sparse Matrix of class "dgCMatrix"                     s0Feature1   0.5114403878Feature2  -0.1954830557Feature3  -0.2405974927Feature4   0.1957977902Feature5  -0.2074132864Feature6  -0.5056236002Feature7   0.3552934341Feature8   0.1057532384Feature9   0.4648827155Feature10  0.1375101610Feature11 -0.0194344395Feature12  0.0047078816Feature13  0.0410461245Feature14  0.0021407848Feature15 -0.0009016771Feature16  0.0001994629Feature17 -0.0372298071Feature18 -0.0100297637Feature19  0.0080887542Feature20 -0.0001315014Feature21 -0.0218432109Feature22 -0.0238062971Feature23 -0.0054209272Feature24 -0.0067311829Feature25 -0.0488808852Feature26 -0.0040028665Feature27  0.0286530927Feature28 -0.0136744813Feature29  0.0086380442Feature30 -0.0336064180

可根据系数选择重要的features进行后续Cox-PH分析。

计算样本C-index

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

pred <- predict(fit, newx = prof)apply(pred, 2, Cindex, y=phen)
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

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

赞 (0)

相关推荐