Skip to contents

Generalizes lassosumRss() to support LASSO, MCP, SCAD, L0, L0L1, and L0L2 penalties. Uses coordinate descent on the objective $$\beta^T R \beta - 2 \beta^T z + \mathrm{penalty}(\beta)$$ where \(R\) is a (possibly pre-shrunk) LD matrix and \(z = \hat\beta / \sqrt{n}\).

Usage

penalizedRss(
  bhat,
  R,
  n,
  penalty = c("lasso", "MCP", "SCAD", "L0", "L0L1", "L0L2"),
  lambda = exp(seq(log(1e-04), log(0.1), length.out = 20)),
  gamma = NULL,
  alpha = 1,
  lambda0 = 0,
  lambda2 = 0,
  thr = 1e-04,
  maxiter = 10000,
  maxSwaps = 100
)

Arguments

bhat

Numeric vector of marginal effect estimates (length p).

R

The LD correlation matrix (a single matrix over the analysed window), as in lassosumRss().

n

GWAS sample size (positive scalar).

penalty

Penalty type: "lasso", "MCP", "SCAD", "L0", "L0L1", or "L0L2".

lambda

Numeric vector of regularization parameter values along which to trace a solution path (warm-started, largest-first). For LASSO/MCP/SCAD this is the primary penalty strength; for L0 variants it controls the L1 component.

gamma

Concavity parameter for MCP (default 3) or SCAD (default 3.7). Ignored for LASSO and L0 variants.

alpha

Elastic-net mixing for MCP/SCAD: \(l_1 = \lambda \alpha\), \(l_2 = \lambda (1-\alpha)\). Default 1 (pure L1, no ridge).

lambda0

L0 penalty weight (number of non-zeros). Required for L0 variants; ignored otherwise. Default 0.

lambda2

L2 penalty weight for L0L2 variant. Default 0.

thr

Convergence threshold. Default 1e-4.

maxiter

Maximum coordinate descent iterations per lambda. Default 10000.

maxSwaps

Maximum swap rounds for L0 variants. Default 100. Set to 0 to disable swaps.

Value

A list with components:

beta

p x length(lambda) matrix of coefficient estimates.

lambda

The lambda values used.

conv

Convergence indicators (1 = converged).

loss

Quadratic loss at each lambda.

fbeta

Full penalized objective at each lambda.

nparams

Number of non-zero coefficients at each lambda.

betaEst

Coefficient vector at the lambda minimizing fbeta.

Examples

