Lecture 28: Stepwise and Penalised Regression

161.251 Regression Modelling

From manual model building to automated search

In Lecture 27, we still chose the candidate models ourselves. Once the number of plausible predictors grows, there are two broad ways to stop doing all of that by hand.

We will use both approaches in this lecture.

Forward, backward, and stepwise search

The basic idea is simple.

• Forward selection starts small and adds one term at a time.
• Backward elimination starts large and removes one term at a time.
• Stepwise selection allows moves in either direction.

. . .

The unit being added or removed is a model term, not an arbitrary coefficient. A factor can generate several coefficients and still count as one term.

For a prespecified nested comparison, a partial \(F\) test is still useful. For a larger automated search, we need a criterion that can be evaluated repeatedly. Here we will use AIC, introduced in Lecture 27.

In practice, step() looks at models one move away from the current model and takes a move that lowers AIC. That makes the search local, so the path and starting model can matter.

Climate data: small enough to watch the search

The climate data are small enough that we can actually watch the search process move from one model to the next.

climate <- read.csv("../data/climate.csv", header = TRUE, row.names = 1) |>
    drop_na(Sun, Lat, Long, MnJanTemp, MnJlyTemp, Rain, Height, Sea, NorthIsland)

climate_null <- lm(Sun ~ 1, data = climate)
climate_full <- lm(Sun ~ Lat + Long + MnJanTemp + MnJlyTemp + Rain + Height + Sea +
    NorthIsland, data = climate)
climate_scope <- list(lower = formula(climate_null), upper = formula(climate_full))

. . .

The full model for annual sunshine is

. . .

\[ \text{Sun} \sim \text{Lat} + \text{Long} + \text{MnJanTemp} + \text{MnJlyTemp} + \text{Rain} + \text{Height} + \text{Sea} + \text{NorthIsland}. \]

Forward selection

Start with the intercept-only model. At each step, ask which single addition gives the best reduction in AIC. We will follow each step and examine how the candidate models are considered by the algorithm.

climate_forward <- step(climate_null, scope = climate_scope, direction = "forward",
    trace = 1)

Start:  AIC=388.08
Sun ~ 1

              Df Sum of Sq     RSS    AIC
+ MnJanTemp    1    506057 1130350 376.76
+ Long         1    344588 1291820 381.57
+ Rain         1    316258 1320150 382.35
+ Lat          1    220053 1416355 384.88
+ MnJlyTemp    1    182734 1453674 385.82
+ Sea          1    148252 1488156 386.66
+ NorthIsland  1     95788 1540620 387.91
+ Height       1     90941 1545467 388.02
<none>                     1636408 388.08

Step:  AIC=376.76
Sun ~ MnJanTemp

              Df Sum of Sq     RSS    AIC
+ Rain         1    132170  998181 374.29
+ NorthIsland  1     76944 1053407 376.22
<none>                     1130350 376.76
+ Sea          1     50284 1080066 377.12
+ Lat          1     43712 1086639 377.34
+ MnJlyTemp    1     22898 1107452 378.03
+ Long         1     14102 1116248 378.31
+ Height       1        40 1130310 378.76

Step:  AIC=374.29
Sun ~ MnJanTemp + Rain

              Df Sum of Sq    RSS    AIC
+ Sea          1     87195 910985 373.00
<none>                     998181 374.29
+ NorthIsland  1     53140 945040 374.32
+ Long         1     12081 986099 375.85
+ Height       1      1193 996988 376.24
+ Lat          1       330 997850 376.27
+ MnJlyTemp    1        16 998165 376.29

Step:  AIC=373
Sun ~ MnJanTemp + Rain + Sea

              Df Sum of Sq    RSS    AIC
+ Height       1     94788 816197 371.04
+ MnJlyTemp    1     52902 858083 372.84
<none>                     910985 373.00
+ Long         1     38986 871999 373.42
+ NorthIsland  1     16393 894592 374.34
+ Lat          1      1945 909040 374.92

