NonConvexPenalizedRegression.jl

Regularization Paths for SCAD and MCP Penalized Regression Models

This is a quick, naive and partial Julia translation of the R package ncvreg. Only gaussian family is translated but if you need more, just ask, i can do it if it is useful for someone.

Algorithm is described in Breheny P and Huang J (2011) "Coordinate descent algorithms for nonconvex penalized regression, with applications to biological feature selection". Annals of Applied Statistics, 5: 232–253

I needed to do regression with SCAD penalty and I can't find it in any Julia package. Perharps it is now implemented in MLJLinearModels.jl.

julia> using LinearAlgebra 
julia> using Random
julia> using NonConvexPenalizedRegression 
julia> using RCall
julia> rng = MersenneTwister(1234);
julia> n, p = 50, 5
(50, 5)

julia> X = randn(rng, n, p)              # feature matrix
50×5 Array{Float64,2}:
  0.867347   -1.22672     0.183976    0.0377058   1.48494
 -0.901744   -0.541716   -1.27635    -0.490009    1.23969
  ⋮
 -1.00978    -1.66323    -0.744522    0.427383   -1.37986
 -0.543805   -0.521229   -0.191176   -0.492253   -0.984217

julia> a0 = collect(1:p)                # ground truths
5-element Array{Int64,1}:
 1
 2
 3
 4
 5

julia> y = X * a0 + 0.1 * randn(n) # generate response
50-element Array{Float64,1}:
   6.411769869798991
  -1.59817739694925
   ⋮
 -11.769008550241434
  -9.107294931708777

julia> XX = hcat(X, randn(rng, n, p))
50×10 Array{Float64,2}:
  0.867347   -1.22672     0.183976    0.0377058  …   1.69129     0.969694    0.222167    -0.953909
  ⋮                                              ⋱
 -0.543805   -0.521229   -0.191176   -0.492253      -2.5788     -0.329958    0.00775707  -0.370354

julia> @rput XX
50×10 Array{Float64,2}:
  0.867347   -1.22672     0.183976    0.0377058  …   1.69129     0.969694    0.222167    -0.953909
  ⋮                                              ⋱
 -0.543805   -0.521229   -0.191176   -0.492253      -2.5788     -0.329958    0.00775707  -0.370354

julia> @rput y
50-element Array{Float64,1}:
   6.411769869798991
  -1.59817739694925
   ⋮
 -11.769008550241434
  -9.107294931708777

julia> R"library(ncvreg)"
RObject{StrSxp}
[1] "ncvreg"    "stats"     "graphics"  "grDevices" "utils"     "datasets"
[7] "methods"   "base"

julia> R"scad <- coef(ncvreg(XX, y, lambda=0.2, penalty='SCAD', eps=.0001))"
RObject{RealSxp}
 (Intercept)           V1           V2           V3           V4           V5
-0.003322903  1.025666431  2.000108987  2.983498545  3.997543804  4.982264144
          V6           V7           V8           V9          V10
 0.000000000  0.000000000  0.000000000  0.000000000  0.000000000

julia> @rget scad
11-element Array{Float64,1}:
 -0.0033229032078964105
  1.0256664305062173
  2.000108987156345
  2.9834985454557157
  3.997543803872521
  4.982264143916537
  0.0
  0.0
  0.0
  0.0
  0.0

julia> println( " R scad = $scad")
 R scad = [-0.0033229032078964105, 1.0256664305062173, 2.000108987156345, 2.9834985454557157, 3.997543803872521, 4.982264143916537, 0.0, 0.0, 0.0, 0.0, 0.0]

julia> λ = [0.2]
1-element Array{Float64,1}:
 0.2

julia> scad = NonConvexPenalizedRegression.coef(SCAD(XX, y, λ))
SCAD([-0.003322960709765954; 1.0256660512338405; … ; 0.0; 0.0])

julia> println( " Julia scad = $scad")
 Julia scad = SCAD([-0.003322960709765954; 1.0256660512338405; 2.00010933635426; 2.983498839847109; 3.99754375703709; 4.982264100245242; 0.0; 0.0; 0.0; 0.0; 0.0])