目录
最优子集回归是从一个较多自变量的模型中，选出优秀的自变量集合做模型的一个方法。
什么是最优子集回归?
- 当我们进行回归分析时，我们可能会碰到自变量过多, 且存在多重共线性, 或者全模型可能是过度拟合的.
- 针对这种情况，可以人为根据经验判断筛选对因变量有影响的自变量，比如距离周边学校的距离去预测房价.
- 但通常我们并不是相关领域的专家, 对可能影响因变量的自变量并不了解，于是我们需要运用算法获得预测效果最优的模型, 进而接近真实模型，比如最优子集回归.
- 最优子集回归, 即对$p$个预测变量的所有可能组合分别使用最小二乘回归进行拟合
- 对含1变量的模型，拟合$p$个模型; 对含两个变量的模型，拟合$p(p-1)/2$个模型，以此类推，总共拟合$2^p$个模型.
- 按照一定的比较准则（如AIC）从中选择一个最优模型.
如何获得全部模型?
- 最优子集回归需要建立大量的模型，手动输入显然不可行，可以使用二进制向量来表示一个模型，值为1表示选择该变量，值为0表示不选择该变量.
- 比如总共有3个自变量时，$(0,1,1)$表示使用第2和第3个自变量建立模型.
- 将所有可能的向量按行排列，得到$2^p\times p$维的矩阵$Z$.
- 矩阵$Z$的行转化为十进制即为$[0,2^p-1]$的所有整数，所以，$[0,2^p-1]$中的所有整数转化为二进制向量并按行排列即可表示所有可能的模型.
循环法获得矩阵$Z$
- 解方程$(1)$，得到$(k_0,k_1,…,k_{p-1})$, 即可将十进制的$n$转化为二进制$p$维向量。
- turnbits_cir 函数即可实现此功能。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
| ##turnbits_cir函数用于将十进制转化为p维二进制向量
##参数n为十进制数，p为二进制向量维数
##输出p维逻辑行向量
turnbits_cir = function(n,p){
z = rep(0,p) #预留空间
tn = n #tn为十进制数值
z[1] = tn%%2 #tn/2的余数为z的第一个值
for(j in 2:p){
tn = (tn-z[j-1])/2 #更新tn的值
if(tn == 0) break #tn=0时跳出循环
z[j] = tn%%2 #z的第j个值为tn/2的余数
}
as.logical(z) #将向量z转化为逻辑值
}
##以p=3为例
p = 3
Z = matrix(unlist(lapply(0:(2^p-1) , turnbits_cir , p)) , #将0：2^p-1转化为p维二进制向量，并按行排列
ncol = p ,byrow = T)
Z
|
这个矩阵的作用是穷举全部自变量的子集。
递归法获得矩阵$Z$
1
2
3
4
5
6
7
8
9
10
11
12
13
14
| ##turnbits_rec函数用于将[0,2^p-1]中的整数转化为(2^p)*p维逻辑矩阵
##输出(2^p)*p维逻辑矩阵
turnbits_rec=function(p){
if(p==1) return (matrix(c(FALSE,TRUE),ncol=1)) #p=1时返回初始值
else return
(cbind(rbind(turnbits_rec(p-1),turnbits_rec(p-1)), #p≠1时给第p个矩阵左侧分别添加一列0和1
#再按行拼接起来得到第p+1个矩阵
rep(c(FALSE,TRUE),rep(2^(p-1),2))))
}
Z=turnbits_rec(3)
Z
|
R语言实现最优子集回归
以AIC为比较准测使用 bestlm 函数实现最优子集回归。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
| ##bestlm函数用于进行最优子集回归
##参数z为二进制行向量，代表不同的模型，参数dataTR为训练集，dataTE为测试集
##dataTR和dataTE应有相同的结构，且因变量在最后一列
##输出2*2^p维矩阵，第一行为测试误差，第二行为AIC值，不同列代表不同模型
bestlm = function(z, dataTR,dataTE){
p = dim(dataTR)[2] #求数据集共有多少列
yxn = names(dataTR) #向量yxn为所有变量的名称
yn = yxn[p] #定义yn为因变量名称
xn = yxn[1:p-1] #定义向量xn为所有自变量名称
tm1 = paste(yn,"~",sep="") #定义tm1为yn~
{if(sum(z) == 0) tm2 = "1" #如果自变量个数为0，tm2为1
else tm2 = paste(xn[z] , collapse = "+")} #自变量个数不为0，tm2为各自变量名称相加
fam = formula(paste(tm1, tm2, sep = "")) #将tm1和tm2拼接并转化为公式
lm1 = lm(fam , dataTR) #建立回归模型
TE = mean((dataTE[,p] - predict(lm1,dataTE))^2) #计算测试误差
AIC = extractAIC(lm(fam,dataTR))[2] #计算模型AIC值
c(TE,AIC) #拼接测试误差和AIC值
}
|
例子
- 接下来，我们用 ElemStatLearn 包中的 prostate 数据集进行最优子集回归
- prostate 是关于前列腺切除手术的数据集，共包含十个变量，其中前八列为自变量，第九列为因变量，第十列是逻辑变量，用于区分训练集与测试集。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
| library(ElemStatLearn) #加载ElemStatLearn包
data(prostate) #加载prostate数据集
data = prostate
p = dim(data)[2]-2 #计算自变量个数
datatr = data[data[,p+2],1:(p+1)] #获得训练集
datate = data[!data[,p+2],1:(p+1)] #获得测试集
Z = turnbits_rec(p) #计算二进制矩阵
mAIC = apply(Z,1,bestlm,datatr,datate) #进行最优子集回归
names(data)[1:p][Z[which.min(mAIC[1,]),]] #按测试误差最小选出最优模型
names(data)[1:p][Z[which.min(mAIC[2,]),]] #按AIC最小选择最优模型
#min(mAIC[2,])
plot(mAIC[1,], mAIC[2,])
|
1
2
3
4
| head(data)
lm1 <- lm(lpsa~.,data[,-10])
extractAIC(lm1)
AIC(lm1)
|
和以前一样，结果省略不贴了，大家自己跑一下。
最小角回归
这里过程比较抽象。
最小角回归是一种变量选择方法，它大大减少了逐步回归过程中的迭代次数，整个求解过程可以被控制在m步内完成.
在传统的向前逐步回归中，我们将每一个估计系数初始值设为0，找到和响应变量相关性最大的自变量，然后尽可能地在这个方向上走最远的距离，直到和当前残差值有很大相关性的另一个自变量的出现。最小角回归基于这种思想，它每一次沿着当前所有模型中变量的角平分线走，直到另一个变量的相关性大到足够进入模型.
现在具体阐述最小角回归的几何原理和迭代步骤：
我们的目标是找一个$\hat{\beta}$用于拟合,拟合值$\hat{\mu}=X\hat{\beta},$使得平方误差最小且满足条件$\sum_{j=1}^{m} \hat{\beta_j}<t.$
首先,引入一些符号表示,设$X$中的变量都是独立的,$X = (x_1,x_2,\cdots,x_m),A$是${1,2,\cdots,m}$的子集,用于记录模型中已有的变量,一开始是空的,$X_A=(\cdots s_jx_j \cdots)_{j\in A},s_j=\pm1,$其中$s_j$由该变量与残差的相关性决定,$G_A =X_A^TX_A,A_A=({1^T}G_A^{-1}{1})^{-1/2},W_A=A_AG_A^{-1}1,$单位角平分线向量$U_A=X_AW_A.$
简单给出角平分线向量的求解过程,由单位角平分线的性质$X_A^TU_A=A_A1,||U_A||^2=1,$求解$W_A$和$A_A,X_A^TX_AW_A=A_A1,W_A=(X_A^TX_A)^{-1}A_A1,$代入$U_A^TU_A=1$得$A_A^21^T(X_A^TX_A)^{-1}(X_A^TX_A)(X_A^TX_A)^{-1}1=1,$故$A_A=({1^T}G_A^{-1}{1})^{-1/2}.$
算法完整迭代步骤,在计算之前，对$X$进行中心化标准化,对$y$进行中心化操作,拟合值$\mu_A$初始值为$0,$估计系数$\hat{\beta}$初始值为$0,$计算所有变量与当前残差的相关性,即$\hat C=X^Ty,$每一个分量分别代表每一个变量与残差的相关性,找出其中最大的$C_j,$记为$C_{max},$将所有相关性等于$C_{max}$的变量编号放进集合A中，其对应的符号记录在集合$S$中,然后计算当前所有集合A中变量的角平分线,即$X_A$列向量的角平分线,这里需要注意的是之所以引入符号$s_j$是为了是列向量经过符号调整后与残差方向的夹角小于$90^。,$这样才能保证我们前进的方向是正确的,才能使拟合值不断靠近y而不是偏离y,计算角平分线到每一个变量的投影,由于数据经过标准化,也即夹角余弦值,当前角平分线的方向决定我们即将要前进的方向,而要走多远由$\hat\gamma$决定,$\hat\gamma$是使得走完这一步以后更新的残差能够与加入系新的变量以后重新计算出的角平分线平行的数值,即现有A集合中的变量以及加入的新变量与残差的相关性相同,$\hat\gamma=min_{j\in A^c}^+(\frac{C_max-C_j}{A_A-a_j},\frac{C_max+C_j}{A_A+a_j}),$可以看出$\hat\gamma$不仅告诉了我们这一步要走多远,还告诉了我们下一步加入A集的变量有哪些.依次更新$\mu_A,\hat{\beta},\hat C$,其中$\mu_A=\mu_A+\hat\gamma U_A,\hat{\beta} = \hat{\beta}+\hat\gamma \beta^+,\beta^+$为$W_A$填充$0$使维度与$\hat{\beta}$相同,$\hat C = \hat C - \hat\gamma a,$如此迭代下去,m步之内完成求解.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
| mchol <- function(x)
{
mn <- dim(x)
m <- mn[1]
n <- mn[2]
#检验传入参数是否符合要求
if(m != n)
{
return ("Wrong dimensions of matrix!")
}
if(sum(t(x) != x) > 0)
{
return ("Input matrix is not symmetrical!")
}
L <- matrix(0, m, m)
for(i in 1:m) #以列为单位求解L
{
L[i,i] <- sqrt(x[i,i])
if(i < m)
{
L[(i+1):m,i] <- x[(i+1):m,i]/L[i,i]
#更新矩阵，便于下一次同样方法计算
TLV <- L[(i+1):m,i]
TLM <- matrix(TLV, m-i, m-i)
TLM <- sweep(TLM, 2, TLV, "*")
x[(i+1):m,(i+1):m] <- x[(i+1):m,(i+1):m] - TLM
}
}
L
}
mforwardsolve=function(L,b){
mn=dim(L); m=mn[1]; n=mn[2]
if(m!=n) return ("Wrong dimensions of matrix L!")
if(m!=length(b)) return ("Wrong dimensions of matrix L or vector b!") #检查输入参数是否符合要求
x=rep(0,m)
for(i in 1:m){
x[i]=b[i]/L[i,i]
if(i<m) b[(i+1):m]=b[(i+1):m]-x[i]*L[(i+1):m,i] #更新右端项，逐步减去已求出的未知量对右端项的贡献
}
x
}
mbacksolve=function(L,b){
mn=dim(L); m=mn[1]; n=mn[2]
if(m!=n) return ("Wrong dimensions of matrix L!")
if(m!=length(b)) return ("Wrong dimensions of matrix L or vector b!") #检查输入参数是否符合要求
x=rep(0,m)
for(i in m:1){
x[i]=b[i]/L[i,i]
if(i>1) b[(i-1):1]=b[(i-1):1]-x[i]*L[(i-1):1,i] #更新右端项
}
x
}
######################################################以上为一些需要用到的函数
mylars = function(x,y,t){
eps = 1e-8
xnames = names(x)
x = as.matrix(cbind(x))
nm = dim(x)
n = nm[1]
m = nm[2]
one = rep(1, n)#中心化
meanx = drop(one %*% x)/n
x = scale(x, meanx, FALSE)
mu = mean(y)
y = drop(y - mu)
normx = sqrt(drop(one %*% (x^2)))#标准化
x = scale(x, FALSE, normx)
xtx = t(x)%*%x
muA = rep(0,n)#初始化拟合值
betahat = rep(0,m)#初始化估计系数
Chat = t(x)%*%y#初始化变量和残差的相关性,因为当前拟合值为0,所以残差就是y
A = NULL
S = NULL
k = 1
while(sum(abs(betahat))< t&&k <= m){
Cmax = max(Chat)
SS = sign(Chat)
if(is.null(A)){
joinA = Chat >= Cmax - eps
A = (1:m)[joinA]
S = SS[joinA]
LA = mchol(t(scale(t(scale(xtx[A,A], FALSE, S)), FALSE, S)))
}
else{
joinA[A] = FALSE
joinA[-A] = Chat[-A] >= Cmax - eps#已经在模型里的变量不再重复加入
for(i in (1:m)[joinA]){#利用原有的矩阵更新cholesky分解矩阵
sign_i = SS[i]
Lvec = mforwardsolve(LA, xtx[A, i]*sign_i)
Lvec_ = sqrt(xtx[i,i] - sum(Lvec*Lvec))
LA = as.matrix(rbind(cbind(LA,0),c(Lvec,Lvec_)))
A = c(A, i)#记录加入模型的变量编号
S = c(S, sign_i)#记录对应变量的符号
}
GA1 = mbacksolve(t(LA), mforwardsolve(LA, rep(1, length(A))))#求解GA1
AA = 1/sqrt(sum(GA1))
WA = AA*GA1
uA = scale(x[,A], FALSE, S)%*%WA#当前角平分线
a = scale(xtx[,A], FALSE, S)%*%WA#各个分量在角平分线的投影
gammahat = min(ifelse(Chat[-A]>0,(Cmax-Chat[-A])/(AA-a[-A]),(Cmax+Chat[-A])/(AA+a[-A])))#找出最优的步长
muA = muA + gammahat*uA#更新拟合值
betaplus = rep(0, m)
betaplus[A] = WA
betahat = betahat + gammahat*betaplus#更新系数估计
Chat = Chat - gammahat*a#更新变量与残差的相关性
k = k+1
}
}
residuals = y - muA
names(betahat) = xnames
betahat_ = scale(t(betahat), FALSE, normx)#还原估计系数
RSS = apply(residuals^2, 2, sum)#计算残差平方和
cat('RSS:', RSS,' ')
cat('mu:',mu)
betahat_
}
library(ElemStatLearn)
data(prostate)
x = prostate[,1:8]
y = prostate[,9]
mylars(x, y, 1000)
|
本文链接：
https://konelane.github.io/2019/03/26/190326数据挖掘最优子集回归&最小角回归/
-- EOF --
转载请注明出处 署名-非商业性使用-禁止演绎 3.0 国际（CC BY-NC-ND 3.0）
