R语言中实现线性回归+l_0范数变量选择的包abess

﹏ヽ暗。殇╰゛Y 2022-09-02 01:29 129阅读 0赞

### 文章目录 ###

*  模型理论
 *  函数介绍
 *  实例1：指定稀疏参数 s s s
 *  实例2：不指定稀疏参数 s s s
 *  实例3：根据边际最大效应进行预筛选

# 模型理论 #

该包实现了一种多项式复杂度的算法来利用下面的模型进行最优子集选择：  
 min ⁡ β ∈ R p 1 2 n ∥ y − X β ∥ 2 2 , s . t . ∥ β ∥ 0 ≤ s , \\min\_\{\\beta\\in\\mathcal\{R\}^p\}\\frac\{1\}\{2n\}\\|y-X\\beta\\|^2\_2, \\quad s.t.\\quad \\|\\beta\\|\_0\\leq s, β∈Rpmin2n1∥y−Xβ∥22,s.t.∥β∥0≤s,  
其中 ∥ ⋅ ∥ 2 \\|\\cdot\\|\_2 ∥⋅∥2是 l 2 l\_2 l2范数， ∥ β ∥ 0 = ∑ i = 1 p I ( β i ≠ 0 ) \\|\\beta\\|\_0=\\sum\_\{i=1\}^p I(\\beta\_i\\neq 0) ∥β∥0=∑i=1pI(βi=0)是 β \\beta β的 l 0 l\_0 l0范数，稀疏水平 s s s通常是一个未知的非负参数。

# 函数介绍 #

> abess(formula,data=,screening.num)：实现算法的函数  
> support.size :系数性参数  
> extract():从abess对象中提取最优模型的信息  
> coef(abess\_fit,support.size = best\_size, sparse = FALSE):获得abess\_fit模型中的系数估计，sparse=T时0系数用.表示，sparse=F时0系数仍然用0表示

# 实例1：指定稀疏参数 s s s #

产生一个样本量为300，协变量个数为1000的数据集，并服从模型  
 y = 3 X 1 + 1.5 X 2 + 2 X 5 + ϵ y=3X\_1+1.5X\_2+2X\_5+\\epsilon y=3X1\+1.5X2\+2X5\+ϵ， 然后利用`abess`进行求解。

library(abess)
    synthetic_data <- generate.data(n = 300, p = 1000, 
                                    beta = c(3, 1.5, 0, 0, 2, rep(0, 995)))
    dim(synthetic_data[["x"]])
    ## [1]  300 1000
    head(synthetic_data[["y"]])
    # [,1]
    # [1,] -4.063922
    # [2,]  3.855246
    # [3,] -3.041391
    # [4,] -1.081257
    # [5,]  4.986772
    # [6,]  4.470901
    
    dat <- cbind.data.frame("y" = synthetic_data[["y"]], 
                            synthetic_data[["x"]])
    abess_fit <- abess(y ~ ., data = dat, support.size = 3)
    abess_fit 
    # 
    # Call:
    #   abess.formula(formula = y ~ ., data = dat, support.size = 3)
    # 
    # support.size      dev      GIC
    # 1            3 1.308784 116.8114
    head(coef(abess_fit, sparse = T))
    # 6 x 1 sparse Matrix of class "dgCMatrix"
    # 3
    # (intercept) -0.01802179
    # x1           2.96418205
    # x2           1.45090693
    # x3           .         
    # x4           .         
    # x5           1.90592036
    head(coef(abess_fit, sparse = F))
    #               3
    # (intercept) -0.01802179
    # x1           2.96418205
    # x2           1.45090693
    # x3           0.00000000
    # x4           0.00000000
    # x5           1.90592036

注意，这里非稀疏的 β \\beta β的估计值和利用`lm()`的估计值相同.

# 实例2：不指定稀疏参数 s s s #

如果不指定 s s s该函数将根据GIC准则找到最优的模型。