Step:  AIC=371.04
Sun ~ MnJanTemp + Rain + Sea + Height

              Df Sum of Sq    RSS    AIC
<none>                     816197 371.04
+ NorthIsland  1     28556 787641 371.76
+ MnJlyTemp    1     22773 793424 372.02
+ Long         1     19420 796777 372.17
+ Lat          1     10788 805409 372.56

. . .

Finally stopped.Let’s check what we ended up with again.

. . .

formula(climate_forward)

Sun ~ MnJanTemp + Rain + Sea + Height

AIC(climate_forward)

[1] 475.2037

Backward elimination

This is the same overall idea as forward selection, however we start with all terms in the model and kick terms out based on the best reduction in AIC.

climate_backward <- step(climate_full, scope = climate_scope, direction = "backward",
    trace = 1)

Start:  AIC=372.48
Sun ~ Lat + Long + MnJanTemp + MnJlyTemp + Rain + Height + Sea + 
    NorthIsland

              Df Sum of Sq    RSS    AIC
- Lat          1       760 680941 370.52
- MnJlyTemp    1      9027 689208 370.95
- Height       1     21370 701551 371.59
<none>                     680181 372.48
- NorthIsland  1     47272 727453 372.90
- MnJanTemp    1     62239 742420 373.63
- Rain         1     79398 759579 374.45
- Long         1    104127 784308 375.61
- Sea          1    112082 792263 375.97

Step:  AIC=370.52
Sun ~ Long + MnJanTemp + MnJlyTemp + Rain + Height + Sea + NorthIsland

              Df Sum of Sq    RSS    AIC
- MnJlyTemp    1      8381 689322 368.96
<none>                     680941 370.52
- Height       1     39668 720609 370.56
- NorthIsland  1     48846 729788 371.01
- Rain         1     86487 767428 372.82
- Long         1    104518 785460 373.66
- Sea          1    114410 795351 374.11
- MnJanTemp    1    147165 828106 375.56

Step:  AIC=368.96
Sun ~ Long + MnJanTemp + Rain + Height + Sea + NorthIsland

              Df Sum of Sq    RSS    AIC
<none>                     689322 368.96
- Height       1     76044 765366 370.73
- Long         1     98319 787641 371.76
- NorthIsland  1    107455 796777 372.17
- Sea          1    125566 814888 372.98
- MnJanTemp    1    140484 829807 373.64
- Rain         1    142964 832286 373.74

. . .

Done.

formula(climate_backward)

Sun ~ Long + MnJanTemp + Rain + Height + Sea + NorthIsland

AIC(climate_backward)

[1] 473.1216

. . .

Now where did we end up? Not the same model as the forward selection. Now let’s try stepwise!?

Stepwise selection

Stepwise selection allows both kinds of move. We start small, but after adding a term the algorithm can reconsider terms already in the model.

climate_both <- step(climate_null, scope = climate_scope, direction = "both", trace = 1)

Start:  AIC=388.08
Sun ~ 1

              Df Sum of Sq     RSS    AIC
+ MnJanTemp    1    506057 1130350 376.76
+ Long         1    344588 1291820 381.57
+ Rain         1    316258 1320150 382.35
+ Lat          1    220053 1416355 384.88
+ MnJlyTemp    1    182734 1453674 385.82
+ Sea          1    148252 1488156 386.66
+ NorthIsland  1     95788 1540620 387.91
+ Height       1     90941 1545467 388.02
<none>                     1636408 388.08

Step:  AIC=376.76
Sun ~ MnJanTemp

              Df Sum of Sq     RSS    AIC
+ Rain         1    132170  998181 374.29
+ NorthIsland  1     76944 1053407 376.22
<none>                     1130350 376.76
+ Sea          1     50284 1080066 377.12
+ Lat          1     43712 1086639 377.34
+ MnJlyTemp    1     22898 1107452 378.03
+ Long         1     14102 1116248 378.31
+ Height       1        40 1130310 378.76
- MnJanTemp    1    506057 1636408 388.08

