一些变量筛选方法——4、模拟实验

本系列博客聚焦于变量筛选的方法，所以前文中提及PCR与PLSR由于只能使数据进行降维，而不能进行变量选择，所以下面的模拟不使用这两种方法。

---

模拟实验

为了比较算法的优劣，我们构造了下列模拟（前两个是线性情况，后三个是非线性情况），其中部分参考现有文献中的构造，包含了线性和非线性的情况：

$\text{1)}Y=c_1{\beta }_1{X}_1+c_2{\beta }_2{X}_2+c_3{\beta }_3\text{I}(X_3<0)+c_4{\beta }_4{X}_4+ϵ,$
$\text{2)}Y=c_1{\beta }_1{X}_1{X}_2+c_3{\beta }_2\text{I}(X_3<0)+c_4{\beta }_3{X}_4+ϵ,$
$\text{3)}Y=c_1{\beta }_1{X}_1+c_2{\beta }_2{X}_2+c_3{\beta }_3\text{I}(X_3<0)+\exp (c_4|X_4|)+ϵ,$
$\text{4)}Y=b_1|X_1|+b_2\frac{X_2}{2-X_2}+b_3\frac{\exp (X_3)}{3}+ϵ,$
$\text{5)}Y=b_1\sin (X_1)+b_2\tanh (X_2)+b_3\frac{X_1}{X_2}+ϵ,$

其中，$I(X_3<0)$为示性函数，当$X_3<0$时为1，否则为0；${\beta }_j=(-1{)}^U(a+|Z|),j=1,2,3,4$。$a=4\log n/\sqrt{n}\text{，}U\sim \text{Bernoulli}(0.4)$，$b_k\sim \text{Uniform}(-1,1)$， $k=1,2,3$ 以及$Z\sim N(0,1)$。同时设置$(c_1,c_2,c_3,c_4)=(2,0.5,3,2)$。

另外，$X_i\sim N(0,1),i=1,2,\cdots ,p$，$ϵ\sim N(0,1)$。

由于最优子集选择等方法适用于维度不是非常高的情况，也为了算法能在交论文之前提交出结果。所以我们选择$n=200$，分别进行模拟。

---

模拟实验一：高维情况

下面取数据维数$p=10,20,30$，分别使用：最优子集选择，向前逐步回归，向后逐步回归（前三种方法我们都使用AIC指标），LASSO（使用CV，最优$\lambda$的选择我们选择1倍标准误的$\lambda$，即${\lambda }_{1se}$，其效果要由于${\lambda }_{min}$，而由于选用了${\lambda }_{1se}$，会导致标准较为严苛），SCAD，PDAS。由于基于Pearson相关系数的SIS，基于距离相关性的SIS，随机森林方法来进行变量选择，这三种方法只能人为指定选择出来的变量个数所以在这个模拟实验中不进行比较。

模拟重复100次，记录选出相关变量的个数平均值与方差，以及选到但选错变量的个数平均值与方差。例如：1）式表示$X_1,X_2,X_3,X_4$ 被选出的个数，以及所有选出的变量中选错的个数。单独写出，便得到如下定义：

• True Positive = 选出的所有变量中，正确的相关变量个数
• False Positive = 选出的所有变量中，错误的相关变量个数
• False Positive + True Positive = 选出的所有变量

公式一、二、三中，相关变量个数均为4个；公式四：3个；公式五：2个。下面展示模拟结果。（注：由于写课程报告的时间实在有限，模拟维数一旦过高，最优子集选择、向前逐步回归、向后逐步回归三种方法的速度会呈指数上升，所以本模拟只使用$p=10,20,30$三种情况）

