變量選擇之SCAD算法

本文提出了一種用于同時(shí)達(dá)到選擇變量和預(yù)測(cè)模型系數(shù)的目的的方法——SCAD。這種方法的罰函數(shù)是對(duì)稱(chēng)且非凹的，并且可處理奇異陣以產(chǎn)生稀疏解。此外，本文提出了一種算法用于優(yōu)化對(duì)應(yīng)的帶懲罰項(xiàng)的似然函數(shù)。這種方法具有廣泛的適用性，可以應(yīng)用于廣義線性模型，強(qiáng)健的回歸模型。借助于波和樣條，還可用于非參數(shù)模型。更進(jìn)一步地，本文證明該方法具有Oracle性質(zhì)。模擬的結(jié)果顯示該方法相比主流的變量選擇模型具有優(yōu)勢(shì)。并且，模型的預(yù)測(cè)誤差公式顯示，該方法實(shí)用性較強(qiáng)。

SCAD的理論理解

在總結(jié)了現(xiàn)有模型的一些缺點(diǎn)之后，本文提出構(gòu)造罰函數(shù)的一些目標(biāo)：

罰函數(shù)是奇異的(singular)

連續(xù)地壓縮系數(shù)

對(duì)較大的系數(shù)產(chǎn)生無(wú)偏的估計(jì)

SCAD模型的Oracle性質(zhì)，使得它的預(yù)測(cè)效果跟真實(shí)模型別無(wú)二致。

并且，這種方法可以應(yīng)用于高維非參數(shù)建模。

SCAD的目標(biāo)函數(shù)如下：

SCAD的罰函數(shù)與$theta$的(近似)關(guān)系如下圖所示。

可見(jiàn)，罰函數(shù)可以用二階泰勒展開(kāi)逼近。

Hard Penality,lasso,SCAD的系數(shù)壓縮情況VS系數(shù)真實(shí)值的情況如下圖所示。

可以看到，lasso壓縮系數(shù)是始終有偏的，Hard penality是無(wú)偏的，但壓縮系數(shù)不連續(xù)。而SCAD既能連續(xù)的壓縮系數(shù)，也能在較大的系數(shù)取得漸近無(wú)偏的估計(jì)。

這使得SCAD具有Oracle性質(zhì)。

SCAD的缺點(diǎn)

模型形式過(guò)于復(fù)雜

迭代算法運(yùn)行速度較慢

在low noise level的情況下表現(xiàn)較優(yōu)，但在high noise level的情況下表現(xiàn)較差。

SCAD的實(shí)現(xiàn)

SCAD迭代公式

SCAD的目標(biāo)函數(shù)如下：

時(shí)，罰函數(shù)可以用二階泰勒展開(kāi)逼近。

從而，有如下迭代公式：

根據(jù)以上公式，代入迭代步驟，即可實(shí)現(xiàn)算法。

SCAD的R實(shí)現(xiàn)

##------數(shù)據(jù)模擬--------

library(MASS)
##mvrnorm()

##定義一個(gè)產(chǎn)生多元正態(tài)分布的隨機(jī)向量協(xié)方差矩陣
Simu_Multi_Norm<-function(x_len, sd = 1, pho = 0.5){
  #初始化協(xié)方差矩陣
  V <- matrix(data = NA, nrow = x_len, ncol = x_len)
 
  #mean及sd分別為隨機(jī)向量x的均值和方差
 
  #對(duì)協(xié)方差矩陣進(jìn)行賦值pho(i,j) = pho^|i-j|
  for(i in 1:x_len){ ##遍歷每一行
    for(j in 1:x_len){ ##遍歷每一列
      V[i,j] <- pho^abs(i-j)
    }
  }
 
  V<-(sd^2) * V
  return(V)
}


##產(chǎn)生模擬數(shù)值自變量X
set.seed(123)
X<-mvrnorm(n = 200, mu = rep(0,10), Simu_Multi_Norm(x_len = 10,sd  = 1, pho = 0.5))

##產(chǎn)生模擬數(shù)值：響應(yīng)變量y
beta<-c(1,2,0,0,3,0,0,0,-2,0)

#alpha<-0

#prob<-exp(alpha + X %*% beta)/(1+exp(alpha + X %*% beta))