Step:  AIC=374.29
Sun ~ MnJanTemp + Rain

              Df Sum of Sq     RSS    AIC
+ Sea          1     87195  910985 373.00
<none>                      998181 374.29
+ NorthIsland  1     53140  945040 374.32
+ Long         1     12081  986099 375.85
+ Height       1      1193  996988 376.24
+ Lat          1       330  997850 376.27
+ MnJlyTemp    1        16  998165 376.29
- Rain         1    132170 1130350 376.76
- MnJanTemp    1    321969 1320150 382.35

Step:  AIC=373
Sun ~ MnJanTemp + Rain + Sea

              Df Sum of Sq     RSS    AIC
+ Height       1     94788  816197 371.04
+ MnJlyTemp    1     52902  858083 372.84
<none>                      910985 373.00
+ Long         1     38986  871999 373.42
- Sea          1     87195  998181 374.29
+ NorthIsland  1     16393  894592 374.34
+ Lat          1      1945  909040 374.92
- Rain         1    169081 1080066 377.12
- MnJanTemp    1    214016 1125001 378.59

Step:  AIC=371.04
Sun ~ MnJanTemp + Rain + Sea + Height

              Df Sum of Sq     RSS    AIC
<none>                      816197 371.04
+ NorthIsland  1     28556  787641 371.76
+ MnJlyTemp    1     22773  793424 372.02
+ Long         1     19420  796777 372.17
+ Lat          1     10788  805409 372.56
- Height       1     94788  910985 373.00
- Sea          1    180791  996988 376.24
- Rain         1    226572 1042769 377.86
- MnJanTemp    1    284136 1100333 379.79

. . .

Done.

. . .

formula(climate_both)

Sun ~ MnJanTemp + Rain + Sea + Height

AIC(climate_both)

[1] 475.2037

Does the bidirectional search finish at the forward model, the backward model, or somewhere else?

Comparing the climate searches

Search	Final formula	AIC
Forward	Sun ~ MnJanTemp + Rain + Sea + Height	475.20
Backward	Sun ~ Long + MnJanTemp + Rain + Height + Sea + NorthIsland	473.12
Stepwise: both	Sun ~ MnJanTemp + Rain + Sea + Height	475.20

The lower AIC is preferred among the models actually reached by the search. Different starting points can finish at different local solutions. Agreement between searches is reassuring for this dataset, but it does not make the selected formula uniquely correct.

A larger problem: Seoul bike demand

The Seoul data contain hourly bike rentals together with weather and calendar variables. We have already watched forward, backward and bidirectional search on the smaller climate problem, so here we will use one bidirectional search and focus on defining the model space.

seoul_hourly <- read_csv("../data/seoul_bike_hourly.csv", show_col_types = FALSE) |>
    mutate(hour = factor(hour))

seoul_subset <- seoul_hourly |>
    select(rented_bike_count, hour, season, holiday, functioning_day, temperature,
        humidity, wind_speed, visibility, dew_point_temperature, solar_radiation,
        rainfall, snowfall) |>
    drop_na()

Note that we are treating hour records as a factor rather than forcing bike demand to change linearly from midnight to 11 pm.

. . .

seoul_null <- lm(rented_bike_count ~ 1, data = seoul_subset)

seoul_full <- lm(rented_bike_count ~ hour + season + holiday + functioning_day +
    temperature + humidity + wind_speed + visibility + dew_point_temperature + solar_radiation +
    rainfall + snowfall, data = seoul_subset)

seoul_scope <- list(lower = formula(seoul_null), upper = formula(seoul_full))

seoul_scope

$lower
rented_bike_count ~ 1

$upper
rented_bike_count ~ hour + season + holiday + functioning_day + 
    temperature + humidity + wind_speed + visibility + dew_point_temperature + 
    solar_radiation + rainfall + snowfall