方法	公式1	公式2	公式3	公式4	公式5
Best SS	4 ± 0	2.36 ± 0.89	2.04 ± 0.9	1.45 ± 0.81	0.73 ± 0.74
Forward	4 ± 0	2.9 ± 0.75	2.27 ± 1.13	1.7 ± 0.92	0.72 ± 0.73
Backward	4 ± 0	2.9 ± 0.76	2.26 ± 1.12	1.69 ± 0.92	0.72 ± 0.73
LASSO	4 ± 0	1.07 ± 0.41	0.35 ± 0.54	0.29 ± 0.46	0.26 ± 0.54
SCAD	4 ± 0	1.62 ± 0.94	1.95 ± 1.1	1.55 ± 0.69	1.28 ± 0.45
PDAS	4 ± 0	2.36 ± 0.89	2.04 ± 0.9	1.44 ± 0.81	0.73 ± 0.74

Table 1: True Positive p = 10

方法	公式1	公式2	公式3	公式4	公式5
Best SS	4 ± 0	2.51 ± 0.87	2.09 ± 0.99	1.37 ± 0.84	0.75 ± 0.78
Forward	4 ± 0	2.5 ± 0.88	1.81 ± 1.04	1.22 ± 0.84	0.61 ± 0.72
Backward	4 ± 0	2.48 ± 0.9	1.82 ± 1.05	1.21 ± 0.84	0.6 ± 0.72
LASSO	3.99 ± 0.1	1 ± 0.2	0.29 ± 0.46	0.26 ± 0.44	0.11 ± 0.42
SCAD	4 ± 0	1.58 ± 0.85	1.67 ± 0.91	1.56 ± 0.69	1.26 ± 0.44
PDAS	4 ± 0	2.5 ± 0.87	2.09 ± 0.99	1.37 ± 0.84	0.75 ± 0.78

Table 2: True Positive p = 20

方法	公式1	公式2	公式3	公式4	公式5
Best SS	4 ± 0	2.48 ± 0.87	2.16 ± 0.91	1.48 ± 0.83	0.71 ± 0.76
Forward	3.99 ± 0.1	2.18 ± 0.82	1.79 ± 1.03	1.22 ± 0.84	0.52 ± 0.7
Backward	3.99 ± 0.1	2.2 ± 0.83	1.77 ± 1.01	1.2 ± 0.8	0.51 ± 0.69
LASSO	3.99 ± 0.1	1 ± 0.14	0.23 ± 0.42	0.19 ± 0.39	0.18 ± 0.48
SCAD	3.99 ± 0.1	1.54 ± 0.81	1.66 ± 0.88	1.52 ± 0.61	1.26 ± 0.44
PDAS	4 ± 0	2.46 ± 0.85	2.14 ± 0.92	1.46 ± 0.85	0.71 ± 0.76

Table 3: True Positive p = 30

方法	公式1	公式2	公式3	公式4	公式5
Best SS	0.94 ± 0.96	1.2 ± 0.98	1.01 ± 0.89	1.27 ± 1.03	1.52 ± 0.93
Forward	0.55 ± 1.16	2.21 ± 1.55	1.59 ± 1.41	1.72 ± 1.6	1.23 ± 1.46
Backward	0.55 ± 1.16	2.21 ± 1.55	1.6 ± 1.41	1.73 ± 1.6	1.23 ± 1.46
LASSO	0.45 ± 0.74	0 ± 0	0.04 ± 0.2	0.07 ± 0.26	0.12 ± 0.33
SCAD	1.67 ± 1.32	1.64 ± 1.37	1.16 ± 1.73	1.07 ± 1.7	0.93 ± 1.67
PDAS	0.94 ± 0.96	1.2 ± 0.98	1.01 ± 0.89	1.25 ± 1.02	1.52 ± 0.93

Table 4: False Positive p = 10

方法	公式1	公式2	公式3	公式4	公式5
Best SS	2.64 ± 1.65	3.06 ± 1.45	2.88 ± 1.47	3.05 ± 1.68	3.28 ± 1.74
Forward	1.08 ± 2.45	2.97 ± 2.21	2.53 ± 2.88	2.05 ± 2.16	1.31 ± 1.91
Backward	1.08 ± 2.45	2.99 ± 2.21	2.52 ± 2.88	2.06 ± 2.15	1.32 ± 1.94
LASSO	0.78 ± 1.1	0 ± 0	0.07 ± 0.26	0.15 ± 0.36	0.16 ± 0.37
SCAD	2 ± 1.6	1.95 ± 2	1.17 ± 1.92	1.57 ± 2.49	1.07 ± 1.89
PDAS	2.63 ± 1.65	3.05 ± 1.44	2.88 ± 1.45	3.06 ± 1.7	3.28 ± 1.74

