Getting started with lasso2

For demonstration purposes we use the prostate cancer data set, which has been widely applied to demonstrate the lasso and related techniques.

To load prostate cancer data:

. insheet using ///
    https://web.stanford.edu/~hastie/ElemStatLearn/datasets/prostate.data,  ///
    clear tab

General demonstration

By default, lasso2 uses the lasso estimator (i.e., alpha(1)). Like Stata’s regress, lasso2 expects the dependent variable to be named first followed by a list of predictors.

. lasso2 lpsa lcavol lweight age lbph svi lcp gleason pgg45

  Knot|  ID     Lambda    s      L1-Norm        EBIC     R-sq   | Entered/removed
------+---------------------------------------------------------+----------------
     1|   1  163.62492     1     0.00000     31.41226   0.0000  | Added _cons.
     2|   2  149.08894     2     0.06390     26.66962   0.0916  | Added lcavol.
     3|   9   77.73509     3     0.40800    -12.63533   0.4221  | Added svi.
     4|  11   64.53704     4     0.60174    -18.31145   0.4801  | Added lweight.
     5|  21   25.45474     5     1.35340    -42.20238   0.6123  | Added pgg45.
     6|  22   23.19341     6     1.39138    -38.93672   0.6175  | Added lbph.
     7|  29   12.09306     7     1.58269    -39.94418   0.6389  | Added age.
     8|  35    6.92010     8     1.71689    -38.84649   0.6516  | Added gleason.
     9|  41    3.95993     9     1.83346    -35.69248   0.6567  | Added lcp.
Use 'long' option for full output. Type e.g. 'lasso2, lic(ebic)' to run the model selected by EBIC.

lasso2 obtains the solution for a list of values (see third column). As decreases, predictors are added to the model. The last column on the right indicates which predictors enter or leave the active set.

Plotting

We can plot the coefficient path using plotpath(). The option accepts lnlambda, lambda and norm, which asks Stata to plot the coefficient estimates against , or the -norm, respectively.

For example, to plot the coefficients against :

. lasso2 lpsa lcavol lweight age lbph svi lcp gleason pgg45, ///
         plotpath(lnlambda) ///
         plotopt(legend(off)) ///
         plotlabel

plotopt() allows to specify additional plotting options that are passed on to Stata’s line command. In the example, legend(off) is used to suppress the legend. plotlabel triggers the display of variable labels next to the line.

The resulting graph looks as follows:

To plot the coefficients against the -norm, we can use:

. lasso2, plotpath(norm) ///
	  plotlabel ///
	  plotopt(legend(off))

By omitting the variable names (before the comma), we make use of the replay syntax. lasso2 uses the previous lasso2 results which avoids time-consuming re-estimation. This only works if lasso2 results are in memory.

This creates the following graph:

The lambda() option

We can use the lambda() option to estimate the model that corresponds to a specific value of .

. lasso2 lpsa lcavol lweight age lbph svi lcp gleason pgg45, ///
	 lambda(10)

---------------------------------------------------
	 Selected |           Lasso   Post-est OLS
------------------+--------------------------------
	   lcavol |       0.5000819      0.5234981
	  lweight |       0.5144276      0.6152349
	      age |      -0.0036627     -0.0190343
	     lbph |       0.0468469      0.0954908
	      svi |       0.5695171      0.6358643
	    pgg45 |       0.0017981      0.0035248
------------------+--------------------------------
   Partialled-out*|
------------------+--------------------------------
	    _cons |      -0.0014767      0.5214696
---------------------------------------------------

The output shows the lasso and post-estimation OLS estimates corresponding to . However, note that there is no justification for using over any other positive value. To pick the “optimal” value for , we can use cross-validation (see cvlasso), theory-driven penalization (rlasso) or information criteria as discussed below.

Estimators

The default estimator of lasso2 is the lasso, which corresponds to alpha(1). The alpha() option controls the elastic net parameter, which determines the relative contribution of (lasso-type) to (ridge-type) penalization.

For ridge regression, use alpha(0):

. lasso2 lpsa lcavol lweight age lbph svi lcp gleason pgg45, ///
         alpha(0) ///
	 lambda(500)
---------------------------------------------------
	 Selected |           Ridge   Post-est OLS
------------------+--------------------------------
	   lcavol |       0.1497346      0.5643413
	  lweight |       0.2497274      0.6220198
	      age |       0.0016636     -0.0212482
	     lbph |       0.0294061      0.0967125
	      svi |       0.2913161      0.7616733
	      lcp |       0.0687956     -0.1060509
	  gleason |       0.0771692      0.0492279
	    pgg45 |       0.0023278      0.0044575