The R term scope tells step() the smallest and largest formulas it is allowed to consider. This is handy if you always want a certain term in the model.

Stepwise search on the Seoul data

seoul_step <- step(seoul_null, scope = seoul_scope, direction = "both", trace = 1)

Start:  AIC=113342.2
rented_bike_count ~ 1

                        Df  Sum of Sq        RSS    AIC
+ temperature            1 1056904520 2587029843 110343
+ hour                  23 1064201100 2579733263 110363
+ season                 3  765709006 2878225356 111282
+ dew_point_temperature  1  525597511 3118336852 111980
+ solar_radiation        1  249823064 3394111299 112722
+ functioning_day        1  151560828 3492373534 112972
+ humidity               1  145437128 3498497235 112987
+ visibility             1  144710241 3499224121 112989
+ snowfall               1   73273235 3570661128 113166
+ rainfall               1   55195401 3588738962 113210
+ wind_speed             1   53446519 3590487844 113215
+ holiday                1   19067701 3624866661 113298
<none>                                3643934363 113342

Step:  AIC=110343.4
rented_bike_count ~ temperature

                        Df  Sum of Sq        RSS    AIC
+ hour                  23  843729089 1743300754 106932
+ humidity               1  304994771 2282035071 109247
+ dew_point_temperature  1  273094201 2313935642 109368
+ functioning_day        1  194870595 2392159247 109659
+ visibility             1  118918824 2468111018 109933
+ rainfall               1   82364833 2504665010 110062
+ wind_speed             1   72162334 2514867509 110098
+ season                 3   29334094 2557695748 110250
+ solar_radiation        1   21261579 2565768264 110273
+ holiday                1    6514428 2580515414 110323
+ snowfall               1    2237212 2584792630 110338
<none>                                2587029843 110343
- temperature            1 1056904520 3643934363 113342

Step:  AIC=106931.6
rented_bike_count ~ temperature + hour

                        Df Sum of Sq        RSS    AIC
+ functioning_day        1 186915633 1556385121 105940
+ humidity               1 173760455 1569540299 106014
+ dew_point_temperature  1 153236015 1590064739 106128
+ rainfall               1  90363923 1652936831 106467
+ solar_radiation        1  75740488 1667560266 106544
+ visibility             1  71178026 1672122728 106568
+ season                 3  38391455 1704909298 106743
+ holiday                1   7370701 1735930053 106896
+ snowfall               1   2460256 1740840498 106921
+ wind_speed             1   1623298 1741677456 106925
<none>                               1743300754 106932
- hour                  23 843729089 2587029843 110343
- temperature            1 836432509 2579733263 110363

Step:  AIC=105940.1
rented_bike_count ~ temperature + hour + functioning_day

                        Df Sum of Sq        RSS    AIC
+ humidity               1 171653389 1384731732 104918
+ dew_point_temperature  1 149699208 1406685913 105056
+ rainfall               1  91573998 1464811123 105411
+ season                 3  83935134 1472449987 105460
+ visibility             1  77443817 1478941304 105495
+ solar_radiation        1  75199277 1481185844 105508
+ holiday                1   5271935 1551113186 105912
+ snowfall               1   3453574 1552931547 105923
+ wind_speed             1   1795566 1554589555 105932
<none>                               1556385121 105940
- functioning_day        1 186915633 1743300754 106932
- hour                  23 835774126 2392159247 109659
- temperature            1 875980432 2432365553 109849

Step:  AIC=104918.4
rented_bike_count ~ temperature + hour + functioning_day + humidity

                        Df  Sum of Sq        RSS    AIC
