大招：召唤最优的多因素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值不断减小，这意味着从开始到最后，模型是越来越好，所以输出的那一大片信息里，最后一个就是最优模型啦。