Table 5: False Positive p = 20

方法	公式1	公式2	公式3	公式4	公式5
Best SS	4.55 ± 2.29	4.99 ± 1.76	4.67 ± 1.93	5.38 ± 2.43	5.16 ± 2.19
Forward	1.31 ± 3.22	3.26 ± 2.35	3.14 ± 4.15	2.47 ± 2.71	1.76 ± 3.64
Backward	1.31 ± 3.22	3.24 ± 2.29	3.16 ± 4.11	2.49 ± 2.72	1.77 ± 3.63
LASSO	1.16 ± 1.77	0 ± 0	0.15 ± 0.36	0.14 ± 0.35	0.23 ± 0.42
SCAD	2.8 ± 2.37	2.31 ± 2.61	1.53 ± 2.2	1.76 ± 2.44	1.18 ± 2.01
PDAS	4.59 ± 2.34	5.03 ± 1.78	4.66 ± 1.88	5.35 ± 2.36	5.18 ± 2.25

Table 6: False Positive p = 30
从上面的结果中可以看出，最优子集选择的结果和PDAS的方法是几乎完全一样的，但前者根本跑不了50维以上的数据，这个从算法理论上也可以看出。这两种方法虽然选对的很多，但是选错的也相较于其他方法会更多，所以可以理解为是比较“激进”的方法。而逐步回归法在维度低的时候效果还不错，但是一旦维度有所提高，效果就渐渐下降。而LASSO可以发现，它是一种比较保守的方法，选择的总变量都比其他方法要少很多，最后选出来的变量稀疏性非常强（也因为我们通过CV选择的$\lambda$是$\lambda_{1se}$）。而最后SCAD效果就要明显好于LASSO，基本上算是$L^0$和$L^1$的权衡。

下面来比较一下六种方法的耗时：

方法	公式1	公式2	公式3	公式4	公式5
Best SS	0.0007	0.0019	0.0013	0.0017	0.0009
Forward	0.0009	0.0012	0.0007	0.0011	0
Backward	0.0003	0.0006	0.0015	0.0011	0.0005
LASSO	0.0418	0.0458	0.0429	0.0473	0.0432
SCAD	0.0322	0.0323	0.035	0.0321	0.0329
PDAS	0.1238	0.1243	0.1263	0.1306	0.1262

Table 7: 方法耗时，p = 10

方法	公式1	公式2	公式3	公式4	公式5
Best SS	0.0006	0.003	0.0023	0.0033	0.003
Forward	0.0002	0.0004	0.0015	0.0004	0.0002
Backward	0.0016	0.0002	0.0006	0.0006	0
LASSO	0.0444	0.0451	0.0466	0.0441	0.0444
SCAD	0.0421	0.0432	0.0441	0.044	0.0431
PDAS	0.1357	0.135	0.1352	0.1352	0.1351

Table 8: 方法耗时，p = 20

方法	公式1	公式2	公式3	公式4	公式5
Best SS	0.0699	0.2519	0.238	0.3365	0.3291
Forward	0.0008	0.0007	0.0005	0.0001	0.0012
Backward	0.0005	0.0007	0.0012	0.0007	0.0002
LASSO	0.0464	0.0489	0.0481	0.0478	0.0481
SCAD	0.0531	0.0572	0.0594	0.0611	0.0608
PDAS	0.1459	0.1459	0.1479	0.1461	0.1461

Table 9: 方法耗时，p = 30

其实从这个耗时中，只能初步发现最优子集选择和逐步回归，它们在维数低的时候最快，但当维数提升了之后，速度就大大下降（由于时间原因，应该还要做$p=50,100$才更有比较价值）。而LASSO和SCAD两种方法速度相当，PDAS在高维时候，速度要明显快于最优子集选择。

