大招：召唤最优的多因素cox模型
AIC信息准则即Akaike information criterion [1] ，是衡量统计模型拟合优良性(Goodness of fit)的一种标准，由于它为日本统计学家赤池弘次创立和发展的，因此又称赤池信息量准则。它建立在熵的概念基础上，可以权衡所估计模型的复杂度和此模型拟合数据的优良性。
赤池信息准则的方法是寻找可以最好地解释数据但包含最少自由参数的模型。
-摘自百度百科
可以使用这一准则得到最优的多因素cox模型，以下是代码实现
1.输入数据
由生存信息与基因表达量组成的数据框。
load("dat.Rdata")
head(dat)
## sample event time
## TCGA-06-0878-01A TCGA-06-0878-01A 0 218
## TCGA-26-5135-01A TCGA-26-5135-01A 1 270
## TCGA-06-5859-01A TCGA-06-5859-01A 0 139
## TCGA-06-2563-01A TCGA-06-2563-01A 0 932
## TCGA-41-2571-01A TCGA-41-2571-01A 1 26
## TCGA-28-5207-01A TCGA-28-5207-01A 1 343
## CYP2E1 MYH7 CASP4
## TCGA-06-0878-01A 0.11449559 0.8049116 2.7651858
## TCGA-26-5135-01A 0.05534342 2.1560713 2.1794106
## TCGA-06-5859-01A 0.03202620 1.2187957 2.3940125
## TCGA-06-2563-01A 0.31261543 0.9101483 2.7048505
## TCGA-41-2571-01A 0.32700857 2.2878184 0.7237877
## TCGA-28-5207-01A 0.15228157 1.6271620 1.9986231
## SIRT1 ATP2A2 PPM1L
## TCGA-06-0878-01A 2.132754 4.931262 1.405607
## TCGA-26-5135-01A 1.947894 5.054068 1.661491
## TCGA-06-5859-01A 2.272596 4.516211 1.345576
## TCGA-06-2563-01A 1.990387 4.879086 1.283925
## TCGA-41-2571-01A 1.655961 4.753264 1.265152
## TCGA-28-5207-01A 2.469023 4.530683 2.079903
## CP BAG1 DNAJB11
## TCGA-06-0878-01A 4.729828 2.981576 4.489935
## TCGA-26-5135-01A 1.514229 3.299911 3.934268
## TCGA-06-5859-01A 3.009394 3.197839 3.360617
## TCGA-06-2563-01A 1.262133 3.533318 3.828734
## TCGA-41-2571-01A 1.410165 3.408114 3.083744
## TCGA-28-5207-01A 1.442487 3.488568 3.558798
## MMP2 SLN FKBP14
## TCGA-06-0878-01A 5.196985 3.2520898 2.710485
## TCGA-26-5135-01A 5.267932 3.3372085 2.295893
## TCGA-06-5859-01A 3.364497 3.6417199 3.461916
## TCGA-06-2563-01A 6.471976 0.9712369 2.666084
## TCGA-41-2571-01A 8.445514 2.8842578 2.158597
## TCGA-28-5207-01A 4.552452 0.9862039 2.322108
## ERLIN1 GPX7 EPM2A
## TCGA-06-0878-01A 2.963291 4.161208 1.824520
## TCGA-26-5135-01A 2.027252 3.812280 1.513929
## TCGA-06-5859-01A 2.702454 4.032823 1.696385
## TCGA-06-2563-01A 2.690188 4.270150 2.174987
## TCGA-41-2571-01A 1.892868 3.539653 1.091940
## TCGA-28-5207-01A 2.802503 3.302680 2.323285
## BRCA1 HFE CDKN3
## TCGA-06-0878-01A 1.929324 2.3744248 2.134364
## TCGA-26-5135-01A 1.535608 0.7702822 2.929064
## TCGA-06-5859-01A 3.118322 2.5060929 2.464366
## TCGA-06-2563-01A 1.324959 1.6765673 2.930530
## TCGA-41-2571-01A 1.082552 0.7707545 3.408663
## TCGA-28-5207-01A 1.958569 1.7947694 2.014698
## SERPINE1 ERP27 NOL3
## TCGA-06-0878-01A 7.693198 0.7583165 3.954610
## TCGA-26-5135-01A 6.229662 0.3775101 3.521318
## TCGA-06-5859-01A 2.814878 0.1522165 2.896874
## TCGA-06-2563-01A 4.167370 0.4890798 1.931305
## TCGA-41-2571-01A 6.717429 0.3638483 3.889378
## TCGA-28-5207-01A 2.902920 0.4894604 2.013993
## BET1 SNCA DNAJB12
## TCGA-06-0878-01A 2.649736 1.750464 2.653257
## TCGA-26-5135-01A 1.836402 2.655945 2.878896
## TCGA-06-5859-01A 2.601570 1.341060 2.979048
## TCGA-06-2563-01A 2.543613 2.339013 2.874104
## TCGA-41-2571-01A 2.465486 2.155570 2.898581
## TCGA-28-5207-01A 3.207056 3.431305 3.037193
## CASP2 GLA RNF185
## TCGA-06-0878-01A 3.056804 3.699737 3.554046
## TCGA-26-5135-01A 2.837734 3.385206 4.087988
## TCGA-06-5859-01A 2.662097 4.166703 4.399242
## TCGA-06-2563-01A 3.175957 3.787763 4.394133
## TCGA-41-2571-01A 3.262038 3.945768 4.358992
## TCGA-28-5207-01A 3.147226 3.779798 3.717130
## MX1 MAPK3 ANXA5
## TCGA-06-0878-01A 3.147325 4.639285 7.777013
## TCGA-26-5135-01A 2.687709 5.064507 7.067037
## TCGA-06-5859-01A 4.173986 4.214156 6.306859
## TCGA-06-2563-01A 4.910899 4.527579 7.995120
## TCGA-41-2571-01A 1.514851 4.887269 6.413893
## TCGA-28-5207-01A 3.614883 4.950130 7.653079
## SHISA5 ERP29 SREBF2
## TCGA-06-0878-01A 4.891004 5.090190 3.646879
## TCGA-26-5135-01A 6.280796 5.088310 4.656178
## TCGA-06-5859-01A 4.708155 5.133025 4.619645
## TCGA-06-2563-01A 5.343170 5.071684 4.689714
## TCGA-41-2571-01A 5.184783 4.741366 4.736130
## TCGA-28-5207-01A 5.152544 5.321610 4.002934
## DDIT3 MR1 PRNP
## TCGA-06-0878-01A 5.599779 2.8585494 6.131223
## TCGA-26-5135-01A 5.330168 1.5743079 6.384507
## TCGA-06-5859-01A 5.621570 3.3338802 6.554235
## TCGA-06-2563-01A 4.388975 2.0166614 6.122956
## TCGA-41-2571-01A 5.840359 0.8671249 6.432797
## TCGA-28-5207-01A 7.208496 2.8783165 6.527731
使用这36个基因，用基于AIC的逐步回归法挑选基因构建cox模型
2.stepwise走起来
library(My.stepwise)
vl <- colnames(dat)[-(1:3)]
My.stepwise.coxph(Time = "time",
Status = "event",
variable.list = vl,
data = dat)
会看到输出了一大片，是进行了多次迭代，目的是为了得出最好的模型，输出结果有几千行，我就不展示了。
3.输出结果
关于这个包的教程资料很少，我没有找到对结果如何理解的描述，所以就自己理解了一下，到底是所有模型都能用，还是最后一个是最好的呢？试一试呗。
这是从输出的一大片信息中复制的最后几个模型公式
library(survival)
model1 = coxph(formula = Surv(time, event) ~ HFE + SHISA5 + BRCA1 + RNF185 +
EPM2A + ERLIN1 + GPX7 + SLN + DNAJB11 + MMP2 + NOL3 + CP +
CYP2E1 + ATP2A2 + GLA + MAPK3, data = dat)
model2 =coxph(formula = Surv(time, event) ~ HFE + SHISA5 + BRCA1 + EPM2A +
ERLIN1 + GPX7 + SLN + DNAJB11 + MMP2 + NOL3 + CP + ATP2A2 +
GLA + MAPK3 + SREBF2 + CASP2 + SNCA + DDIT3 + BAG1, data = dat)
model3 = coxph(formula = Surv(time, event) ~ HFE + SHISA5 + BRCA1 + EPM2A +
ERLIN1 + GPX7 + SLN + DNAJB11 + MMP2 + NOL3 + CP + ATP2A2 +
GLA + MAPK3 + SREBF2 + CASP2 + SNCA + DDIT3 + BAG1 + ANXA5,
data = dat)
可以看看他们的c-index值和AIC值
summary(model1)$concordance
## C se(C)
## 0.86042967 0.01027966
summary(model2)$concordance
## C se(C)
## 0.864850308 0.009782714
summary(model3)$concordance
## C se(C)
## 0.865988496 0.009754695
broom::glance(model1)$AIC
## [1] 2615.754
broom::glance(model2)$AIC
## [1] 2611.273
broom::glance(model3)$AIC
## [1] 2611.064
随着迭代次数的增加，模型的c-index值不断增加，AIC值不断减小，这意味着从开始到最后，模型是越来越好，所以输出的那一大片信息里，最后一个就是最优模型啦。
4.问题解决
1.发现问题
收到反馈说这个代码跑不动了。会出warning且模型无法得到。
2. 解决问题
代码不曾修改，发现R包也不曾更新。这就奇怪了，原来能跑动的代码现在瓦特了嘛？
我的讲课电脑重新安装了R和R包，跑代码时出现了和学员一样的问题。然后拿出另一台电脑，却能跑通。然后对比了一下这几个依赖包的版本，发现只有survival是不同的。
所以解决办法是：安装3.3-1版本的survival。重装之后恢复正常啦。
https://cran.r-project.org/web/packages/survival/index.html
在这里面找3.3-1，下载到工作目录下，装即可。
install.packages("survival_3.3-1.tar.gz",repos = NULL)
共有 0 条评论