prob<-exp( X %*% beta)/(1+exp( X %*% beta))

y<-rbinom(n = 200, size = 1,p = prob)

##產(chǎn)生model matrix
mydata<-data.frame(X = X, y = y)

#X<-model.matrix(y~., data = mydata)

##包含截矩項(xiàng)的系數(shù)
#b_real<-c(alpha,beta)
b_real<-beta

########----定義懲罰項(xiàng)相關(guān)的函數(shù)-----------------


##定義懲罰項(xiàng)

####運(yùn)行發(fā)現(xiàn)，若lambda設(shè)置為2，則系數(shù)全被壓縮為0.

####本程序根據(jù)rcvreg用CV選出來(lái)的lambda設(shè)置一個(gè)較為合理的lambda。


p_lambda<-function(theta,lambda = 0.025){
  p_lambda<-sapply(theta, function(x){
    if(abs(x)< lambda){
      return(lambda^2 - (abs(x) - lambda)^2)
    }else{
      return(lambda^2)
    }
  }
 
  )
 
  return(p_lambda)
}


##定義懲罰項(xiàng)導(dǎo)數(shù)
p_lambda_d<-function(theta,a = 3.7,lambda = 0.025){
  if(abs(theta) > lambda){
    if(a * lambda > theta){
      return((a * lambda - theta)/(a - 1))
    }else{
      return(0)
    }
  }else{
    return(lambda)
  }
}

# ##當(dāng)beta_j0不等于0，定義懲罰項(xiàng)導(dǎo)數(shù)近似
# p_lambda_d_apro<-function(beta_j0,beta_j,a = 3.7, lambda = 2){
#   return(beta_j * p_lambda_d(beta = beta_j0,a = a, lambda = lambda)/abs(beta_j0))
# }
#
#
# ##當(dāng)beta_j0 不等于0，指定近似懲罰項(xiàng),使用泰勒展開(kāi)逼近
# p_lambda_apro<-function(beta_j0,beta_j,a = 3.7, lambda = 2){
#   if(abs(beta_j0)< 1e-16){
#     return(0)
#   }else{
#     p_lambda<-p_lambda(theta = beta_j0, lambda = lambda) +
#       0.5 * (beta_j^2 - beta_j0^2) * p_lambda_d(theta = beta_j0, a = a, lambda = lambda)/abs(beta_j0)
#   }
# }



#define the log-likelihood function
loglikelihood_SCAD<-function(X, y, b){
  linear_comb<-as.vector(X %*% b)
  ll<-sum(y*linear_comb) + sum(log(1/(1+exp(linear_comb)))) - nrow(X)*sum(p_lambda(theta = b))
  return (ll)
}


##初始化系數(shù)
#b0<-rep(0,length(b_real))

#b0<- b_real+rnorm(length(b_real), mean = 0, sd = 0.1)

##將無(wú)懲罰時(shí)的優(yōu)化結(jié)果作為初始值
b.best_GS<-b.best

b0<-b.best_GS

##b1用于記錄更新系數(shù)
b1<-b0

##b.best用于存放歷史最大似然值對(duì)應(yīng)系數(shù)
b.best_SCAD<-b0

# the initial value of loglikelihood

ll.old<-loglikelihood_SCAD(X = X,y = y, b = b0)


# initialize the difference between the two steps of theta
diff<-1  
#record the number of iterations
iter<-0
#set the threshold to stop iterations
epsi<-1e-10
#the maximum iterations  
max_iter<-100000
#初始化一個(gè)列表用于存放每一次迭代的系數(shù)結(jié)果
b_history<-list(data.frame(b0))

#初始化列表用于存放似然值
ll_list<-list(ll.old)