所以综上可以发现，如果要做最优子集选择，我们可以直接使用PDAS方法，然后面对维数较高的情形时，或者$p>n$的时候，我们可以使LASSO，SCAD与PDAS方法。这里值得一提的是，LASSO在模拟的表现中相对不突出，但在平时做项目的过程中发现，针对实际数据，其实LASSO的效果也非常突出。

---

模拟实验二：超高维情况

下面取数据维数$p=100,200$，由于最优子集选择，向前逐步回归，向后逐步回归做不动超高维的情况，LASSO、SDAR与PDAS，用下述模拟方法无法进行比较。所以这里使用：基于Pearson相关系数的SIS，基于距离相关性的SIS，随机森林方法三种方法来进行变量选择。

模拟重复100次，记录每次将所有相关变量全部选出，相关变量重要性（相关性）排最后的在第几位。例如：1）式表示$X_1,X_2,X_3,X_4$全部被选出，假设$X_1$重要性（相关性）排在第3位置；$X_2$排第1；$X_2$排第1；$X_3$排第15；$X_4$排第10；$X_5$排第9，则此时记录15，这个值越小，表示模型越出色。

公式	Cor-SIS	DC-SIS	RF
公式1	63.82 ± 27.67	37.3 ± 31.81	36.55 ± 31.43
公式2	75.07 ± 18.62	70.48 ± 22.12	55.6 ± 27.41
公式3	75.85 ± 20.4	52.75 ± 28.42	74.68 ± 19.97
公式4	71.19 ± 23.66	38.37 ± 35.82	65.89 ± 25.22
公式5	59.02 ± 25.7	14.7 ± 27.02	57.42 ± 27.1

Table 10: 全部相关变量被选出所需最少个数，p = 100

公式	Cor-SIS	DC-SIS	RF
公式1	140.64 ± 46.97	86.95 ± 68.71	85.57 ± 67.67
公式2	150.41 ± 37.84	140.68 ± 45.46	116.96 ± 60.33
公式3	148.57 ± 39.09	104.16 ± 56.48	151.01 ± 39.22
公式4	140.25 ± 45.72	71.08 ± 63.78	136.04 ± 48.36
公式5	125.97 ± 53.4	37.4 ± 60.39	118.68 ± 60.96

Table 11: 全部相关变量被选出所需最少个数，p = 200

通过上表发现，复杂非线性的情况，DC-SIS的表现要明显好于其他两种方法，这是由于DC统计量的性质决定的。而在理论模型相对简单的情况下，RF的效果也表现的不错。

但光看效果，我们不能对一个算法做出评价，还需要看其运算耗时，其耗时如下表所示：

公式	公式1	公式2	公式3	公式4	公式5
Cor-SIS	0.0026	0.0031	0.0019	0.0036	0.0033
DC-SIS	0.2613	0.2616	0.2629	0.2608	0.2635
RF	0.3556	0.3635	0.3951	0.3992	0.4097

Table 12: 三种方法耗时，p = 100

公式	公式1	公式2	公式3	公式4	公式5
Cor-SIS	0.0046	0.0051	0.0073	0.0044	0.0074
DC-SIS	0.5624	0.5606	0.5651	0.5635	0.5606
RF	0.4652	0.4782	0.5382	0.5354	0.558

Table 13: 三种方法耗时，p = 200

通过比较可以发现，最简单的SIS的所用时是最短的，但其实DC-SIS处理数据虽然比较长，也在能接受的范围内。相较于前面的实际选取变量的效果优势，其实是可以适当牺牲一些运算性能。而DC-SIS类似算法都是进行变量初步筛选，所以不能漏掉相关变量，这点非常重要！所以后面针对超高维的数据情况，我们会先用DC-SIS先进行变量初步筛选，将维数降下来，然后再使用LASSO进行回归。

---

后面我们会进行真实案例的分析。传送门：一些变量筛选方法——5、真实数据与总结