+ season                 3   95589314 1289142418 104298
+ rainfall               1   39929433 1344802299 104664
+ solar_radiation        1    9463826 1375267907 104860
+ holiday                1    8093992 1376637740 104869
+ visibility             1    1888005 1382843728 104908
+ dew_point_temperature  1    1752158 1382979575 104909
+ snowfall               1     479728 1384252004 104917
<none>                                1384731732 104918
+ wind_speed             1     153460 1384578273 104919
- humidity               1  171653389 1556385121 105940
- functioning_day        1  184808566 1569540299 106014
- hour                  23  708681498 2093413230 108493
- temperature            1 1032873206 2417604939 109798

Step:  AIC=104297.8
rented_bike_count ~ temperature + hour + functioning_day + humidity + 
    season

                        Df Sum of Sq        RSS    AIC
+ rainfall               1  39434247 1249708172 104028
+ solar_radiation        1  10944435 1278197983 104225
+ holiday                1   5717923 1283424496 104261
+ dew_point_temperature  1   2599798 1286542620 104282
+ snowfall               1   1875639 1287266780 104287
<none>                               1289142418 104298
+ wind_speed             1    209094 1288933324 104298
+ visibility             1    126700 1289015718 104299
- season                 3  95589314 1384731732 104918
- temperature            1 155591977 1444734395 105294
- humidity               1 183307569 1472449987 105460
- functioning_day        1 229470047 1518612465 105731
- hour                  23 717078111 2006220529 108126

Step:  AIC=104027.6
rented_bike_count ~ temperature + hour + functioning_day + humidity + 
    season + rainfall

                        Df Sum of Sq        RSS    AIC
+ solar_radiation        1   8812374 1240895798 103968
+ holiday                1   5603476 1244104696 103990
+ snowfall               1   1362942 1248345230 104020
+ dew_point_temperature  1    936585 1248771587 104023
+ wind_speed             1    303290 1249404882 104028
<none>                               1249708172 104028
+ visibility             1     55499 1249652673 104029
- rainfall               1  39434247 1289142418 104298
- season                 3  95094127 1344802299 104664
- humidity               1 130345225 1380053397 104895
- temperature            1 149274881 1398983053 105014
- functioning_day        1 229378123 1479086294 105502
- hour                  23 730246773 1979954945 108013

Step:  AIC=103967.7
rented_bike_count ~ temperature + hour + functioning_day + humidity + 
    season + rainfall + solar_radiation

                        Df Sum of Sq        RSS    AIC
+ holiday                1   5630805 1235264993 103930
+ dew_point_temperature  1   2347099 1238548700 103953
+ snowfall               1   1054360 1239841439 103962
<none>                               1240895798 103968
+ visibility             1    172237 1240723561 103968
+ wind_speed             1      7361 1240888437 103970
- solar_radiation        1   8812374 1249708172 104028
- rainfall               1  37302185 1278197983 104225
- humidity               1  80806547 1321702345 104518
- season                 3  96310930 1337206728 104616
- temperature            1 118326884 1359222682 104764
- functioning_day        1 230699299 1471595097 105459
- hour                  23 703948892 1944844690 107858

Step:  AIC=103929.8
rented_bike_count ~ temperature + hour + functioning_day + humidity + 
    season + rainfall + solar_radiation + holiday

                        Df Sum of Sq        RSS    AIC
+ dew_point_temperature  1   2409125 1232855868 103915
+ snowfall               1    924749 1234340244 103925
<none>                               1235264993 103930
+ visibility             1    208411 1235056582 103930
+ wind_speed             1     18006 1235246987 103932
- holiday                1   5630805 1240895798 103968
- solar_radiation        1   8839703 1244104696 103990
- rainfall               1  37187683 1272452676 104188
- humidity               1  81886272 1317151266 104490
- season                 3  93917308 1329182302 104566
- temperature            1 121768802 1357033795 104751
- functioning_day        1 229113857 1464378850 105418
- hour                  23 701986875 1937251868 107826

Step:  AIC=103914.7
rented_bike_count ~ temperature + hour + functioning_day + humidity + 
    season + rainfall + solar_radiation + holiday + dew_point_temperature

                        Df Sum of Sq        RSS    AIC