#######-------SCAD迭代---------
while(diff > epsi & iter < max_iter){
  for(j in 1:length(b_real)){
    if(abs(b0[j]) < 1e-06){
      next()
    }else{
      
      #線性部分
      linear_comb<-as.vector(X %*% b0)
      
      #分子
      nominator<-sum(y*X[,j] - X[,j] * exp(linear_comb)/(1+exp(linear_comb))) +
        nrow(X)*b0[j]*p_lambda_d(theta = b0[j])/abs(b0[j])
      
      
      #分母,即二階導(dǎo)部分
      denominator<- -sum(X[,j]^2 * exp(linear_comb)/(1+exp(linear_comb))^2) +
        nrow(X)*p_lambda_d(theta = b0[j])/abs(b0[j])
      
      #2-(3) :更新b0[j]
      b0[j]<-b0[j] - nominator/denominator
      
      #2-(4)
      if(abs(b0[j]) < 1e-06){
        b0[j] <- 0
      }
      
      #       #更新似然值
      #       ll.new<- loglikelihood_SCAD(X = X, y = y, b = b0)
      #       
      #       
      #       
      #       #若似然值有所增加，則將當(dāng)前系數(shù)保存
      #       if(ll.new > ll.old){
      #         #更新系數(shù)
      #         b.best_SCAD[j]<-b0[j]
      #       }
      #       
      #       #求差異
      #       diff<- abs((ll.new - ll.old)/ll.old)
      #       ll.old <- ll.new
      #       iter<- iter+1
      #       b_history[[iter]]<-data.frame(b0)
      #       ll_list[[iter]]<-ll.old
      #       ##當(dāng)達(dá)到停止條件時(shí)，跳出循環(huán)
      #       if(diff < epsi){
      #         break
      #       }
      #       
      
    }
  }
 
  #更新似然值
  ll.new<- loglikelihood_SCAD(X = X, y = y, b = b0)
 
 
 
  #若似然值有所增加，則將當(dāng)前系數(shù)保存
  if(ll.new > ll.old){
    #更新系數(shù)
    b.best_SCAD<-b0
  }
 
  #求差異
  diff<- abs((ll.new - ll.old)/ll.old)
  ll.old <- ll.new
  iter<- iter+1
  b_history[[iter]]<-data.frame(b0)
  ll_list[[iter]]<-ll.old
 
 
}


b_hist<-do.call(rbind,b_history)
#b_hist

ll_hist<-do.call(rbind,ll_list)
#ll_hist

#
iter

##
ll.best<-max(ll_hist)
ll.best

##
b.best_SCAD

##對(duì)比
cbind(coeff_glm,b.best,b.best_SCAD,b_real)


##----------ncvreg驗(yàn)證-----------
library(ncvreg)
my_ncvreg<-ncvreg(X,y,family = c("binomial"),penalty = c("SCAD"),lambda = 2)
my_ncvreg$beta


my_ncvreg<-ncvreg(X,y,family = c("binomial"),penalty = c("SCAD"))
summary(my_ncvreg)
my_ncvreg$beta

###用cv找最優(yōu)的lambda

scad_cv<-cv.ncvreg(X,y,family = c("binomial"),penalty='SCAD')
scad_cv$lambda.min
mySCAD=ncvreg(X,y,family = c("binomial"),penalty='SCAD',lambda=scad_cv$lambda.min)

summary(mySCAD)

ncv_SCAD<-mySCAD$beta[-1]

##對(duì)比
myFinalResults<-cbind(無(wú)懲罰項(xiàng)回歸=coeff_glm, GS迭代 = b.best,
                            GS_SCAD迭代 = b.best_SCAD, ncvreg = ncv_SCAD,真實(shí)值 = b_real)
save(myFinalResults,file = "myFinalResults.rda")

想深入學(xué)習(xí)統(tǒng)計(jì)學(xué)知識(shí)，為數(shù)據(jù)分析筑牢根基？那快來(lái)看看統(tǒng)計(jì)學(xué)極簡(jiǎn)入門(mén)課程！

學(xué)習(xí)入口：https://edu.cda.cn/goods/show/3386?targetId=5647&preview=0

課程由專(zhuān)業(yè)數(shù)據(jù)分析師打造，完全免費(fèi)，60 天有效期且隨到隨學(xué)。它用獨(dú)特思路講重點(diǎn)，從數(shù)據(jù)種類(lèi)到統(tǒng)計(jì)學(xué)體系，內(nèi)容通俗易懂。學(xué)完它，能讓你輕松入門(mén)統(tǒng)計(jì)學(xué)，還能提升數(shù)據(jù)分析能力。趕緊點(diǎn)擊鏈接開(kāi)啟學(xué)習(xí)，讓自己在數(shù)據(jù)領(lǐng)域更上一層樓！