------------------+--------------------------------
   Partialled-out*|
------------------+--------------------------------
	    _cons |       0.6322244      0.1815609
---------------------------------------------------

In contrast to the lasso, the ridge includes all predictors in the model. It does not perform variable selection, even if is large, whereas the elastic net can yield sparse solutions even when is small:

. lasso2 lpsa lcavol lweight age lbph svi lcp gleason pgg45, ///
         alpha(0.1) ///
	 lambda(500)
---------------------------------------------------
	 Selected |     Elastic net   Post-est OLS
	          |   (alpha=0.100)
------------------+--------------------------------
	   lcavol |       0.1222143      0.5378499
	  lweight |       0.1236962      0.6620155
	      svi |       0.1854247      0.6991923
	      lcp |       0.0424339     -0.0813594
	  gleason |       0.0116789      0.0322875
	    pgg45 |       0.0008686      0.0036868
------------------+--------------------------------
   Partialled-out*|
------------------+--------------------------------
	    _cons |       1.7319371     -1.1240093
---------------------------------------------------

Lastly, to employ square-root lasso, use the sqrt option:

. lasso2 lpsa lcavol lweight age lbph svi lcp gleason pgg45, ///
         sqrt ///
	 lambda(20)

Note that the alpha() and sqrt option are incompatible.

Information criteria

To select the “best” value of lambda, lasso2 offers four information criteria:

  • Akaike Information Criterion (),
  • Bayesian Information Criterion (),
  • Extended Bayesian Information Criterion ()
  • and Corrected Akaike Information Criterion ()

By default, the is shown in the output.

To estimate the model corresponding to the minimum information criterion, the lic() is used which accepts aic, bic, ebic and aicc. lic() can either be specified in the first lasso2 call or using the replay syntax (to avoid re-estimation).

. lasso2 lpsa lcavol lweight age lbph svi lcp gleason pgg45, ///
	 lic(aic)

The same can be achieved in two steps using the replay syntax:

. lasso2 lpsa lcavol lweight age lbph svi lcp gleason pgg45
. lasso2, lic(aic)

Use lambda=7.594796178345335 (selected by AIC).

---------------------------------------------------
	 Selected |           Lasso   Post-est OLS
------------------+--------------------------------
	   lcavol |       0.5057140      0.5234981
	  lweight |       0.5386738      0.6152349
	      age |      -0.0073599     -0.0190343
	     lbph |       0.0585468      0.0954908
	      svi |       0.5854749      0.6358643
	    pgg45 |       0.0022134      0.0035248
------------------+--------------------------------
   Partialled-out*|
------------------+--------------------------------
	    _cons |       0.1243026      0.5214696
---------------------------------------------------

As indicated in the output, the AIC selects lambda=7.59.

Predicted values

The predict post-estimation command allows to return residuals and predicted values. In the output below, xbhat1 is generated by re-estimating the model for . The noisily option triggers the display of the estimation results. xbhat2 is generated by linear approximation using the two beta estimates closest to lambda=10.

. lasso2 lpsa lcavol lweight age lbph svi lcp gleason pgg45
. cap drop xbhat1
. predict double xbhat1, xb l(10) noisily
. cap drop xbhat2
. predict double xbhat2, xb l(10) approx

Alternatively, we can explicitly run the model using lambda(10). If lasso2 is called with a scalar value, the subsequent predict command requires no lambda() option.

. lasso2 lpsa lcavol lweight age lbph svi lcp gleason pgg45, ///
	 lambda(10)
. cap drop xbhat3
. predict double xbhat3, xb

All three methods yield the same results. However note that the linear approximation is only exact for the lasso which is piecewise linear.

Adaptive lasso

The adaptive lasso relies on an initial estimator to calculate the penalty loadings. By default, lasso2 uses OLS as the initial estimator as originally suggested by Zou (2006). If the number of parameters exceeds the numbers of observations, univariate OLS is used; see Huang et al. (2008).

. lasso2 lpsa lcavol lweight age lbph svi lcp gleason pgg45, ///
         adaptive
. mat list e(Ups)

See the OLS estimates for comparison.

. reg lpsa lcavol lweight age lbph svi lcp gleason pgg45

Other initial estimators such as ridge regression are possible.

. lasso2 lpsa lcavol lweight age lbph svi lcp gleason pgg45, ///
	 l(10) alpha(0)
. mat bhat_ridge = e(b)
. lasso2 lpsa lcavol lweight age lbph svi lcp gleason pgg45, 
	 adaptive adaloadings(bhat_ridge)
. mat list e(Ups)

More

Please check the help file for more information and examples.

. help lasso2