+ snowfall               1   1172965 1231682903 103908
<none>                               1232855868 103915
+ visibility             1    141721 1232714147 103916
+ wind_speed             1     15752 1232840116 103917
- temperature            1   1252858 1234108726 103922
- dew_point_temperature  1   2409125 1235264993 103930
- holiday                1   5692831 1238548700 103953
- solar_radiation        1  10275359 1243131227 103985
- humidity               1  18314836 1251170704 104042
- rainfall               1  34564709 1267420577 104155
- season                 3  94736908 1327592776 104557
- functioning_day        1 230297464 1463153332 105413
- hour                  23 704315728 1937171597 107827

Step:  AIC=103908.4
rented_bike_count ~ temperature + hour + functioning_day + humidity + 
    season + rainfall + solar_radiation + holiday + dew_point_temperature + 
    snowfall

                        Df Sum of Sq        RSS    AIC
<none>                               1231682903 103908
+ visibility             1    154303 1231528601 103909
+ wind_speed             1     16316 1231666587 103910
- temperature            1   1134633 1232817536 103914
- snowfall               1   1172965 1232855868 103915
- dew_point_temperature  1   2657341 1234340244 103925
- holiday                1   5549734 1237232637 103946
- solar_radiation        1  10033817 1241716721 103977
- humidity               1  19243708 1250926611 104042
- rainfall               1  34029648 1265712551 104145
- season                 3  95724751 1327407655 104558
- functioning_day        1 229762878 1461445781 105405
- hour                  23 701524297 1933207200 107811

. . .

Done.

. . .

formula(seoul_step)

rented_bike_count ~ temperature + hour + functioning_day + humidity + 
    season + rainfall + solar_radiation + holiday + dew_point_temperature + 
    snowfall

seoul_step$anova

                      Step  Df   Deviance Resid. Df Resid. Dev      AIC
1                           NA         NA      8759 3643934363 113342.2
2            + temperature  -1 1056904520      8758 2587029843 110343.4
3                   + hour -23  843729089      8735 1743300754 106931.6
4        + functioning_day  -1  186915633      8734 1556385121 105940.1
5               + humidity  -1  171653389      8733 1384731732 104918.4
6                 + season  -3   95589314      8730 1289142418 104297.8
7               + rainfall  -1   39434247      8729 1249708172 104027.6
8        + solar_radiation  -1    8812374      8728 1240895798 103967.7
9                + holiday  -1    5630805      8727 1235264993 103929.8
10 + dew_point_temperature  -1    2409125      8726 1232855868 103914.7
11              + snowfall  -1    1172965      8725 1231682903 103908.4

This is where automated search becomes convenient. We do not want to fit every plausible formula ourselves. But the procedure still makes a consistent decision at each step: a term is either in the model or out.

Penalized Regression Models

Until now, we have taken discrete steps when selecting terms: include a term, or leave it out. Just as moving from categories to numbers gave us a smoother way to represent a variable, there is a smoother approach to building multivariable models: shrinkage.

When shrinkage is applied to regression, we call it a penalised regression model. Rather than immediately removing terms, we can keep the larger model but discourage large coefficient values. Coefficients are pulled towards zero, with the amount of shrinkage controlled by a penalty.

Different penalties let us express different kinds of model structure. We might expect the response to reflect many small influences, or instead expect that only a smaller number of predictors have substantial effects.

Remember that Ordinary Least Squares (OLS) minimises

\[ RSS = \sum_{i=1}^{n}(y_i-\hat y_i)^2. \]

Penalised regression minimises

\[ RSS + \lambda \times \text{penalty}. \]

The two penalties we will use are L2 for ridge and L1 for lasso.

• \(\lambda = 0\) gives ordinary least squares.
• Larger \(\lambda\) means a stronger penalty.
• The intercept is normally left unpenalised.

Ridge regression: L2 penalty

Ridge regression minimises

