rm(list = ls())
load(file = "datasets/lnc_expr_clin.RData")
#去掉没有生存信息的样本
lnc_expr_clin1 <- lnc_expr_clin[!is.na(lnc_expr_clin$time_months),]
lnc_expr_clin1 <- lnc_expr_clin1[lnc_expr_clin1$time_months>0,]
#选择其中一部分数据
dat.cox <- lnc_expr_clin1[,c(72:73,1:59)]
dim(dat.cox)
## [1] 297 61
dat.cox[1:4,1:6]
## event time_months PGM5-AS1 LINC01082 AC005180.2 AC005180.1
## 1 0 36.33 0.15064007 0.2642238 0.0000000 0.1547768
## 2 0 13.87 0.06309362 0.1666554 0.3105983 0.2436603
## 3 1 21.83 2.16399508 3.5662920 2.2454129 2.0073496
## 4 0 18.20 2.73075081 1.7314314 0.8609916 0.732301416 变量选择之逐步回归
筛选变量的方法多如牛毛，我们今天介绍逐步选择法。逐步选择有3种策略，分别是向前（forward）、向后（backward）、逐步法（stepwise）。
向前逐步选择从一个零特征模型开始，然后每次添加一个特征，直到所有特征添加完毕。在这个过程中，被添加的选定特征建立的模型具有最小的RSS（残差平方和）。所以理论上，第一个选定的特征应该能最好地解释响应变量，依此类推。
注意
添加一个特征一定会使RSS减少，使R方增加，但不一定能提高模型的拟合度和可解释性。
向后逐步回归从一个包含所有特征的模型开始，每次删除一个起最小作用的特征。
逐步法是一种混合方法，这种算法先通过向前逐步回归添加特征，然后检查是否有特征不再对提高模型拟合度起作用，如果有则删除。每次建模之后，分析者都可以检查模型输出，并使用各种统计量选择能提供最佳拟合的特征。
注意
逐步回归技术会遇到非常严重的问题。对于一个数据集，你先用向前逐步回归，然后再用向后逐步回归，可能会得到两个完全矛盾的模型。最重要的一点是，逐步回归会使回归系数发生偏离，换句话说，会使回归系数的值过大，置信区间过窄（Tibrishani，1996）。
逐步选择可以使用MASS包的stepAIC实现，如R语言和医学统计学系列（7）：多元线性回归所介绍的，也可以直接使用R语言自带的step函数。stepAIC比step功能多一点，但是对于简单的变量筛选来说这俩没区别。
16.1 加载数据
我们使用TCGA-BLCA的lncRNA数据（数据已放在粉丝qq群文件），其中包括408个样本，time_months是生存时间，event是生存状态，1代表死亡，0代表生存，其余变量都是自变量。
现在这个数据一共59个自变量，我们先使用所有自变量建立cox回归模型。
16.2 建立模型
我们这个是生存数据，使用cox回归。如果你的数据是其他类型，使用逻辑回归或者线性回归都是可以的。
library(survival)
fit.cox <- coxph(Surv(time_months,event)~.,data = dat.cox)
fit.cox
## Call:
## coxph(formula = Surv(time_months, event) ~ ., data = dat.cox)
##
## coef exp(coef) se(coef) z p
## `PGM5-AS1` -0.008183 0.991850 0.222738 -0.037 0.97069
## LINC01082 0.345614 1.412858 0.403674 0.856 0.39190
## AC005180.2 0.977584 2.658027 0.906672 1.078 0.28094
## AC005180.1 0.846348 2.331118 1.193460 0.709 0.47823
## FENDRR -0.653451 0.520247 0.548215 -1.192 0.23328
## AC053503.3 -0.589548 0.554578 0.553104 -1.066 0.28647
## MIR100HG 0.902471 2.465687 0.386229 2.337 0.01946
## AP001107.5 -0.812922 0.443560 0.677700 -1.200 0.23032
## `C5orf66-AS1` 0.094286 1.098873 0.295767 0.319 0.74989
## NR4A1AS 0.874631 2.397990 0.336703 2.598 0.00939
## AL162424.1 0.049223 1.050455 0.324490 0.152 0.87943
## AF001548.1 0.949594 2.584659 0.788709 1.204 0.22860
## AC099850.4 0.205430 1.228053 0.252745 0.813 0.41634
## `MBNL1-AS1` -0.798659 0.449932 0.675864 -1.182 0.23733
## `ADAMTS9-AS1` -0.065160 0.936917 1.776167 -0.037 0.97074
## MIR22HG 0.108393 1.114486 0.228433 0.475 0.63514
## MIR200CHG -0.070914 0.931542 0.215534 -0.329 0.74214
## AC093010.3 -0.658117 0.517825 0.347044 -1.896 0.05791
## LINC00865 -0.282616 0.753809 0.284716 -0.993 0.32089
## AP003071.4 -0.506228 0.602765 0.675734 -0.749 0.45377
## PCAT6 0.347869 1.416047 0.199663 1.742 0.08146
## LINC02657 -0.175082 0.839388 0.122482 -1.429 0.15287
## `PPP1R14B-AS1` 0.144657 1.155643 0.238269 0.607 0.54377
## AC012085.2 -2.267686 0.103552 1.671128 -1.357 0.17479
## `ACTA2-AS1` -0.084816 0.918681 0.760755 -0.111 0.91123
## AC036108.3 1.405644 4.078154 1.413552 0.994 0.32003
## AC079313.2 0.167743 1.182633 1.020173 0.164 0.86940
## AC020916.1 0.037921 1.038649 0.201470 0.188 0.85070
## SNHG25 -0.151295 0.859594 0.330536 -0.458 0.64715
## AL049555.1 0.398962 1.490277 0.207738 1.921 0.05479
## `MIR1-1HG-AS1` -1.851131 0.157059 1.991975 -0.929 0.35274
## AC018904.1 0.026484 1.026838 0.233038 0.114 0.90952
## SNHG12 0.144329 1.155264 0.344430 0.419 0.67519
## `SPINT1-AS1` 0.676775 1.967522 0.374980 1.805 0.07110
## `KRT7-AS` -0.137828 0.871248 0.129775 -1.062 0.28821
## MIR205HG -0.076826 0.926051 0.161261 -0.476 0.63378
## `HAND2-AS1` 1.853233 6.380415 1.742430 1.064 0.28751
## AL445524.1 -0.255737 0.774345 0.220099 -1.162 0.24527
## LINC01980 -0.171282 0.842584 0.128515 -1.333 0.18260
## `ZNF710-AS1` -0.959290 0.383165 0.459985 -2.085 0.03703
## AC092718.4 0.010577 1.010633 0.272577 0.039 0.96905
## AC008735.2 -0.002696 0.997308 0.291545 -0.009 0.99262
## LINC01133 0.122659 1.130499 0.120170 1.021 0.30739
## AC025575.2 0.158544 1.171804 0.135523 1.170 0.24205
## `MAFG-DT` 0.343519 1.409901 0.242305 1.418 0.15627
## CASC9 -0.118071 0.888633 0.154923 -0.762 0.44598
## AL390719.2 0.177187 1.193854 0.225404 0.786 0.43182
## AC002398.2 -1.089318 0.336446 1.662224 -0.655 0.51225
## AC008736.1 -0.117863 0.888818 0.163847 -0.719 0.47193
## AL161431.1 0.158805 1.172110 0.112004 1.418 0.15623
## `PCCA-DT` -0.456381 0.633572 0.254320 -1.795 0.07273
## AC245041.2 0.243686 1.275944 0.200371 1.216 0.22392
## U62317.1 -0.162513 0.850005 0.205601 -0.790 0.42928
## U62317.2 -0.131903 0.876426 0.315579 -0.418 0.67597
## `VPS9D1-AS1` -0.044547 0.956431 0.174172 -0.256 0.79813
## AL023284.4 -0.335339 0.715095 0.245641 -1.365 0.17220
## AATBC 0.136654 1.146432 0.180919 0.755 0.45005
## LINC00641 0.383262 1.467062 0.474750 0.807 0.41950
## AC015912.3 -0.562875 0.569569 0.296322 -1.900 0.05749
##
## Likelihood ratio test=78.86 on 59 df, p=0.04311
## n= 297, number of events= 71下面就是用逐步法选择变量。
16.3 逐步选择
我们使用逐步选择法进行变量筛选：
fit.step <- step(fit.cox,direction = "both")
#save(fit.step,file = "./datasets/fit.step.edata")
Start: AIC=794.12
Surv(time_months, event) ~ `PGM5-AS1` + LINC01082 + AC005180.2 +
AC005180.1 + FENDRR + AC053503.3 + MIR100HG + AP001107.5 +
`C5orf66-AS1` + NR4A1AS + AL162424.1 + AF001548.1 + AC099850.4 +
`MBNL1-AS1` + `ADAMTS9-AS1` + MIR22HG + MIR200CHG + AC093010.3 +
LINC00865 + AP003071.4 + PCAT6 + LINC02657 + `PPP1R14B-AS1` +
AC012085.2 + `ACTA2-AS1` + AC036108.3 + AC079313.2 + AC020916.1 +
SNHG25 + AL049555.1 + `MIR1-1HG-AS1` + AC018904.1 + SNHG12 +
`SPINT1-AS1` + `KRT7-AS` + MIR205HG + `HAND2-AS1` + AL445524.1 +
LINC01980 + `ZNF710-AS1` + AC092718.4 + AC008735.2 + LINC01133 +
AC025575.2 + `MAFG-DT` + CASC9 + AL390719.2 + AC002398.2 +
AC008736.1 + AL161431.1 + `PCCA-DT` + AC245041.2 + U62317.1 +
U62317.2 + `VPS9D1-AS1` + AL023284.4 + AATBC + LINC00641 +
AC015912.3
Df AIC
- AC008735.2 1 792.12
- `ADAMTS9-AS1` 1 792.12
- `PGM5-AS1` 1 792.12
- AC092718.4 1 792.12
- `ACTA2-AS1` 1 792.13
- AC018904.1 1 792.13
- AL162424.1 1 792.14
- AC079313.2 1 792.15
- AC020916.1 1 792.15
- `VPS9D1-AS1` 1 792.18
- `C5orf66-AS1` 1 792.22
- MIR200CHG 1 792.23
- U62317.2 1 792.29
- SNHG12 1 792.29
- SNHG25 1 792.33
- MIR22HG 1 792.34
- MIR205HG 1 792.35
- `PPP1R14B-AS1` 1 792.49
- AC002398.2 1 792.56
- AC005180.1 1 792.63
- AC008736.1 1 792.65
- AATBC 1 792.69
- AP003071.4 1 792.69
- CASC9 1 792.70
- AL390719.2 1 792.74
- U62317.1 1 792.75
- LINC00641 1 792.76
- AC099850.4 1 792.79
- LINC01082 1 792.86
- `MIR1-1HG-AS1` 1 793.03
- AC036108.3 1 793.10
- LINC00865 1 793.13
- LINC01133 1 793.13
- `HAND2-AS1` 1 793.23
- `KRT7-AS` 1 793.25
- AC005180.2 1 793.25
- AC053503.3 1 793.29
- AL445524.1 1 793.49
- AC025575.2 1 793.50
- AF001548.1 1 793.55
- FENDRR 1 793.61
- AC245041.2 1 793.61
- `MBNL1-AS1` 1 793.62
- AP001107.5 1 793.84
- LINC01980 1 793.87
- AL023284.4 1 794.01
- AL161431.1 1 794.11
<none> 794.12
- `MAFG-DT` 1 794.13
- AC012085.2 1 794.14
- LINC02657 1 794.20
- PCAT6 1 795.17
- `PCCA-DT` 1 795.43
- `SPINT1-AS1` 1 795.45
- AC093010.3 1 795.84
- AL049555.1 1 795.88
- AC015912.3 1 795.94
- `ZNF710-AS1` 1 796.86
- MIR100HG 1 797.13
- NR4A1AS 1 798.62
Step: AIC=792.12
## 省略巨多中间过程
Step: AIC=790.12
## 省略巨多中间过程
Step: AIC=734.72
Surv(time_months, event) ~ AC005180.2 + MIR100HG + AP001107.5 +
NR4A1AS + AC093010.3 + PCAT6 + AC036108.3 + AL049555.1 +
`MIR1-1HG-AS1` + `SPINT1-AS1` + LINC01980 + `ZNF710-AS1` +
AL161431.1 + `PCCA-DT` + U62317.1 + AL023284.4 + AC015912.3
Df AIC
<none> 734.72
+ LINC00641 1 734.91
+ AC012085.2 1 735.01
- AL049555.1 1 735.02
+ AC002398.2 1 735.06
- AC036108.3 1 735.12
+ `MAFG-DT` 1 735.17
- `ZNF710-AS1` 1 735.28
+ AF001548.1 1 735.44
+ AL445524.1 1 735.44
- `MIR1-1HG-AS1` 1 735.54
- U62317.1 1 735.76
+ `C5orf66-AS1` 1 735.89
+ `MBNL1-AS1` 1 735.91
+ MIR205HG 1 735.92
+ AP003071.4 1 735.99
- `PCCA-DT` 1 736.08
+ AATBC 1 736.14
+ LINC01133 1 736.25
+ AC099850.4 1 736.27
- AL161431.1 1 736.29
+ AC245041.2 1 736.30
+ AC008735.2 1 736.34
+ AC025575.2 1 736.36
+ SNHG12 1 736.36
+ MIR200CHG 1 736.42
+ LINC02657 1 736.46
- AL023284.4 1 736.47
+ `KRT7-AS` 1 736.49
+ SNHG25 1 736.50
+ `PPP1R14B-AS1` 1 736.52
+ `ADAMTS9-AS1` 1 736.53
+ U62317.2 1 736.53
+ FENDRR 1 736.56
+ `ACTA2-AS1` 1 736.60
+ AC008736.1 1 736.60
+ `HAND2-AS1` 1 736.60
+ `PGM5-AS1` 1 736.61
+ MIR22HG 1 736.61
+ AL390719.2 1 736.66
+ `VPS9D1-AS1` 1 736.67
+ CASC9 1 736.67
+ AC053503.3 1 736.67
+ AC005180.1 1 736.68
+ AL162424.1 1 736.68
+ LINC01082 1 736.71
+ AC079313.2 1 736.71
+ LINC00865 1 736.71
+ AC092718.4 1 736.71
+ AC020916.1 1 736.72
+ AC018904.1 1 736.72
- NR4A1AS 1 736.74
- MIR100HG 1 736.79
- LINC01980 1 736.86
- `SPINT1-AS1` 1 736.89
- AP001107.5 1 737.71
- PCAT6 1 738.00
- AC015912.3 1 738.51
- AC093010.3 1 739.96
- AC005180.2 1 739.97查看下结果：
summary(fit.step)
Call:
coxph(formula = Surv(time_months, event) ~ AC005180.2 + MIR100HG +
AP001107.5 + NR4A1AS + AC093010.3 + PCAT6 + AC036108.3 +
AL049555.1 + `MIR1-1HG-AS1` + `SPINT1-AS1` + LINC01980 +
`ZNF710-AS1` + AL161431.1 + `PCCA-DT` + U62317.1 + AL023284.4 +
AC015912.3, data = dat.cox)
n= 297, number of events= 71
coef exp(coef) se(coef) z Pr(>|z|)
AC005180.2 0.83728 2.31007 0.30399 2.754 0.00588 **
MIR100HG 0.59783 1.81818 0.28762 2.079 0.03766 *
AP001107.5 -1.63414 0.19512 0.80039 -2.042 0.04118 *
NR4A1AS 0.48095 1.61761 0.22753 2.114 0.03453 *
AC093010.3 -0.65126 0.52139 0.24377 -2.672 0.00755 **
PCAT6 0.34371 1.41017 0.14756 2.329 0.01985 *
AC036108.3 1.48666 4.42231 0.96894 1.534 0.12495
AL049555.1 0.23994 1.27117 0.15931 1.506 0.13203
`MIR1-1HG-AS1` -1.95396 0.14171 1.21419 -1.609 0.10756
`SPINT1-AS1` 0.48314 1.62116 0.23919 2.020 0.04340 *
LINC01980 -0.17331 0.84087 0.08564 -2.024 0.04300 *
`ZNF710-AS1` -0.55489 0.57414 0.35644 -1.557 0.11953
AL161431.1 0.17192 1.18758 0.08854 1.942 0.05217 .
`PCCA-DT` -0.35372 0.70207 0.19499 -1.814 0.06966 .
U62317.1 -0.20670 0.81326 0.12168 -1.699 0.08937 .
AL023284.4 -0.30170 0.73956 0.15674 -1.925 0.05425 .
AC015912.3 -0.43266 0.64878 0.18428 -2.348 0.01888 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
AC005180.2 2.3101 0.4329 1.27311 4.1916
MIR100HG 1.8182 0.5500 1.03470 3.1949
AP001107.5 0.1951 5.1251 0.04065 0.9367
NR4A1AS 1.6176 0.6182 1.03563 2.5266
AC093010.3 0.5214 1.9180 0.32334 0.8407
PCAT6 1.4102 0.7091 1.05601 1.8831
AC036108.3 4.4223 0.2261 0.66204 29.5404
AL049555.1 1.2712 0.7867 0.93025 1.7370
`MIR1-1HG-AS1` 0.1417 7.0566 0.01312 1.5308
`SPINT1-AS1` 1.6212 0.6168 1.01443 2.5908
LINC01980 0.8409 1.1892 0.71094 0.9946
`ZNF710-AS1` 0.5741 1.7417 0.28551 1.1546
AL161431.1 1.1876 0.8420 0.99839 1.4126
`PCCA-DT` 0.7021 1.4244 0.47908 1.0288
U62317.1 0.8133 1.2296 0.64070 1.0323
AL023284.4 0.7396 1.3522 0.54395 1.0055
AC015912.3 0.6488 1.5413 0.45211 0.9310
Concordance= 0.735 (se = 0.029 )
Likelihood ratio test= 54.26 on 17 df, p=9e-06
Wald test = 47.59 on 17 df, p=1e-04
Score (logrank) test = 51.16 on 17 df, p=3e-05最终59个变量剩下17个，筛选效果还不错。
这个筛选过程是根据AIC进行的，一般会选择AIC最小的结果。AIC全称赤池信息量准则(Akaike information criterion，AIC)，是评估统计模型的复杂度和衡量统计模型”拟合优度”（Goodness of Fit）的一种标准，是由日本统计学家赤池弘次创立和发展的。赤池信息量准则的方法是寻找可以最好地解释数据但包含最少自由参数的模型。
查看最终的AIC和BIC：
# 初始模型的AIC
AIC(fit.cox)
## [1] 794.1182
# 筛选后的AIC和BIC
AIC(fit.step)
## [1] 734.7167
BIC(fit.step)
## [1] 773.1823查看回归系数：
step.coef <- coef(fit.step)
step.coef
## AC005180.2 MIR100HG AP001107.5 NR4A1AS AC093010.3
## 0.8372759 0.5978335 -1.6341403 0.4809485 -0.6512620
## PCAT6 AC036108.3 AL049555.1 `MIR1-1HG-AS1` `SPINT1-AS1`
## 0.3437133 1.4866631 0.2399385 -1.9539623 0.4831391
## LINC01980 `ZNF710-AS1` AL161431.1 `PCCA-DT` U62317.1
## -0.1733137 -0.5548853 0.1719177 -0.3537250 -0.2067033
## AL023284.4 AC015912.3
## -0.3017031 -0.4326573提取这17个变量的名字：
step.lnc <- names(coef(fit.step))
step.lnc
## [1] "AC005180.2" "MIR100HG" "AP001107.5" "NR4A1AS"
## [5] "AC093010.3" "PCAT6" "AC036108.3" "AL049555.1"
## [9] "`MIR1-1HG-AS1`" "`SPINT1-AS1`" "LINC01980" "`ZNF710-AS1`"
## [13] "AL161431.1" "`PCCA-DT`" "U62317.1" "AL023284.4"
## [17] "AC015912.3"简单。
当然使用broom也是可以提取这个结果的：
broom::tidy(fit.step)
## # A tibble: 17 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 AC005180.2 0.837 0.304 2.75 0.00588
## 2 MIR100HG 0.598 0.288 2.08 0.0377
## 3 AP001107.5 -1.63 0.800 -2.04 0.0412
## 4 NR4A1AS 0.481 0.228 2.11 0.0345
## 5 AC093010.3 -0.651 0.244 -2.67 0.00755
## 6 PCAT6 0.344 0.148 2.33 0.0198
## 7 AC036108.3 1.49 0.969 1.53 0.125
## 8 AL049555.1 0.240 0.159 1.51 0.132
## 9 `MIR1-1HG-AS1` -1.95 1.21 -1.61 0.108
## 10 `SPINT1-AS1` 0.483 0.239 2.02 0.0434
## 11 LINC01980 -0.173 0.0856 -2.02 0.0430
## 12 `ZNF710-AS1` -0.555 0.356 -1.56 0.120
## 13 AL161431.1 0.172 0.0885 1.94 0.0522
## 14 `PCCA-DT` -0.354 0.195 -1.81 0.0697
## 15 U62317.1 -0.207 0.122 -1.70 0.0894
## 16 AL023284.4 -0.302 0.157 -1.92 0.0542
## 17 AC015912.3 -0.433 0.184 -2.35 0.0189一目了然，简洁清晰，broom真的是神包~
还可以1行代码查看模型的各种统计值，包括P值、R2、AIC、BIC、C-index等等：
broom::glance(fit.step)
## # A tibble: 1 × 18
## n nevent statistic.log p.value.log statistic.sc p.value.sc statistic.wald
## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 297 71 54.3 0.00000900 51.2 0.0000279 47.6
## # ℹ 11 more variables: p.value.wald <dbl>, statistic.robust <dbl>,
## # p.value.robust <dbl>, r.squared <dbl>, r.squared.max <dbl>,
## # concordance <dbl>, std.error.concordance <dbl>, logLik <dbl>, AIC <dbl>,
## # BIC <dbl>, nobs <int>16.4 自助法stepwise
这里再给大家介绍下一种自助法stepwise，可以通过自助法重抽样进行逐步筛选变量，比如进行1000次bootstrap。
该方法借助bootStepAIC实现，是基于stepAIC()函数的，支持”lm”, “aov”,“glm”, “negbin”, “polr”, “survreg”,以及”coxph”。
使用方法也很简单，下面是一个10次bootstrap的逐步选择法(因为耗时太长了，我只用了10次) ：
library(bootStepAIC)
# 10次bootstrap
fit.boot <- boot.stepAIC(fit.cox,data=dat.cox,direction="both",B=10,seed=123)
#fit.boot
Final Model:
Surv(time_months, event) ~ AC005180.2 + MIR100HG + AP001107.5 +
NR4A1AS + AC093010.3 + PCAT6 + AC036108.3 + AL049555.1 +
`MIR1-1HG-AS1` + `SPINT1-AS1` + LINC01980 + `ZNF710-AS1` +
AL161431.1 + `PCCA-DT` + U62317.1 + AL023284.4 + AC015912.3提取变量名字：
names(coef(fit.boot$OrigStepAIC))
[1] "AC005180.2" "MIR100HG" "AP001107.5" "NR4A1AS"
[5] "AC093010.3" "PCAT6" "AC036108.3" "AL049555.1"
[9] "`MIR1-1HG-AS1`" "`SPINT1-AS1`" "LINC01980" "`ZNF710-AS1`"
[13] "AL161431.1" "`PCCA-DT`" "U62317.1" "AL023284.4"
[17] "AC015912.3" 和不进行bootstrap的方法得到的结果是一样的：
identical(sort(step.lnc),sort(names(coef(fit.boot$OrigStepAIC))))
[1] TRUE