set.seed(42)
p <- 10; n <- 100
bhat <- rnorm(p, sd = 0.1)
R <- diag(p)
# MCP
penalizedRss(bhat, R, n, penalty = "MCP")
#> $beta
#>                [,1]          [,2]          [,3]          [,4]          [,5]
#>  [1,]  0.0137095845  0.0137095845  0.0137095845  0.0137095845  0.0137095845
#>  [2,] -0.0056469817 -0.0056469817 -0.0056469817 -0.0056469817 -0.0056469817
#>  [3,]  0.0036312841  0.0036312841  0.0036312841  0.0036312841  0.0036312841
#>  [4,]  0.0063286260  0.0063286260  0.0063286260  0.0063286260  0.0063286260
#>  [5,]  0.0040426832  0.0040426832  0.0040426832  0.0040426832  0.0040426832
#>  [6,] -0.0010612452 -0.0010612452 -0.0010612452 -0.0010612452 -0.0009496679
#>  [7,]  0.0151152200  0.0151152200  0.0151152200  0.0151152200  0.0151152200
#>  [8,] -0.0009465904 -0.0009465904 -0.0009465904 -0.0009465904 -0.0007776857
#>  [9,]  0.0201842371  0.0201842371  0.0201842371  0.0201842371  0.0201842371
#> [10,] -0.0006271410 -0.0006271410 -0.0006271410 -0.0004942588 -0.0002985116
#>                [,6]          [,7]         [,8]         [,9]        [,10]
#>  [1,]  1.370958e-02  1.370958e-02  0.013709584  0.013709584  0.013709584
#>  [2,] -5.646982e-03 -5.646982e-03 -0.005646982 -0.005646982 -0.004515496
#>  [3,]  3.631284e-03  3.631284e-03  0.003535514  0.002697455  0.001491950
#>  [4,]  6.328626e-03  6.328626e-03  0.006328626  0.006328626  0.005537963
#>  [5,]  4.042683e-03  4.042683e-03  0.004042683  0.003314554  0.002109048
#>  [6,] -6.680954e-04 -2.630676e-04  0.000000000  0.000000000  0.000000000
#>  [7,]  1.511522e-02  1.511522e-02  0.015115220  0.015115220  0.015115220
#>  [8,] -4.961133e-04 -9.108539e-05  0.000000000  0.000000000  0.000000000
#>  [9,]  2.018424e-02  2.018424e-02  0.020184237  0.020184237  0.020184237
#> [10,] -1.693917e-05  0.000000e+00  0.000000000  0.000000000  0.000000000
#>               [,11]         [,12]       [,13]       [,14]       [,15] [,16]
#>  [1,]  0.0137095845  0.0123809845 0.008792977 0.003631808 0.000000000     0
#>  [2,] -0.0027814373 -0.0002870804 0.000000000 0.000000000 0.000000000     0
#>  [3,]  0.0000000000  0.0000000000 0.000000000 0.000000000 0.000000000     0
#>  [4,]  0.0038039038  0.0013095469 0.000000000 0.000000000 0.000000000     0
#>  [5,]  0.0003749896  0.0000000000 0.000000000 0.000000000 0.000000000     0
#>  [6,]  0.0000000000  0.0000000000 0.000000000 0.000000000 0.000000000     0
#>  [7,]  0.0151152200  0.0144894378 0.010901430 0.005740262 0.000000000     0
#>  [8,]  0.0000000000  0.0000000000 0.000000000 0.000000000 0.000000000     0
#>  [9,]  0.0201842371  0.0201842371 0.018504956 0.013343787 0.005919705     0
#> [10,]  0.0000000000  0.0000000000 0.000000000 0.000000000 0.000000000     0
#>       [,17] [,18] [,19] [,20]
#>  [1,]     0     0     0     0
#>  [2,]     0     0     0     0
#>  [3,]     0     0     0     0
#>  [4,]     0     0     0     0
#>  [5,]     0     0     0     0
#>  [6,]     0     0     0     0
#>  [7,]     0     0     0     0
#>  [8,]     0     0     0     0
#>  [9,]     0     0     0     0
#> [10,]     0     0     0     0
#> 
#> $lambda
#>  [1] 0.0001000000 0.0001438450 0.0002069138 0.0002976351 0.0004281332
#>  [6] 0.0006158482 0.0008858668 0.0012742750 0.0018329807 0.0026366509
#> [11] 0.0037926902 0.0054555948 0.0078475997 0.0112883789 0.0162377674
#> [16] 0.0233572147 0.0335981829 0.0483293024 0.0695192796 0.1000000000
#> 
#> $conv
#>  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
#> 
#> $loss
#>  [1] -0.0009277110 -0.0009277110 -0.0009277110 -0.0009276934 -0.0009275620
#>  [6] -0.0009269812 -0.0009259487 -0.0009252863 -0.0009238932 -0.0009150743
#> [11] -0.0008840717 -0.0008396894 -0.0007790770 -0.0005875828 -0.0002039265
#> [16]  0.0000000000  0.0000000000  0.0000000000  0.0000000000  0.0000000000
#> 
#> $fbeta
#>  [1] -9.274110e-04 -9.270903e-04 -9.264266e-04 -9.250887e-04 -9.225100e-04
#>  [6] -9.177926e-04 -9.088671e-04 -8.912145e-04 -8.575437e-04 -7.997471e-04
#> [11] -7.092633e-04 -5.614662e-04 -3.590607e-04 -1.494649e-04 -2.336194e-05
#> [16]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
#> 
#> $nparams
#>  [1] 10 10 10 10 10 10  9  7  7  7  6  5  3  3  1  0  0  0  0  0
#> 
#> $betaEst
#>  [1]  0.0137095845 -0.0056469817  0.0036312841  0.0063286260  0.0040426832
#>  [6] -0.0010612452  0.0151152200 -0.0009465904  0.0201842371 -0.0006271410
#> 
# SCAD
penalizedRss(bhat, R, n, penalty = "SCAD")
#> $beta
#>                [,1]          [,2]          [,3]          [,4]          [,5]
#>  [1,]  0.0137095845  0.0137095845  0.0137095845  0.0137095845  0.0137095845
#>  [2,] -0.0056469817 -0.0056469817 -0.0056469817 -0.0056469817 -0.0056469817
#>  [3,]  0.0036312841  0.0036312841  0.0036312841  0.0036312841  0.0036312841
#>  [4,]  0.0063286260  0.0063286260  0.0063286260  0.0063286260  0.0063286260
#>  [5,]  0.0040426832  0.0040426832  0.0040426832  0.0040426832  0.0040426832
#>  [6,] -0.0010612452 -0.0010612452 -0.0010612452 -0.0010377129 -0.0007536876
#>  [7,]  0.0151152200  0.0151152200  0.0151152200  0.0151152200  0.0151152200
#>  [8,] -0.0009465904 -0.0009465904 -0.0009465904 -0.0008556141 -0.0005715889
#>  [9,]  0.0201842371  0.0201842371  0.0201842371  0.0201842371  0.0201842371
#> [10,] -0.0006271410 -0.0006271410 -0.0005457056 -0.0003482533 -0.0001990078
#>                [,6]          [,7]         [,8]         [,9]         [,10]
#>  [1,]  1.370958e-02  1.370958e-02  0.013709584  0.013709584  0.0137095845
#>  [2,] -5.646982e-03 -5.646982e-03 -0.005646982 -0.004979307 -0.0032301425
#>  [3,]  3.631284e-03  3.631284e-03  0.002993912  0.001798303  0.0009946332
#>  [4,]  6.328626e-03  6.328626e-03  0.006328626  0.006061919  0.0043127541
#>  [5,]  4.042683e-03  4.042683e-03  0.003647310  0.002431304  0.0014060323
#>  [6,] -4.453969e-04 -1.753784e-04  0.000000000  0.000000000  0.0000000000
#>  [7,]  1.511522e-02  1.511522e-02  0.015115220  0.015115220  0.0151152200
#>  [8,] -3.307422e-04 -6.072359e-05  0.000000000  0.000000000  0.0000000000
#>  [9,]  2.018424e-02  2.018424e-02  0.020184237  0.020184237  0.0201842371
#> [10,] -1.129278e-05  0.000000e+00  0.000000000  0.000000000  0.0000000000
#>              [,11]         [,12]       [,13]       [,14]      [,15] [,16] [,17]
#>  [1,]  0.013519367  0.0099001043 0.005861985 0.002421206 0.00000000     0     0
#>  [2,] -0.001854292 -0.0001913869 0.000000000 0.000000000 0.00000000     0     0
#>  [3,]  0.000000000  0.0000000000 0.000000000 0.000000000 0.00000000     0     0
#>  [4,]  0.002535936  0.0008730313 0.000000000 0.000000000 0.00000000     0     0
#>  [5,]  0.000249993  0.0000000000 0.000000000 0.000000000 0.00000000     0     0
#>  [6,]  0.000000000  0.0000000000 0.000000000 0.000000000 0.00000000     0     0
#>  [7,]  0.015115220  0.0121325843 0.007267620 0.003826841 0.00000000     0     0
#>  [8,]  0.000000000  0.0000000000 0.000000000 0.000000000 0.00000000     0     0
#>  [9,]  0.020184237  0.0201833762 0.014977248 0.008895858 0.00394647     0     0
#> [10,]  0.000000000  0.0000000000 0.000000000 0.000000000 0.00000000     0     0
#>       [,18] [,19] [,20]
#>  [1,]     0     0     0
#>  [2,]     0     0     0
#>  [3,]     0     0     0
#>  [4,]     0     0     0
#>  [5,]     0     0     0
#>  [6,]     0     0     0
#>  [7,]     0     0     0
#>  [8,]     0     0     0
#>  [9,]     0     0     0
#> [10,]     0     0     0
#> 
#> $lambda
#>  [1] 0.0001000000 0.0001438450 0.0002069138 0.0002976351 0.0004281332
#>  [6] 0.0006158482 0.0008858668 0.0012742750 0.0018329807 0.0026366509
#> [11] 0.0037926902 0.0054555948 0.0078475997 0.0112883789 0.0162377674
#> [16] 0.0233572147 0.0335981829 0.0483293024 0.0695192796 0.1000000000
#> 
#> $conv
#>  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
#> 
#> $loss
#>  [1] -0.0009277110 -0.0009277110 -0.0009277044 -0.0009276244 -0.0009272925
#>  [6] -0.0009265732 -0.0009257482 -0.0009247329 -0.0009188221 -0.0009014867
#> [11] -0.0008689195 -0.0008128306 -0.0006735436 -0.0004415435 -0.0001437383
#> [16]  0.0000000000  0.0000000000  0.0000000000  0.0000000000  0.0000000000
#> 
#> $fbeta
#>  [1] -9.272410e-04 -9.267385e-04 -9.257101e-04 -9.236947e-04 -9.200037e-04
#>  [6] -9.131254e-04 -8.995113e-04 -8.728295e-04 -8.258891e-04 -7.522001e-04
#> [11] -6.309980e-04 -4.447533e-04 -2.512276e-04 -9.964324e-05 -1.557462e-05
#> [16]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
#> 
#> $nparams
#>  [1] 10 10 10 10 10 10  9  7  7  7  6  5  3  3  1  0  0  0  0  0
#> 
#> $betaEst
#>  [1]  0.0137095845 -0.0056469817  0.0036312841  0.0063286260  0.0040426832
#>  [6] -0.0010612452  0.0151152200 -0.0009465904  0.0201842371 -0.0006271410
#> 
# L0
penalizedRss(bhat, R, n, penalty = "L0", lambda0 = 0.01,
              lambda = c(0))
#> $beta
#>       [,1]
#>  [1,]    0
#>  [2,]    0
#>  [3,]    0
#>  [4,]    0
#>  [5,]    0
#>  [6,]    0
#>  [7,]    0
#>  [8,]    0
#>  [9,]    0
#> [10,]    0
#> 
#> $lambda
#> [1] 0
#> 
#> $conv
#> [1] 1
#> 
#> $loss
#> [1] 0
#> 
#> $fbeta
#> [1] 0
#> 
#> $nparams
#> [1] 0
#> 
#> $betaEst
#>  [1] 0 0 0 0 0 0 0 0 0 0
#>