\[ RSS + \lambda\sum_{j=1}^{p}\beta_j^2. \]

• Large coefficients are penalised more strongly than small coefficients.
• As \(\lambda\) increases, the fitted coefficient vector is pulled towards zero.
• Ridge generally keeps all predictors rather than setting coefficients exactly to zero.

. . .

In matrix form,

\[ \hat{\boldsymbol\beta}_{\text{ridge}} = (X^\top X+\lambda I)^{-1}X^\top y. \]

A useful way to read this is that ridge adds \(\lambda\) to the diagonal of \(X^\top X\). That stabilises the inversion when predictors contain overlapping information.

Lasso regression: L1 penalty

Lasso minimises

\[ RSS + \lambda\sum_{j=1}^{p}|\beta_j|. \]

• Lasso also shrinks coefficients as \(\lambda\) increases.
• Unlike ridge, some coefficients can become exactly zero.
• Lasso can therefore perform shrinkage and variable selection at the same time.

. . .

Ridge and lasso both control coefficient size. The difference is what the penalty allows to happen at zero.

Why standardisation matters

The penalty acts on coefficient size, and coefficient size depends on the units used for the predictor.

Suppose the same rainfall variable is measured in millimetres or metres.

Rainfall units	Equivalent coefficient	Contribution to an L2 penalty
millimetres	0.002	4e-06
metres	2.000	4e+00

The fitted relationship is identical, but the raw coefficient penalty is not. Without standardisation, changing units can change how strongly a predictor is penalised.

Using glmnet

For a Gaussian response variable, glmnet minimises

\[ \text{Objective}(\beta_0,\boldsymbol\beta) = \frac{RSS}{2n} + \lambda\left[ \frac{1-\alpha}{2}\sum_j\beta_j^2 + \alpha\sum_j|\beta_j| \right]. \]

The first term measures fit, with RSS divided by \(2n\). The bracketed term penalises coefficient size.

Here \(\alpha\) is the mixing parameter, with \(0 \leq \alpha \leq 1\). We choose its value to set the balance between the L2 and L1 parts of the penalty.

• alpha = 0: the L1 term disappears, giving ridge (L2).
• alpha = 1: the L2 term disappears, giving lasso (L1).
• \(0 < \alpha < 1\): elastic net, a mixture of the two. We will not develop elastic net further here.

. . .

Note: glmnet standardises predictors by default when standardize = TRUE.

. . .

The remaining tuning parameter is \(\lambda\), which we choose by cross-validation.

library(glmnet)

X <- model.matrix(rented_bike_count ~ hour + season + holiday + functioning_day +
    temperature + humidity + wind_speed + visibility + dew_point_temperature + solar_radiation +
    rainfall + snowfall, data = seoul_subset)[, -1]

y <- seoul_subset$rented_bike_count

set.seed(42)
foldid <- sample(rep(seq_len(10), length.out = nrow(X)))

Ridge: choosing \(\lambda\)

lambda.min gives the lowest mean cross-validation error.

lambda.1se is the largest \(\lambda\) whose estimated error is within one standard error of that minimum. It gives stronger regularisation while remaining close to the best estimated cross-validation performance.

Note that the intercept is not penalised, but the slopes are. The penalty is applied to the coefficients of the model matrix, so it can cause the intercept to change as the slopes are shrunk.

Ridge: lambda.min versus lambda.1se

Let’s look at both the min and 1se in an example dataset.

ridge_coef_compare <- cbind(lambda.min = as.matrix(coef(seoul_ridge, s = "lambda.min"))[,
    1], lambda.1se = as.matrix(coef(seoul_ridge, s = "lambda.1se"))[, 1])

round(ridge_coef_compare, 3)

                      lambda.min lambda.1se
