* certification script for * lassopack package 1.4.3 5jan2024, aa/ms * parts of the script use R's glmnet and Matlab code "SqrtLassoIterative.m". set more off cscript "lasso2" adofile lasso2 lasso2_p lassoutils clear all capture log close set rmsg on program drop _all log using cs_lasso2,replace about which lasso2 which lasso2_p which lassoutils * data source //global prostate prostate.data global prostate https://web.stanford.edu/~hastie/ElemStatLearn/datasets/prostate.data * simple ridge regression program cap program drop estridge program define estridge, rclass syntax varlist , Lambda(real) [NOCONStant] local yvar : word 1 of `varlist' local xvars : list varlist - yvar qui putmata X=(`xvars') y=(`yvar'), replace mata: n=rows(y) if ("`noconstant'"=="") { mata: X=X:-mean(X) mata: y=y:-mean(y) } mata: p=cols(X) mata: beta=lusolve(X'X+(`lambda')/2*I(p),X'y) tempname bhat mata: st_matrix("`bhat'",beta') mat list `bhat' return matrix bhat = `bhat' end cap program drop comparemat program define comparemat , rclass syntax anything [, tol(real 10e-3)] local A : word 1 of `0' local B : word 2 of `0' tempname Amat Bmat mat `Amat' = `A' mat `Bmat' = `B' local diff=mreldif(`Amat',`Bmat') di as text "mreldif=`diff'. tolerance = `tol'" mat list `Amat' mat list `Bmat' return scalar mreldif = `diff' assert `diff'<`tol' end * program to compare two vectors using col names cap program drop comparevec program define comparevec , rclass syntax anything [, tol(real 10e-3)] local A : word 1 of `0' local B : word 2 of `0' tempname Amat Bmat mat `Amat' = `A' mat `Bmat' = `B' local Anames: colnames `Amat' local Bnames: colnames `Bmat' local maxdiff = 0 local num = 0 foreach var of local Anames { local aix = colnumb(`Amat',"`var'") local bix = colnumb(`Bmat',"`var'") //di `aix' //di `bix' local thisdiff=reldif(el(`Amat',1,`aix'),el(`Bmat',1,`bix')) if `thisdiff'>`maxdiff' { local diff = `thisdiff' } local num=`num'+1 } di as text "Max rel dif = `maxdiff'. tolerance = `tol'" mat list `Amat' mat list `Bmat' return scalar maxdiff = `maxdiff' assert `maxdiff'<`tol' end ******************************************************************************** *** replicate glmnet *** ******************************************************************************** sysuse auto, clear drop if rep78==. global model price mpg-foreign // # the following R code was run using ‘glmnet’ version 4.0-2 /* library("glmnet") library("haven") library("tidyverse") auto <- read_dta("http://www.stata-press.com/data/r9/auto.dta") auto <- auto %>% drop_na() n <- nrow(auto) price <- auto$price X <- auto[,c("mpg","rep78","headroom","trunk", "weight","length","turn", "displacement","gear_ratio","foreign")] X <- X %>% mutate(foreign = as.integer(foreign)) %>% as.matrix() */ // single lambda, lasso /* > r<-glmnet(X,price,alpha=1,lambda=1000,thresh=1e-15) > t(coef(r)) 1 x 11 sparse Matrix of class "dgCMatrix" [[ suppressing 11 column names ‘(Intercept)’, ‘mpg’, ‘rep78’ ... ]] s0 4336.84 . . . . 0.3819041 . . 3.289185 . . */ mat G = 0.3819041, 3.289185, 4336.84 lasso2 $model, lambda(1000) lglmnet assert mreldif(e(b),G) <1e-5 // single lambda, ridge /* > r<-glmnet(X,price,alpha=0,lambda=1000,thresh=1e-15) > t(coef(r)) 1 x 11 sparse Matrix of class "dgCMatrix" [[ suppressing 11 column names ‘(Intercept)’, ‘mpg’, ‘rep78’ ... ]] s0 3425.776 -57.0062 296.9702 -392.4219 23.70664 s0 0.9938447 5.562674 -24.18045 8.855065 -536.018 s0 1726.658 */ mat G = -57.0062, 296.9702, -392.4219, 23.70664, 0.9938447, 5.562674, /// -24.18045, 8.855065, -536.018, 1726.65, 3425.776 lasso2 $model, alpha(0) lambda(1000) lglmnet assert mreldif(e(b),G) <1e-5 // single lambda, elastic net /* > r<-glmnet(X,price,alpha=0.5,lambda=1000,thresh=1e-15) > t(coef(r)) 1 x 11 sparse Matrix of class "dgCMatrix" [[ suppressing 11 column names ‘(Intercept)’, ‘mpg’, ‘rep78’ ... ]] s0 2883.129 -5.096381 . . . 0.692915 . . 5.942777 . s0 308.225 */ mat G = -5.096381, 0.692915, 5.942777, 308.225, 2883.129 lasso2 $model, alpha(0.5) lambda(1000) lglmnet assert mreldif(e(b),G) <1e-5 // lasso, lambda grid /* > r<-glmnet(X,price,alpha=1,nlambda=5,thres=1e-15) > r$lambda [1] 1584.1894208 158.4189421 15.8418942 1.5841894 [5] 0.1584189 > r$dev.ratio [1] 0.0000000 0.5272556 0.5970093 0.5988338 0.5988521 */ mat L = 1584.1894208, 158.4189421, 15.8418942, 1.5841894, 0.1584189 mat L = L' mat D = 0.0000000, 0.5272556, 0.5970093, 0.5988338, 0.5988521 mat D = D' lasso2 $model, lglmnet lcount(5) assert mreldif(e(lambdamat0),L) < 1e-5 assert mreldif(e(rsq),D) < 1e-5 // ridge, lambda grid /* > r<-glmnet(X,price,alpha=0,nlambda=5,thres=1e-15) > r$lambda [1] 1584189.4208 158418.9421 15841.8942 1584.1894 [5] 158.4189 > r$dev.ratio [1] 2.777378e-36 4.347508e-02 2.060125e-01 4.453574e-01 [5] 5.744210e-01 */ mat L = 1584189.4208, 158418.9421, 15841.8942, 1584.1894, 158.4189 mat L = L' // first R-sq of glmnet appears to be wrong - see below * mat D = 2.777378e-36, 4.347508e-02, 2.060125e-01, 4.453574e-01, 5.744210e-01 * mat D = D' lasso2 $model, alpha(0) lglmnet lcount(5) long assert mreldif(e(lambdamat0),L) < 1e-5 * assert mreldif(e(rsq),D) < 1e-5 // single lambda, ridge - see above /* > r<-glmnet(X,price,alpha=0,lambda=1584189.4208,thresh=1e-15) > r$dev.ratio [1] 0.004941719 */ lasso2 $model, alpha(0) lglmnet lambda(1584189.4208) assert reldif(e(r2),0.004941719) < 1e-5 // elastic net, grid of 5 /* > r<-glmnet(X,price,alpha=0.5,nlambda=5,thresh=1e-15) > r$lambda [1] 3168.3788416 316.8378842 31.6837884 3.1683788 [5] 0.3168379 > r$dev.ratio [1] 0.0000000 0.5088361 0.5942304 0.5987942 0.5988517 */ mat L = 3168.3788416, 316.8378842, 31.6837884, 3.1683788, 0.3168379 mat L = L' mat D = 0.0000000, 0.5088361, 0.5942304, 0.5987942, 0.5988517 mat D = D' lasso2 $model, alpha(0.5) lglmnet lcount(5) long assert mreldif(e(lambdamat0),L) < 1e-5 assert mreldif(e(rsq),D) < 1e-5 // single lambda, lasso, nocons /* > r<-glmnet(X,price,alpha=1,lambda=1000,intercept=FALSE,thresh=1e-15) > t(coef(r)) 1 x 11 sparse Matrix of class "dgCMatrix" [[ suppressing 11 column names ‘(Intercept)’, ‘mpg’, ‘rep78’ ... ]] s0 . . . . . 0.379965 26.1441 . . . . */ mat G = 0.379965, 26.1441 lasso2 $model, lambda(1000) lglmnet nocons assert mreldif(e(b),G) <1e-5 // single lambda, lasso, no standardisation /* > r<-glmnet(X,price,alpha=1,lambda=1000,standardize=FALSE,thresh=1e-15) > t(coef(r)) 1 x 11 sparse Matrix of class "dgCMatrix" [[ suppressing 11 column names ‘(Intercept)’, ‘mpg’, ‘rep78’ ... ]] s0 10119.16 . . . . 3.549222 -62.23642 -106.0307 s0 6.079465 . . */ mat G = 3.549222, -62.23642, -106.0307, 6.079465, 10119.16 lasso2 $model, lambda(1000) lglmnet nostd assert mreldif(e(b),G) <1e-5 // single lambda, lasso, no standardisation, noconstant /* > r<-glmnet(X,price,alpha=1,lambda=1000,standardize=FALSE,intercept=FALSE,thresh=1e-15) > t(coef(r)) 1 x 11 sparse Matrix of class "dgCMatrix" [[ suppressing 11 column names ‘(Intercept)’, ‘mpg’, ‘rep78’ ... ]] s0 . 6.830147 . . . 1.85882 -2.334519 . 3.886964 . . */ mat G = 6.830147, 1.85882, -2.334519, 3.886964 lasso2 $model, lambda(1000) lglmnet nostd nocons // note looser tolerance assert mreldif(e(b),G) <1e-4 // single lambda, lasso, mpg and foreign unpenalized /* > p.fac = rep(1, 10) > p.fac[c(1, 10)] = 0 > r<-glmnet(X,price,alpha=1,lambda=500,penalty.factor=p.fac, thresh=1e-15) > t(coef(r)) 1 x 11 sparse Matrix of class "dgCMatrix" [[ suppressing 11 column names ‘(Intercept)’, ‘mpg’, ‘rep78’ ... ]] s0 9140.732 -225.1226 . . . . . . 6.065554 . 1962.096 */ mat G = -225.1226, 6.065554, 1962.096, 9140.732 lasso2 $model, lambda(500) lglmnet notpen(mpg foreign) assert mreldif(e(b),G) <1e-5 // single lambda, ridge, mpg and foreign unpenalized /* > p.fac = rep(1, 10) > p.fac[c(1, 10)] = 0 > r<-glmnet(X,price,alpha=0,lambda=500,penalty.factor=p.fac, thresh=1e-15) > t(coef(r)) 1 x 11 sparse Matrix of class "dgCMatrix" [[ suppressing 11 column names ‘(Intercept)’, ‘mpg’, ‘rep78’ ... ]] s0 4420.934 -65.4987 167.7168 -446.8571 22.46099 1.26494 s0 2.118738 -25.1642 10.62293 -910.106 3185.862 */ mat G = -65.4987, 167.7168, -446.8571, 22.46099, 1.26494, /// 2.118738, -25.1642, 10.62293, -910.106, 3185.862, 4420.934 lasso2 $model, lambda(500) alpha(0) lglmnet notpen(mpg foreign) assert mreldif(e(b),G) <1e-5 // single lambda, elastic net, mpg and foreign unpenalized /* > p.fac = rep(1, 10) > p.fac[c(1, 10)] = 0 > r<-glmnet(X,price,alpha=0.5,lambda=500,penalty.factor=p.fac, thresh=1e-15) > t(coef(r)) 1 x 11 sparse Matrix of class "dgCMatrix" [[ suppressing 11 column names ‘(Intercept)’, ‘mpg’, ‘rep78’ ... ]] s0 4089.962 -132.7721 . . . 0.7644731 . . 8.896648 . 2639.545 */ mat G = -132.7721, 0.7644731, 8.896648, 2639.545, 4089.962 lasso2 $model, lambda(500) alpha(0.5) lglmnet notpen(mpg foreign) assert mreldif(e(b),G) <1e-5 // single lambda, lasso, misc penalty loadings /* > p.fac = rep(1, 10) > p.fac[c(1, 10)] = 4 > p.fac[c(2,9)] = 3 > p.fac[c(3,8)] = 2 > r<-glmnet(X,price,alpha=1,lambda=400,penalty.factor=p.fac, thresh=1e-15) > t(coef(r)) 1 x 11 sparse Matrix of class "dgCMatrix" [[ suppressing 11 column names ‘(Intercept)’, ‘mpg’, ‘rep78’ ... ]] s0 1254.559 . . . . 2.136852 . -42.89822 . . 393.3358 */ mat G = 2.136852, -42.89822, 393.3358, 1254.559 mat psi = 4, 3, 2, 1, 1, 1, 1, 2, 3, 4 lasso2 $model, lambda(400) lglmnet ploadings(psi) assert mreldif(e(b),G) <1e-5 // Elastic net, single lambda, misc penalty loadings /* > p.fac = rep(1, 10) > p.fac[c(1, 10)] = 4 > p.fac[c(2,9)] = 3 > p.fac[c(3,8)] = 2 > r<-glmnet(X,price,alpha=0.5,lambda=400,penalty.factor=p.fac,thresh=1e-15) > t(coef(r)) 1 x 11 sparse Matrix of class "dgCMatrix" [[ suppressing 11 column names ‘(Intercept)’, ‘mpg’, ‘rep78’ ... ]] s0 1232.917 . 91.49722 -256.4055 . 2.29446 . -78.49864 5.61806 . s0 1398.059 */ mat G = 91.49722, -256.4055, 2.29446, -78.49864, 5.61806, 1398.059, 1232.917 mat psi = 4, 3, 2, 1, 1, 1, 1, 2, 3, 4 lasso2 $model, lambda(400) alpha(0.5) lglmnet ploadings(psi) assert mreldif(e(b),G) <1e-5 // Elastic net, single lambda, misc penalty loadings, no standardization /* > p.fac = rep(1, 10) > p.fac[c(1, 10)] = 4 > p.fac[c(2,9)] = 3 > p.fac[c(3,8)] = 2 > r<-glmnet(X,price,alpha=0.5,lambda=400,penalty.factor=p.fac,standardize=FALSE,thresh=1e-15) > t(coef(r)) 1 x 11 sparse Matrix of class "dgCMatrix" [[ suppressing 11 column names ‘(Intercept)’, ‘mpg’, ‘rep78’ ... ]] s0 16039.13 -58.51055 322.0282 -218.2771 21.3079 4.093964 -72.31631 s0 -247.7182 8.463531 . . */ mat G = -58.51055, 322.0282, -218.2771, 21.3079, 4.093964, -72.31631, /// -247.7182, 8.463531, 16039.13 mat psi = 4, 3, 2, 1, 1, 1, 1, 2, 3, 4 lasso2 $model, lambda(400) alpha(0.5) lglmnet ploadings(psi) nostd // note looser tolerance assert mreldif(e(b),G) <1e-4 ******************************************************************************** *** lglmnet option - consistency with lasso2 default parameterization *** ******************************************************************************** * uses auto dataset * can change sd(y) and lambdas to suit * note lglmnet implies prestd // any lambda, any alpha sysuse auto, clear gen double y = price/10 local glmnetlambda_a = 100 local glmnetlambda_b = 10 drop if rep78==. qui sum y local sd = r(sd) * 1/sqrt( r(N)/(r(N)-1) ) local glmnetalpha = 0.5 local L1lambda_a = `glmnetlambda_a' * 2 * 69 local L2lambda_a = `L1lambda_a' / `sd' local L1lambda_b = `glmnetlambda_b' * 2 * 69 local L2lambda_b = `L1lambda_b' / `sd' di "glmnetlambda_a=`glmnetlambda_a' L1lambda_a=`L1lambda_a' L2lambda_a=`L2lambda_a'" di "glmnetlambda_b=`glmnetlambda_b' L1lambda_b=`L1lambda_b' L2lambda_b=`L2lambda_b'" // Lasso / alpha=1 lasso2 y mpg-foreign, lambda(`glmnetlambda_a') lglmnet local objfn = e(objfn) storedresults save glmnet e() lasso2 y mpg-foreign, lambda(`L1lambda_a') prestd storedresults compare glmnet e(), tol(1e-8) /// exclude(macros: lasso2opt scalar: niter lambda objfn slambda sobjfn) assert reldif(2*`objfn',e(objfn))<1e-8 // lambda list lasso2 y mpg-foreign, lambda(`glmnetlambda_a' `glmnetlambda_b') lglmnet storedresults save glmnet e() lasso2 y mpg-foreign, lambda(`L1lambda_a' `L1lambda_b') prestd storedresults compare glmnet e(), tol(1e-8) /// exclude(macros: lasso2opt /// scalar: niter lmax lmax0 lmin lmin0 /// laic laicc lbic lebic slaic slaicc slbic slebic /// matrix: lambdamat lambdamat0 slambdamat slambdamat0) // Enet local alpha=(`glmnetalpha'*`sd')/( (1-`glmnetalpha') + `glmnetalpha'*`sd' ) local lambda_a=(1-`glmnetalpha')*`L2lambda_a' + `glmnetalpha'*`L1lambda_a' local lambda_b=(1-`glmnetalpha')*`L2lambda_b' + `glmnetalpha'*`L1lambda_b' di "alpha=`alpha' lambda_a=`lambda_a' lambda_b=`lambda_b'" lasso2 y mpg-foreign, alpha(`glmnetalpha') lambda(`glmnetlambda_a') lglmnet local objfn = e(objfn) storedresults save glmnet e() lasso2 y mpg-foreign, alpha(`alpha') lambda(`lambda_a') prestd storedresults compare glmnet e(), tol(1e-8) /// exclude(macros: lasso2opt scalar: niter alpha lambda objfn slambda sobjfn) di 2*`objfn' di e(objfn) assert reldif(2*`objfn',e(objfn))<1e-8 // lambda list lasso2 y mpg-foreign, alpha(`glmnetalpha') lambda(`glmnetlambda_a' `glmnetlambda_b') lglmnet storedresults save glmnet e() lasso2 y mpg-foreign, alpha(`alpha') lambda(`lambda_a' `lambda_b') prestd storedresults compare glmnet e(), tol(1e-8) /// exclude(macros: lasso2opt /// scalar: niter lmax lmax0 lmin lmin0 /// laic laicc lbic lebic slaic slaicc slbic slebic alpha /// matrix: lambdamat lambdamat0 slambdamat slambdamat0) // Ridge / alpha=0 lasso2 y mpg-foreign, alpha(0) lambda(`glmnetlambda_a') lglmnet local objfn = e(objfn) storedresults save glmnet e() lasso2 y mpg-foreign, alpha(0) lambda(`L2lambda_a') prestd storedresults compare glmnet e(), tol(1e-8) /// exclude(macros: lasso2opt scalar: niter lambda objfn slambda sobjfn) di 2*`objfn' di e(objfn) assert reldif(2*`objfn',e(objfn))<1e-8 // lambda list lasso2 y mpg-foreign, alpha(0) lambda(`glmnetlambda_a' `glmnetlambda_b') lglmnet storedresults save glmnet e() lasso2 y mpg-foreign, alpha(0) lambda(`L2lambda_a' `L2lambda_b') prestd storedresults compare glmnet e(), tol(1e-8) /// exclude(macros: lasso2opt /// scalar: niter lmax lmax0 lmin lmin0 /// laic laicc lbic lebic slaic slaicc slbic slebic /// matrix: lambdamat lambdamat0 slambdamat slambdamat0) ******************************************************************************** *** Validate ploadings(.) option *** ******************************************************************************** * Default lasso2 behavior is to standardize, either on the fly or in advance. * Confirm ploadings(.) work correctly by providing SDs of X. * Behavior should be identical to standardization on-the-fly. sysuse auto, clear drop if rep78==. cap mat drop psi_sd foreach var of varlist mpg-foreign { qui sum `var' mat psi_sd = nullmat(psi_sd) , r(sd) * sqrt( (r(N)-1)/r(N) ) } // lasso lasso2 price mpg-foreign, lambda(1000) alpha(1) storedresults save nopload e() lasso2 price mpg-foreign, lambda(1000) alpha(1) pload(psi_sd) storedresults compare nopload e(), tol(1e-8) /// exclude(macros: lasso2opt) // ridge lasso2 price mpg-foreign, lambda(1000) alpha(0) storedresults save nopload e() lasso2 price mpg-foreign, lambda(1000) alpha(0) pload(psi_sd) storedresults compare nopload e(), tol(1e-8) /// exclude(macros: lasso2opt) // elastic net lasso2 price mpg-foreign, lambda(1000) alpha(0.5) storedresults save nopload e() lasso2 price mpg-foreign, lambda(1000) alpha(0.5) pload(psi_sd) storedresults compare nopload e(), tol(1e-8) /// exclude(macros: lasso2opt) ******************************************************************************** *** Validate adaptive lasso option *** ******************************************************************************** sysuse auto, clear drop if rep78==. cap mat drop sdvec // standardized variables used to get OLS with std coefs foreach var of varlist price mpg-foreign { cap drop `var'_sd qui sum `var', meanonly qui gen double `var'_sd = `var'-r(mean) qui sum `var' qui replace `var'_sd = `var'_sd * 1/r(sd) * sqrt( r(N)/(r(N)-1) ) mat sdvec = nullmat(sdvec), r(sd)/sqrt( r(N)/(r(N)-1) ) } mat ysd = sdvec[1,1] mat xsd = sdvec[1,2..11] // replicate adaptive lasso // use standardized coefficients + prestd qui reg price_sd mpg_sd-foreign_sd mata: st_matrix("psi",1:/abs(st_matrix("e(b)"))) mat psi = psi[1,1..10] mat psi2 = J(1,10,1) // lasso lasso2 price mpg-foreign, lambda(10000) adaptive prestd mat b=e(b) lasso2 price mpg-foreign, lambda(10000) ploadings(psi) prestd assert mreldif(b,e(b)) < 1e-7 // elastic net lasso2 price mpg-foreign, lambda(10000) adaptive prestd alpha(0.5) mat b=e(b) lasso2 price mpg-foreign, lambda(10000) ploadings(psi) ploadings2(psi2) prestd alpha(0.5) assert mreldif(b,e(b)) < 1e-7 // using adaloadings option + standardized coefficients + prestd qui reg price_sd mpg_sd-foreign_sd mata: st_matrix("psi",st_matrix("e(b)")) mat psi = psi[1,1..10] mat psi2 = J(1,10,1) // lasso lasso2 price mpg-foreign, lambda(10000) adaptive prestd mat b=e(b) lasso2 price mpg-foreign, lambda(10000) adaloadings(psi) prestd assert mreldif(b,e(b)) < 1e-7 // elastic net lasso2 price mpg-foreign, lambda(10000) adaptive prestd alpha(0.5) mat b=e(b) lasso2 price mpg-foreign, lambda(10000) adaloadings(psi) prestd alpha(0.5) assert mreldif(b,e(b)) < 1e-7 // use unstandardized coefficients and no prestandardization qui reg price mpg-foreign mata: st_matrix("psi",1:/abs(st_matrix("e(b)"))) // note that psi needs to be rescaled by sd(y) in order for lambda to be the same mat psi = psi[1,1..10] * ysd[1,1] mat psi2 = xsd // lasso lasso2 price mpg-foreign, lambda(10000) adaptive mat b=e(b) lasso2 price mpg-foreign, lambda(10000) ploadings(psi) assert mreldif(b,e(b)) < 1e-7 // elastic net lasso2 price mpg-foreign, lambda(10000) adaptive alpha(0.5) mat b=e(b) lasso2 price mpg-foreign, lambda(10000) ploadings(psi) ploadings2(psi2) alpha(0.5) assert mreldif(b,e(b)) < 1e-7 // use adaloadings option + unstandardized coefficients + no prestandardization qui reg price mpg-foreign mata: st_matrix("psi",st_matrix("e(b)")) // note that psi does NOT need to be rescaled by sd(y) - handled automatically mat psi = psi[1,1..10] mat psi2 = xsd // lasso lasso2 price mpg-foreign, lambda(10000) adaptive mat b=e(b) lasso2 price mpg-foreign, lambda(10000) adaloadings(psi) assert mreldif(b,e(b)) < 1e-7 // elastic net lasso2 price mpg-foreign, lambda(10000) adaptive alpha(0.5) mat b=e(b) lasso2 price mpg-foreign, lambda(10000) adaloadings(psi) alpha(0.5) assert mreldif(b,e(b)) < 1e-7 // confirm prestd and no prestd yield same estimates // lasso, theta=1 lasso2 price mpg-foreign, lambda(10000) adaptive mat b=e(b) lasso2 price mpg-foreign, lambda(10000) adaptive prestd assert mreldif(b,e(b)) < 1e-7 // lasso, theta=2 lasso2 price mpg-foreign, lambda(10000) adaptive adatheta(2) mat b=e(b) lasso2 price mpg-foreign, lambda(10000) adaptive prestd adatheta(2) assert mreldif(b,e(b)) < 1e-7 // elastic net lasso2 price mpg-foreign, lambda(10000) adaptive alpha(0.5) adatheta(2) mat b=e(b) lasso2 price mpg-foreign, lambda(10000) adaptive prestd alpha(0.5) adatheta(2) assert mreldif(b,e(b)) < 1e-7 * lglmnet version * lglmnet standardizes automatically unless overridden by nostd (=unitloadings) // replicate adaptive lasso with lglmnet and no standardization // nostd overrides standardization // use unstandardized coefficients qui reg price mpg-foreign mata: st_matrix("psi",1:/abs(st_matrix("e(b)"))) mat psi = psi[1,1..10] mat psi2 = J(1,10,1) // lasso lasso2 price mpg-foreign, lambda(1000) adaptive lglmnet nostd mat b=e(b) lasso2 price mpg-foreign, lambda(1000) ploadings(psi) nostd lglmnet assert mreldif(b,e(b)) < 1e-7 // elastic net lasso2 price mpg-foreign, lambda(1000) adaptive lglmnet nostd alpha(0.5) mat b=e(b) lasso2 price mpg-foreign, lambda(1000) ploadings(psi) ploadings2(psi2) nostd lglmnet alpha(0.5) assert mreldif(b,e(b)) < 1e-7 // use adaloadings option + unstandardized coefficients + no standardization qui reg price mpg-foreign mata: st_matrix("psi",st_matrix("e(b)")) mat psi = psi[1,1..10] // lasso lasso2 price mpg-foreign, lambda(10000) adaptive lglmnet nostd mat b=e(b) lasso2 price mpg-foreign, lambda(10000) adaloadings(psi) nostd lglmnet assert mreldif(b,e(b)) < 1e-7 // elastic net lasso2 price mpg-foreign, lambda(10000) adaptive lglmnet nostd alpha(0.5) mat b=e(b) lasso2 price mpg-foreign, lambda(10000) adaloadings(psi) nostd lglmnet alpha(0.5) assert mreldif(b,e(b)) < 1e-7 // replicate adaptive lasso with lglmnet and standardization // lglmnet standardizes by default unless overridden by nostd // use standardized coefficients // note dep var doesn't have to be standardized qui reg price mpg_sd-foreign_sd mata: st_matrix("psi",1:/abs(st_matrix("e(b)"))) mat psi = psi[1,1..10] mat psi2 = J(1,10,1) // lasso lasso2 price mpg-foreign, lambda(1000) adaptive lglmnet mat b=e(b) lasso2 price mpg-foreign, lambda(1000) ploadings(psi) lglmnet assert mreldif(b,e(b)) < 1e-7 // elastic net lasso2 price mpg-foreign, lambda(1000) adaptive lglmnet alpha(0.5) mat b=e(b) lasso2 price mpg-foreign, lambda(1000) ploadings(psi) ploadings2(psi2) lglmnet alpha(0.5) assert mreldif(b,e(b)) < 1e-7 // use standardized coefficients + adaloadings option // note dep var doesn't have to be standardized qui reg price mpg_sd-foreign_sd mata: st_matrix("psi",st_matrix("e(b)")) mat psi = psi[1,1..10] // lasso lasso2 price mpg-foreign, lambda(1000) adaptive lglmnet mat b=e(b) lasso2 price mpg-foreign, lambda(1000) adaloadings(psi) lglmnet assert mreldif(b,e(b)) < 1e-7 // elastic net lasso2 price mpg-foreign, lambda(1000) adaptive lglmnet alpha(0.5) mat b=e(b) lasso2 price mpg-foreign, lambda(1000) adaloadings(psi) lglmnet alpha(0.5) assert mreldif(b,e(b)) < 1e-7 ******************************************************************************** *** replicate sqrt-lasso Matlab program *** ******************************************************************************** * load example data insheet using "$prostate", tab clear global model lpsa lcavol lweight age lbph svi lcp gleason pgg45 // uses the Matlab code "SqrtLassoIterative.m" (available on request) lasso2 $model, sqrt l(40) unitload mat a=e(betaAll) /* ans = 0.3627 0 0 0 0 0 0 0.0103 1.7383 */ mat b = (0.3627,0,0,0,0,0,0,0.0103,1.7383) comparemat a b lasso2 $model, sqrt l(10) unitload mat a=e(betaAll) /* ans = 0.5771 0.1965 -0.0092 0.0773 0.0685 0 0 0.0063 1.3946 */ mat b = (0.5771,0.1965,-0.0092,0.0773,0.0685,0,0,0.0063,1.3946) comparemat a b lasso2 $model, sqrt l(1) unitload mat a=e(betaAll) /* ans = 0.5610 0.5774 -0.0196 0.0950 0.6700 -0.0766 0.0088 0.0049 0.5259 */ mat b = (0.5610, 0.5774,-0.0196,0.0950,0.6700,-0.0766,0.0088,0.0049,0.5259) comparemat a b ******************************************************************************** *** validation using Stata's elasticnet *** ******************************************************************************** // ridge (requires 2 grid points for some reason) cap noi elasticnet linear $model, alphas(0) grid(2, min(0.25)) lassoselect alpha = 0 lambda = 0.25 lassocoef, display(coef, penalized) mat b=e(b) // Stata lambda = 2N*lambda global L=2*e(N)*0.25 lasso2 $model, lambda($L) alpha(0) assert mreldif(b,e(b))<1e-7 // lasso cap noi elasticnet linear $model, alphas(1) grid(1, min(0.25)) lassoselect alpha = 1 lambda = 0.25 lassocoef, display(coef, penalized) mat b=e(b) // Stata lambda = 2N*lambda global L=2*e(N)*0.25 lasso2 $model, lambda($L) assert mreldif(b,e(b))<1e-7 // elastic net cap noi elasticnet linear $model, alphas(0.5) grid(1, min(0.25)) lassoselect alpha = 0.5 lambda = 0.25 lassocoef, display(coef, penalized) mat b=e(b) // Stata lambda = 2N*lambda global L=2*e(N)*0.25 lasso2 $model, lambda($L) alpha(0.5) assert mreldif(b,e(b))<1e-7 ******************************************************************************** *** norecover option *** ******************************************************************************** // partial() with constant lasso2 $model, partial(age) l(50 20 10) mat A = e(betas) mat A = A[2,1..9] lasso2 $model, l(20) partial(age) postall mat B = e(b) comparemat A B lasso2 $model, partial(age) l(50 20 10) nor mat A = e(betas) mat A = A[2,1..7] lasso2 $model, l(20) partial(age) nor postall mat B = e(b) comparemat A B // partial() with constant, unitloadings lasso2 $model, partial(age) l(50 20 10) unitl mat A = e(betas) mat A = A[2,1..9] lasso2 $model, l(20) partial(age) postall unitl mat B = e(b) comparemat A B lasso2 $model, partial(age) l(50 20 10) nor unitl mat A = e(betas) mat A = A[2,1..7] lasso2 $model, l(20) partial(age) nor postall unitl mat B = e(b) comparemat A B // no partial() w/ constant, unitloadings lasso2 $model, l(50 20 10) unitl mat A = e(betas) mat A = A[2,1..9] lasso2 $model, l(20) postall unitl mat B = e(b) comparemat A B lasso2 $model, l(50 20 10) nor unitl mat A = e(betas) mat A = A[2,1..9] lasso2 $model, l(20) nor postall unitl mat B = e(b) comparemat A B // no partial() w/o constant, unit loadings lasso2 $model, l(50 20 10) unitl nocons mat A = e(betas) mat A = A[2,1..8] lasso2 $model, l(20) postall unitl nocons mat B = e(b) comparemat A B lasso2 $model, l(50 20 10) nor unitl nocons mat A = e(betas) mat A = A[2,1..8] lasso2 $model, l(20) nor postall unitl nocons mat B = e(b) comparemat A B ******************************************************************************** *** options *** ******************************************************************************** cap lasso2 $model, alpha(0) sqrt if _rc != 198 { exit 1 } * // should say that lcount/lmax/lminr are being ignored lasso2 $model, lambda(10) lcount(10) lasso2 $model, lambda(10) lmax(100) lasso2 $model, lambda(10) lminr(0.01) // plotting only supported for lambda list lasso2 $model, lambda(10) plotpath(lambda) // incompatible options wrt penalty loadings cap lasso2 $model, ploadings(abc) adaptive if _rc != 198 { exit 1 } * cap lasso2 $model, ploadings(abc) adatheta(3) if _rc != 198 { exit 1 } * // var may not appear in partial() and notpen() cap lasso2 $model, partial(age svi lcp) notpen(age svi) if _rc != 198 { exit 1 } * // controls the output and content of e(b) lasso2 $model, l(20) displayall lasso2 $model, l(20) postall mat list e(b) lasso2 $model, l(20) displayall postall mat list e(b) ******************************************************************************** *** verify results are the same for scalar lambda vs lambda list *** ******************************************************************************** global lambdalist 150 130 100 80 60 30 10 5 3 1 * lasso lasso2 $model, l($lambdalist) mat A = e(betas) local j=1 foreach i of numlist $lambdalist { mat a = A[`j',1..9] lasso2 $model, l(`i') mat b = e(betaAll) comparemat a b local j=`j'+1 } * * lasso (w/o constant) lasso2 $model, l($lambdalist) nocons mat A = e(betas) local j=1 foreach i of numlist $lambdalist { mat a = A[`j',1..8] lasso2 $model, l(`i') nocons mat b = e(betaAll) comparemat a b local j=`j'+1 } * * post-lasso lasso2 $model, l($lambdalist) ols mat A = e(betas) local j=1 foreach i of numlist $lambdalist { mat a = A[`j',1..9] lasso2 $model, l(`i') mat b = e(betaAllOLS) comparemat a b local j=`j'+1 } * // replaced lambda=40 with lambda=60 // lambda=40 in path used beta=zeros from lambda=100 // as initial beta; would end up at poor solution global sqrtlambdalist 100 60 20 10 5 1 * sqrt-lasso lasso2 $model, l($sqrtlambdalist) sqrt mat A = e(betas) local j=1 foreach i of numlist $sqrtlambdalist { mat a = A[`j',1..8] di `i' lasso2 $model, l(`i') sqrt mat b = e(betaAll) mat b = b[1,1..8] comparemat a b local j=`j'+1 } * * post-sqrt-lasso ols lasso2 $model, l($sqrtlambdalist) sqrt ols mat A = e(betas) local j=1 foreach i of numlist $sqrtlambdalist { di "this lambda: `i'" mat a = A[`j',1..9] lasso2 $model, l(`i') sqrt ols mat b = e(betaAllOLS) comparemat a b local j=`j'+1 } * * ridge lasso2 $model, l($lambdalist) alpha(0) mat A = e(betas) local j=1 foreach i of numlist $lambdalist { mat a = A[`j',1..9] lasso2 $model, l(`i') alpha(0) mat b = e(betaAll) comparemat a b local j=`j'+1 } * * ols ridge lasso2 $model, l($lambdalist) alpha(0) ols mat A = e(betas) local j=1 foreach i of numlist $lambdalist { mat a = A[`j',1..9] lasso2 $model, l(`i') alpha(0) ols mat b = e(betaAllOLS) comparemat a b local j=`j'+1 } * * elastic net foreach ai of numlist 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 { lasso2 $model, l($lambdalist) alpha(`ai') mat A = e(betas) local j=1 foreach i of numlist $lambdalist { mat a = A[`j',1..9] lasso2 $model, l(`i') alpha(`ai') mat b = e(betaAll) comparemat a b local j=`j'+1 } } * * elastic net with ols foreach ai of numlist 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 { lasso2 $model, l($lambdalist) alpha(`ai') ols mat A = e(betas) local j=1 foreach i of numlist $lambdalist { mat a = A[`j',1..9] lasso2 $model, l(`i') alpha(`ai') ols mat b = e(betaAllOLS) comparemat a b local j=`j'+1 } } * ******************************************************************************** *** pre-estimation standardisation vs std on the fly *** ******************************************************************************** // lasso // standardisation using penalty loadings (default) lasso2 $model, l(10) mat A = e(beta) // pre-estimation standardisation of data lasso2 $model, l(10) prestd mat B = e(beta) comparemat A B , tol(10e-6) // lasso [nocons] // standardisation using penalty loadings (default) lasso2 $model, l(10) nocons mat A = e(beta) // pre-estimation standardisation of data lasso2 $model, l(10) prestd nocons mat B = e(beta) comparemat A B , tol(10e-6) // sqrt lasso // standardisation using penalty loadings (default) lasso2 $model, l(10) sqrt mat A = e(beta) // pre-estimation standardisation of data lasso2 $model, l(10) sqrt prestd mat B = e(beta) comparemat A B , tol(10e-6) // sqrt lasso [nocons] // standardisation using penalty loadings (default) lasso2 $model, l(10) sqrt nocons mat A = e(beta) // pre-estimation standardisation of data lasso2 $model, l(10) sqrt prestd nocons mat B = e(beta) comparemat A B , tol(10e-6) // elastic net foreach lam of numlist 1 10 50 150 160 { foreach ai of numlist 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 { // in original units // standardisation using penalty loadings (default) lasso2 $model, l(`lam') alpha(`ai') mat A = e(beta) // pre-estimation standardisation of data lasso2 $model, l(`lam') prestd alpha(`ai') mat B = e(beta) comparemat A B , tol(10e-6) } } * // elastic net [nocons] foreach lam of numlist 1 10 50 150 160 { foreach ai of numlist 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 { // in original units // standardisation using penalty loadings (default) lasso2 $model, l(`lam') alpha(`ai') nocons mat A = e(beta) // pre-estimation standardisation of data lasso2 $model, l(`lam') prestd alpha(`ai') nocons mat B = e(beta) comparemat A B , tol(10e-6) } } * ******************************************************************************** *** verify ridge regression results *** ******************************************************************************** lasso2 $model, l(150) alpha(0) unitload mat A = e(beta) mat A = A[1,1..8] // excl intercept estridge $model, l(150) mat B = r(bhat) comparemat A B , tol(10e-6) ******************************************************************************** *** verify post-estimation OLS results *** ******************************************************************************** * lasso foreach i of numlist 0.5 1 4 10 15 50 150 { lasso2 $model, l(`i') ols mat A = e(b) reg lpsa `e(selected)' mat B = e(b) comparemat A B } * * lasso (w/o constant) foreach i of numlist 0.5 1 4 10 15 50 150 { lasso2 $model, l(`i') ols nocons mat A = e(b) reg lpsa `e(selected)', nocons mat B = e(b) comparemat A B } * * sqrt lasso foreach i of numlist 0.5 1 4 10 15 50 150 { lasso2 $model, l(`i') ols sqrt mat A = e(b) reg lpsa `e(selected)' mat B = e(b) comparemat A B } * * elastic net foreach i of numlist 0.5 1 4 10 15 50 150 { lasso2 $model, l(`i') ols alpha(.5) mat A = e(b) reg lpsa `e(selected)' mat B = e(b) comparemat A B } * * ridge foreach i of numlist 0.5 1 4 10 15 50 150 { lasso2 $model, l(`i') ols alpha(0) mat A = e(b) reg lpsa `e(selected)' mat B = e(b) comparemat A B } * ******************************************************************************** *** partial() vs notpen() *** ******************************************************************************** lasso2 $model, partial(lcp) l(50) mat A = e(b) lasso2 $model, notpen(lcp) l(50) mat B = e(b) comparevec A B lasso2 $model, partial(lcp) l(50) sqrt mat A = e(b) lasso2 $model, notpen(lcp) l(50) sqrt mat B = e(b) comparevec A B lasso2 $model, partial(lcp) l(50) alpha(0.5) mat A = e(b) lasso2 $model, notpen(lcp) l(50) alpha(0.5) mat B = e(b) comparevec A B lasso2 $model, partial(lcp) l(50) alpha(0) mat A = e(b) lasso2 $model, notpen(lcp) l(50) alpha(0) mat B = e(b) comparevec A B lasso2 $model, lambda(10) partial(age) notpen(svi) mat A = e(b) lasso2 $model, lambda(10) partial(svi) notpen(age) mat B = e(b) comparevec A B ******************************************************************************** *** penalty loadings vs notpen (see help file) *** ******************************************************************************** lasso2 $model, l(10) notpen(lcavol) unitloadings mat A = e(b) mat myloadings = (0,1,1,1,1,1,1,1) lasso2 $model, l(10) ploadings(myloadings) mat B = e(b) comparemat A B ******************************************************************************** *** ic option to control display of output *** ******************************************************************************** di as red "should display EBIC (the default):" lasso2 $model sleep 1000 di as red "should display AIC:" lasso2 $model , ic(aic) sleep 1000 di as red "should display AICc:" lasso2 $model , ic(aicc) sleep 1000 di as red "should display BIC:" lasso2 $model , ic(bic) sleep 1000 di as red "should display EBIC:" lasso2 $model , ic(ebic) ******************************************************************************** *** degrees of freedom calculation *** ******************************************************************************** // replicate dof w/o constant and no standardisation [OK] lasso2 $model , alpha(0) l(20 .1) long unitload nocons nopath mat D= e(dof) mat list e(dof) putmata y=(lpsa) X=(lcavol lweight age lbph svi lcp gleason pgg45 ), replace mata: df=trace(X*invsym((X'X):+20/2*I(8))*X') mata: df mata: st_local("df",strofreal(df)) assert reldif(el(D,1,1),`df')<10^-6 // standardisation w/o constant [OK] lasso2 $model , alpha(0) l(20 .1) long nocons nopath mat D1 = e(dof) mat list e(dof) lasso2 $model , alpha(0) l(20 .1) long prestd nocons nopath mat D2 = e(dof) mat list e(dof) comparemat D1 D2 putmata y=(lpsa) X=(lcavol lweight age lbph svi lcp gleason pgg45 ), replace mata: s = sqrt(mean((X:-mean(X)):^2)) mata: ssq = s:^2 mata: Xs=X:/s mata: df1=trace(X*invsym((X'X):+20/2*diag(ssq))*X') // "on the fly" standardisation mata: df2=trace(Xs*invsym((Xs'Xs):+20/2*I(8))*Xs') // pre-standardisation mata: df1,df2 mata: st_local("df1",strofreal(df1)) mata: st_local("df2",strofreal(df2)) assert reldif(el(D1,1,1),`df1')<10^-6 assert reldif(el(D1,1,1),`df2')<10^-6 assert reldif(`df1',`df2')<10^-6 // dof with constant and no standardisation [OK] lasso2 $model , alpha(0) l(20 .1) long unitload nopath mat list e(dof) mat D = e(dof) putmata y=(lpsa) X=(lcavol lweight age lbph svi lcp gleason pgg45 ), replace mata: Xone=(X,J(97,1,1)) mata: Psicons = I(9) mata: Psicons[9,9]=0 mata: Xdm = X :- mean(X) mata: trace(X*invsym((X'X):+20/2*I(8))*X') // this is wrong (ignores constant) mata: trace(Xdm*invsym((Xdm'Xdm):+20/2*I(8))*Xdm') // this is correct but missing the constant mata: trace(Xone*invsym((Xone'Xone):+20/2*Psicons)*Xone') // this is correct (incl constant) mata: df1=trace(Xone*invsym((Xone'Xone):+20/2*Psicons)*Xone') // this should be correct mata: df2=trace(Xdm*invsym((Xdm'Xdm):+20/2*I(8))*Xdm')+1 // this is correct after +1 for constant mata: df1,df2 mata: st_local("df1",strofreal(df1)) mata: st_local("df2",strofreal(df2)) assert reldif(el(D,1,1),`df1')<10^-6 assert reldif(el(D,1,1),`df2')<10^-6 assert reldif(`df1',`df2')<10^-6 // standardisation w/ constant [OK] lasso2 $model , alpha(0) l(20 10 1 .1) long nopath mat list e(dof) mat D1 = e(dof) lasso2 $model , alpha(0) l(20 10 1 .1) long prestd nopath mat list e(dof) mat D2 = e(dof) comparemat D1 D2 putmata y=(lpsa) X=(lcavol lweight age lbph svi lcp gleason pgg45), replace mata: s = sqrt(mean((X:-mean(X)):^2)) mata: ssq = s:^2 mata: Xs=(X:-mean(X)):/s mata: Xone=(X,J(97,1,1)) mata: df1=trace(Xdm*invsym((Xdm'Xdm):+20/2*diag(ssq))*Xdm') +1 mata: df2=trace(Xs*invsym((Xs'Xs):+20/2*I(8))*Xs') +1 mata: df1,df2 mata: st_local("df1",strofreal(df1)) mata: st_local("df2",strofreal(df2)) assert reldif(el(D1,1,1),`df1')<10^-6 assert reldif(el(D1,1,1),`df2')<10^-6 assert reldif(`df1',`df2')<10^-6 ******************************************************************************** *** lic option *** ******************************************************************************** * check that right lambda is used foreach ic of newlist ebic aic aicc bic { lasso2 $model local optlambda=e(l`ic') lasso2, lic(`ic') postres local thislambda=e(lambda) assert reldif(`optlambda',`thislambda')<10^-8 } * ******************************************************************************** *** sklearn option *** ******************************************************************************** cap python query // run this section only if Stata is Python-aware if _rc==0 { sysuse auto, clear // single lambda // lasso lasso2 mpg rep78-foreign, lglmnet lambda(1) savedresults save nosklearn e() lasso2 mpg rep78-foreign, lglmnet lambda(1) sklearn savedresults comp nosklearn e(), /// exclude(macros: lasso2opt scalars: niter) /// tol(1e-7) verbose // ridge lasso2 mpg rep78-foreign, lglmnet lambda(1) alpha(0) savedresults save nosklearn e() lasso2 mpg rep78-foreign, lglmnet lambda(1) alpha(0) sklearn savedresults comp nosklearn e(), /// exclude(macros: lasso2opt scalars: niter) /// tol(1e-7) verbose // elastic net lasso2 mpg rep78-foreign, lglmnet lambda(1) alpha(0.5) savedresults save nosklearn e() lasso2 mpg rep78-foreign, lglmnet lambda(1) alpha(0.5) sklearn savedresults comp nosklearn e(), /// exclude(macros: lasso2opt scalars: niter) /// tol(1e-7) verbose // unitloadings lasso2 mpg rep78-foreign, lglmnet lambda(1) unitloadings savedresults save nosklearn e() lasso2 mpg rep78-foreign, lglmnet lambda(1) unitloadings sklearn savedresults comp nosklearn e(), /// exclude(macros: lasso2opt scalars: niter) /// tol(1e-7) verbose // lambda path // lasso // nodevcrit required because sklearn does full lambda grid lasso2 mpg rep78-foreign, lglmnet long nodevcrit savedresults save nosklearn e() lasso2 mpg rep78-foreign, lglmnet long sklearn savedresults comp nosklearn e(), /// exclude(macros: lasso2opt scalars: lcount) /// tol(1e-7) verbose // ridge // nodevcrit required because sklearn does full lambda grid lasso2 mpg rep78-foreign, lglmnet long alpha(0) nodevcrit savedresults save nosklearn e() lasso2 mpg rep78-foreign, lglmnet long alpha(0) sklearn savedresults comp nosklearn e(), /// exclude(macros: lasso2opt scalars: lcount) /// tol(1e-7) verbose // elastic net // nodevcrit required because sklearn does full lambda grid lasso2 mpg rep78-foreign, lglmnet long alpha(0.5) nodevcrit savedresults save nosklearn e() lasso2 mpg rep78-foreign, lglmnet long alpha(0.5) sklearn savedresults comp nosklearn e(), /// exclude(macros: lasso2opt scalars: lcount) /// tol(1e-7) verbose // unitloadings // nodevcrit required because sklearn does full lambda grid lasso2 mpg rep78-foreign, lglmnet long unitloadings nodevcrit savedresults save nosklearn e() lasso2 mpg rep78-foreign, lglmnet long sklearn unitloadings savedresults comp nosklearn e(), /// exclude(macros: lasso2opt scalars: lcount) /// tol(1e-7) verbose // error checking // sklearn not supported with custom penalty loadings mat psi = 4, 3, 2, 1, 1, 1, 1, 2, 3 cap lasso2 mpg rep78-foreign, lambda(1) lglmnet sklearn ploadings(psi) assert _rc > 0 } ******************************************************************************** *** predicted values (see help file) *** ******************************************************************************** // xbhat1 is generated by re-estimating the model for lambda=10. 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. * load example data insheet using "$prostate", tab clear global model lpsa lcavol lweight age lbph svi lcp gleason pgg45 lasso2 $model cap drop xbhat1 predict double xbhat1, xb l(10) noisily cap drop xbhat2 predict double xbhat2, xb l(10) approx // The model is estimated explicitly using lambda=100. If lasso2 is //called with a scalar lambda value, the // subsequent predict command requires no lambda() option. lasso2 $model, 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. assert (xbhat1-xbhat2<10e-8) & (xbhat3-xbhat2<10e-8) //It is also possible to obtain predicted values by referencing a specific // lambda ID using the lid() option. lasso2 $model cap drop xbhat4 predict double xbhat4, xb lid(21) cap drop xbhat5 predict double xbhat5, xb l(25.45473900468241) assert (xbhat4-xbhat5<10e-8) ******************************************************************************** *** misc options/syntax checks *** ******************************************************************************** // Support for inrange(.) and similar [if] expressions: lasso2 $model if inrange(age,50,70) ******************************************************************************** *** plotting *** ******************************************************************************** lasso2 $model lasso2, plotpath(lambda) plotlabel plotopt(legend(off)) lasso2, plotpath(lnlambda) plotlabel plotopt(legend(off)) lasso2, plotpath(norm) plotlabel plotopt(legend(off)) lasso2, plotpath(norm) plotlabel plotopt(legend(off)) plotvar(lcavol) ******************************************************************************** *** validate RSS / r-squared *** ******************************************************************************** // loop over three mehods. "nopath" corresponds to default standardisation on the fly. // "nopath" is just a placeholder that doesn't affect calculations. foreach method in prestd unitl nopath { foreach a of numlist 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 { lasso2 $model, alpha(`a') l(10) `method' cap drop xb cap drop r predict double xb, xb predict double r, r sum lpsa di r(sd)^2*98 local tss=r(sd)^2*98 sum xb di r(sd)^2*98 local ess=r(sd)^2*98 sum r di r(sd)^2*98 local rss=r(sd)^2*98 di "r-squared" local rsq = 1-`rss'/`tss' di `rsq' lasso2 $model, alpha(`a') l(10 0.1) `method' mat list e(rsq) mat RSQ = e(rsq) di el(RSQ,1,1) di reldif(`rsq',el(RSQ,1,1)) assert reldif(`rsq',el(RSQ,1,1))<10^(-3) } } * ******************************************************************************** *** validate EBIC default gamma *** ******************************************************************************** webuse air2, clear lasso2 air L(1/24).air local myebicgamma = 1-log(e(N))/(2*log(e(p))) di `myebicgamma' di e(ebicgamma) assert reldif(`myebicgamma',e(ebicgamma))<10^(-3) ******************************************************************************** *** panel example: validate within transformation *** ******************************************************************************** use "http://fmwww.bc.edu/ec-p/data/macro/abdata.dta", clear lasso2 ys l(0/3).k l(0/3).n, fe l(10) mat A = e(b) lasso2 ys l(0/3).k l(0/3).n ibn.id, partial(ibn.id) l(10) nor mat B = e(b) comparemat A B // noftools option lasso2 ys l(0/3).k l(0/3).n, fe l(10) savedresults save ftools e() cap noi assert "`e(noftools)'"=="" // will be error if ftools not installed lasso2 ys l(0/3).k l(0/3).n, fe l(10) noftools assert "`e(noftools)'"=="noftools" savedresults comp ftools e(), exclude(macros: lasso2opt) ******************************************************************************** *** check if partial() works with fe *** *** and followed by lic() *** ******************************************************************************** clear use https://www.stata-press.com/data/r16/nlswork lasso2 ln_w grade age c.age#c.age ttl_exp c.ttl_exp#c.ttl_exp tenure /// c.tenure#c.tenure 2.race not_smsa south i.year, fe ereturn list lasso2 ln_w grade age c.age#c.age ttl_exp c.ttl_exp#c.ttl_exp tenure /// c.tenure#c.tenure 2.race not_smsa south i.year, fe partial(i.year) lasso2, lic(ebic) lasso2 ln_w grade age c.age#c.age ttl_exp c.ttl_exp#c.ttl_exp tenure /// c.tenure#c.tenure 2.race not_smsa south i.year, fe partial(i.year) /// lic(ebic) ******************************************************************************** *** check residuals with fe *** ******************************************************************************** clear use https://www.stata-press.com/data/r16/nlswork foreach opt in xb u e ue xbu { cap drop `opt'hat_lasso2 cap drop `opt'hat_xtreg } cap drop esample **replace ln_w = . if year == 80 lasso2 ln_w grade age c.age#c.age ttl_exp c.ttl_exp#c.ttl_exp tenure /// c.tenure#c.tenure 2.race not_smsa south , fe gen byte esample=e(sample) lasso2, lic(ebic) postres ols mat bl2 = e(b) local selected = e(selected) di "`selected'" //ereturn list // confirm post-lasso LSDV matches xtreg,fe foreach opt in xb u e ue xbu { predict double `opt'hat_lasso2, ols `opt' } xtreg ln_w `selected' if esample, fe foreach opt in xb u e ue xbu { predict double `opt'hat_xtreg, `opt' } foreach opt in xb u e ue xbu { di "checking option `opt' assert reldif(`opt'hat_lasso2, `opt'hat_xtreg) < 1e-6 } ** alt approach from 2nd half clear use https://www.stata-press.com/data/r16/nlswork replace ln_w = . if year == 80 lasso2 ln_w grade age c.age#c.age ttl_exp c.ttl_exp#c.ttl_exp tenure /// c.tenure#c.tenure 2.race not_smsa south , fe lasso2, lic(ebic) postres ols local sel = e(selected) di "`sel'" predict double uehat , ue noi predict double ehat , e noi predict double xbhat , xb noi predict double xbuhat , xbu noi predict double uhat , u noi xtreg ln_w `sel' if e(sample), fe mat bxtreg = e(b) predict double uehat_xtreg , ue predict double ehat_xtreg , e predict double xbhat_xtreg , xb predict double xbuhat_xtreg , xbu predict double uhat_xtreg , u assert abs(ehat_xtreg-ehat)<10e-8 | (missing(ehat_xtreg) | missing(ehat)) assert abs(uehat_xtreg-uehat)<10e-8 | (missing(uehat_xtreg) | missing(uehat)) assert abs(xbhat_xtreg-xbhat)<10e-8 | (missing(xbhat_xtreg) | missing(xbhat)) assert abs(xbuhat_xtreg-xbuhat)<10e-8 | (missing(xbuhat_xtreg) | missing(xbuhat)) assert abs(uhat_xtreg-uhat)<10e-8 | (missing(uhat_xtreg) | missing(uhat)) ******************************************************************************** *** finish *** ******************************************************************************** cap log close set rmsg off