abess_fit <- abess(y ~ ., data = dat)
    abess_fit
    # Call:
    #   abess.formula(formula = y ~ ., data = dat)
    # 
    # support.size        dev      GIC
    # 1             0 14.6129090 804.5716
    # 2             1  6.9808311 594.9777
    # 3             2  3.5881757 407.3478
    # 4             3  1.3087837 116.8114
    # 5             4  1.2706026 119.9566
    # 6             5  1.2314391 122.5915
    # 7             6  1.1937385 125.2908
    # 8             7  1.1612046 129.0284
    # 9             8  1.1315011 133.2819
    # 10            9  1.0970605 136.0360
    # 11           10  1.0638825 138.8504
    # 12           11  1.0349404 142.6034
    # 13           12  1.0089948 147.0139
    # 14           13  0.9828985 151.1800
    # 15           14  0.9586390 155.7099
    # 16           15  0.9361511 160.6159
    # 17           16  0.9142904 165.5546
    # 18           17  0.8936287 170.7246
    # 19           18  0.8713593 175.1811
    # 20           19  0.8508145 180.0503
    # 21           20  0.8300798 184.6759
    # 22           21  0.8109140 189.6952
    # 23           22  0.7894801 193.6863
    # 24           23  0.7719085 198.9611
    # 25           24  0.7542340 204.0394
    # 26           25  0.7362937 208.8446
    best_size <- abess_fit[["best.size"]]
    print(best_size)
    ## [1] 3
    head(coef(abess_fit, support.size = best_size, sparse = FALSE))
    #                3
    # (intercept) -0.01802179
    # x1           2.96418205
    # x2           1.45090693
    # x3           0.00000000
    # x4           0.00000000
    # x5           1.90592036

# 实例3：根据边际最大效应进行预筛选 #

在`abess`中使用参数`screening.num`可以实现预筛选

library(abess)
    working_directory <- getwd()
    if (file.exists("crime.rda")) {
        
      load("crime.rda")
    } else {
        
      crime_data_url <- "https://github.com/abess-team/abess/raw/master/R-package/data-raw/crime.rda"
      download.file(crime_data_url, "crime.rda")
      load(file.path(working_directory, "crime.rda"))
    }
    
    abess_fit <- abess(y ~ ., data = crime, screening.num = 1000)
    
    head(abess_fit[["screening.vars"]])
    # [1] "pctBlack"     "pctWhite"     "medIncome"    "pctWdiv"      "pctPubAsst"  
    # [6] "medFamIncome"
    
    best_model <- extract(abess_fit)
    str(abess_fit)
    # List of 13
    # $ beta          :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
    # .. ..@ i       : int [1:528] 442 442 3950 442 3267 3950 442 1747 3267 3950 ...
    # .. ..@ p       : int [1:34] 0 0 1 3 6 10 15 21 28 36 ...
    # .. ..@ Dim     : int [1:2] 5253 33
    # .. ..@ Dimnames:List of 2
    # .. .. ..$ : chr [1:5253] "pop" "perHoush" "pctBlack" "pctWhite" ...
    # .. .. ..$ : chr [1:33] "0" "1" "2" "3" ...
    # .. ..@ x       : num [1:528] -0.251 -0.204 4.195 -0.161 0.527 ...
    # .. ..@ factors : list()
    # $ intercept     : num [1:33] 599 2115 1774 1270 1045 ...
    # $ dev           : num [1:33] 381863 172847 155975 147048 141513 ...
    # $ tune.value    : num [1:33] 6426 6043 6004 5987 5981 ...
    # $ nobs          : int 500
    # $ nvars         : int 5253
    # $ family        : chr "gaussian"
    # $ tune.path     : chr "sequence"
    # $ support.size  : int [1:33] 0 1 2 3 4 5 6 7 8 9 ...
    # $ tune.type     : chr "GIC"
    # $ best.size     : int 8
    # $ screening.vars: chr [1:1000] "pctBlack" "pctWhite" "medIncome" "pctWdiv" ...
    # $ call          : language abess.formula(formula = y ~ ., data = crime, screening.num = 1000)
    # - attr(*, "class")= chr "abess"
    
    attributes(abess_fit)
    # $names
    # [1] "beta"           "intercept"      "dev"            "tune.value"    
    # [5] "nobs"           "nvars"          "family"         "tune.path"     
    # [9] "support.size"   "tune.type"      "best.size"      "screening.vars"
    # [13] "call"          
    # 
    # $class
    # [1] "abess"