(Intercept)               22.494     25.913
hour1                   -127.853   -137.617
hour2                   -230.654   -234.233
hour3                   -310.709   -309.547
hour4                   -368.618   -364.096
hour5                   -357.351   -354.308
hour6                   -203.497   -211.891
hour7                     84.164     57.467
hour8                    429.401    383.211
hour9                      0.562     -6.609
hour10                  -210.398   -193.483
hour11                  -214.825   -190.558
hour12                  -175.048   -148.785
hour13                  -172.031   -143.176
hour14                  -166.132   -136.794
hour15                   -86.583    -64.576
hour16                    37.001     47.004
hour17                   292.668    278.585
hour18                   709.671    659.393
hour19                   465.680    426.805
hour20                   394.720    356.573
hour21                   385.548    345.673
hour22                   294.332    259.289
hour23                    73.516     52.597
seasonSpring            -120.309    -83.020
seasonSummer             -98.679    -44.672
seasonWinter            -342.588   -314.342
holidayNo Holiday        110.174    102.283
functioning_dayYes       859.762    773.532
temperature               14.155     12.924
humidity                  -7.285     -6.203
wind_speed                 1.095      2.601
visibility                 0.028      0.044
dew_point_temperature      7.330      6.563
solar_radiation           72.552     64.888
rainfall                 -57.972    -55.675
snowfall                  16.581      5.211

The lambda.1se column uses the stronger penalty, so the coefficient vector is more strongly regularised.

For the rest of this lecture we will use lambda.1se.

For ridge that means stronger shrinkage. For lasso it usually also means fewer non-zero coefficients.

Lasso coefficient paths

Each line is one model-matrix coefficient. The colours group related predictors; several lines can share a colour, especially the hour indicators. Moving to the right means increasing \(\lambda\).

For lasso, coefficients do not just get smaller. Some paths hit zero and stay there. At lambda.1se we use the stronger of the two marked penalties.

coef(seoul_lasso, s = "lambda.1se")

37 x 1 sparse Matrix of class "dgCMatrix"
                         lambda.1se
(Intercept)            -55.05261698
hour1                  -67.04615252
hour2                 -174.41218658
hour3                 -257.76351662
hour4                 -318.12606332
hour5                 -305.09850207
hour6                 -142.94326049
hour7                   90.73816413
hour8                  465.81340699
hour9                   33.12534552
hour10                 -94.56642398
hour11                 -84.28396219
hour12                 -33.43775151
hour13                 -30.79494042
hour14                 -33.63804606
hour15                   .         
hour16                  73.99601783
hour17                 323.42845839
hour18                 744.22563094
hour19                 478.45963889
hour20                 401.08296136
hour21                 393.26739265
hour22                 299.67977064
hour23                  68.59146328
seasonSpring           -90.28970236
seasonSummer           -57.83894556
seasonWinter          -335.82286493
holidayNo Holiday       76.37915944
functioning_dayYes     845.74464940
temperature             21.64353137
humidity                -6.52268074
wind_speed               .         
visibility               0.01082108
dew_point_temperature    .         
solar_radiation         13.14976824
rainfall               -56.58813318
snowfall                 .         

Stepwise regression versus shrinkage

Method	Decision being made	Typical result
Stepwise search	Add or remove whole model terms	A smaller selected formula
Ridge, L2	Penalise squared coefficient size	Keep predictors but shrink/stabilise the coefficient vector
Lasso, L1	Penalise absolute coefficient size	Shrink coefficients and set some exactly to zero

Stepwise and lasso can both produce sparse models, but they do not select the same object. Stepwise works with model terms. glmnet works with columns of the model matrix, so a multi-level factor can contribute several penalised coefficients. These procedures are also data-dependent. Change the sample slightly and the selected stepwise model or lasso set can change.

Looking ahead

Automated search and shrinkage give us useful ways to work with a larger predictor set, but they expose another problem.

When predictors carry nearly the same information, ordinary least-squares coefficients can become unstable. A small change in the data or model can produce a large change in the separate coefficient estimates.

That is the multicollinearity problem we turn to next.