1.6.13.1 The data🔗ℹ

glmnet/datasets provides R’s mtcars: 32 cars from the 1974 Motor Trend road tests, among them their miles per gallon ("mpg"), weight in thousands of pounds ("wt") and number of cylinders ("cyl"), 4, 6 or 8. It also provides R’s iris: four measurements of 150 flowers and their species, as strings:

> (require glmnet/datasets)

> (table-column-names mtcars)

'("mpg" "cyl" "disp" "hp" "drat" "wt" "qsec" "vs" "am" "gear" "carb")

> (table-column-names iris)

'("Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width" "Species")

1.6.13.2 A factor🔗ℹ

(factor cyl) codes the cylinder counts by treatment contrasts: four cylinders, the first level, is the baseline, and six and eight get a 0/1 column each:

> (define X (formula-design-matrix (mpg . ~ . wt + (factor cyl)) mtcars))

> (design-matrix-column-names X)

'("wt" "(factor cyl)6" "(factor cyl)8")

> (design-matrix->rows (design-matrix-select-rows X '(0 2 4)))

'((2.62 1.0 0.0) (2.32 0.0 0.0) (3.44 0.0 1.0))

The first, third and fifth cars have six, four and eight cylinders. formula-path fits the lasso path:

> (define cylinders     (formula-path (mpg . ~ . wt + (factor cyl)) mtcars))

> (formula-model-levels cylinders)

'(("(factor cyl)" "4" "6" "8"))

> (coef cylinders #:lambda 1.0)

'(("(Intercept)" . 32.60819854761173)   ("wt" . -3.6480225224836307)   ("(factor cyl)6" . 0.0)   ("(factor cyl)8" . -1.785081342060045))

> (coef cylinders #:lambda 0.1)

'(("(Intercept)" . 33.87817654234954)   ("wt" . -3.2723606946002666)   ("(factor cyl)6" . -3.744964749774043)   ("(factor cyl)8" . -5.577915562591466))

> (coef cylinders #:lambda 0.01)

'(("(Intercept)" . 33.973711143650725)   ("wt" . -3.2208608052011565)   ("(factor cyl)6" . -4.1533045215146736)   ("(factor cyl)8" . -5.970825952082223))

The weight enters the path first, alone. Eight cylinders follow at λ ≈ 2.45, since the eight-cylinder cars differ most from the four-cylinder ones, and six cylinders only at λ ≈ 0.80. At λ = 0.01 the coefficients are close to those of R’s unpenalized lm(mpg ~ wt + factor(cyl), mtcars), 34.0 − 3.21 wt − 4.26 for six cylinders and − 6.07 for eight: at the same weight, six cylinders cost about 4 miles per gallon and eight about 6.

1.6.13.3 A slope for each cylinder count🔗ℹ

(* wt (factor cyl)) crosses the weight with the cylinders. The interactions are the weight in the rows of six or eight cylinders and 0 in the others, so each is a difference between that level’s slope and the four-cylinder cars’ slope:

> (define slopes (formula-path (mpg . ~ . wt * (factor cyl)) mtcars))

> (formula-model-predictor-names slopes)

'("wt" "(factor cyl)6" "(factor cyl)8" "wt:(factor cyl)6" "wt:(factor cyl)8")

> (coef slopes #:lambda 0.01)

'(("(Intercept)" . 38.019185835953586)   ("wt" . -4.9806957794227555)   ("(factor cyl)6" . -6.024488400937778)   ("(factor cyl)8" . -13.398640376657047)   ("wt:(factor cyl)6" . 1.0564291290632093)   ("wt:(factor cyl)8" . 2.603363975356239))

The interactions enter late, eight cylinders at λ ≈ 0.059 and six at λ ≈ 0.016. At λ = 0.01, a thousand pounds cost a four-cylinder car 5.0 miles per gallon, a six-cylinder car 4.98 − 1.06 = 3.9 and an eight-cylinder car 4.98 − 2.60 = 2.4. R’s lm puts the slopes further apart, at 5.6, 2.8 and 2.2, and the lasso has shrunk the differences.

1.6.13.4 New cars🔗ℹ

predict codes new cars with the levels the model was fitted with. These three have only four and eight cylinders, which is fine; a car with five is an error that names the level, as R’s "factor has new levels" is:

> (define new-cars (list (cons "wt" '(2.5 3.5 3.5)) (cons "cyl" '(4 4 8))))

> (predict cylinders new-cars #:lambda 0.1)

'(25.697274805848878 22.42491411124861 16.846998548657147)

> (predict cylinders (list (cons "wt" '(3.0)) (cons "cyl" '(5))))

predict: a factor has levels that the model was not fitted

with

factor: "(factor cyl)"

new levels: '("5")

levels: '("4" "6" "8")

At λ = 0.1, a thousand pounds more cost the four-cylinder car 3.3 miles per gallon, and eight cylinders instead of four another 5.6.

1.6.13.5 Species as classes🔗ℹ

The species are strings, so a multinomial fit on them has the species as classes, in R’s order of levels. coef keys the coefficients by species, and predict with #:type 'class names a species for each new flower:

> (define species     (formula-path (~ Species all) iris #:family 'multinomial))

> (coef species #:lambda 0.05)

'(("setosa"    ("(Intercept)" . 2.869706551251394)    ("Sepal.Length" . 0.0)    ("Sepal.Width" . 0.7350548897518949)    ("Petal.Length" . -1.3604056950808465)    ("Petal.Width" . 0.0))   ("versicolor"    ("(Intercept)" . 1.3944568206526933)    ("Sepal.Length" . 0.0)    ("Sepal.Width" . -0.06639471192067056)    ("Petal.Length" . 0.0)    ("Petal.Width" . 0.0))   ("virginica"    ("(Intercept)" . -4.264163371904088)    ("Sepal.Length" . 0.0)    ("Sepal.Width" . 0.0)    ("Petal.Length" . 0.0)    ("Petal.Width" . 3.329533653482147)))

> (define flowers     (list (cons "Sepal.Length" '(5.0 6.0 6.5))           (cons "Sepal.Width" '(3.4 2.8 3.0))           (cons "Petal.Length" '(1.5 4.5 5.6))           (cons "Petal.Width" '(0.2 1.4 2.1))))

> (predict species flowers #:type 'class #:lambda 0.05)

'("setosa" "versicolor" "virginica")

> (predict species flowers #:type 'response #:lambda 0.05)

'((0.8957714953832591 0.10334940400127686 0.0008791006154639557)   (0.058961068224158836 0.6515762179126585 0.2894627138631826)   (0.004206932896646122 0.17686862244732146 0.8189244446560323))

At λ = 0.05 each species keeps few measurements: short petals and wide sepals for setosa, narrow sepals for versicolor, and wide petals for virginica. The first flower is a setosa with probability 0.90, the second a versicolor with 0.65 and the third a virginica with 0.82.

1.6.13.6 Matching R🔗ℹ

The same fits in R are

and the same for mpg ~ wt * factor(cyl). R 4.5.3 with glmnet 4.1.10 gives the same paths, 61 λ values from 5.147 down to 0.01938 for the cylinders, 92 down to 0.001083 for the slopes and 100 down to 0.0000435 for the species, and the coefficients, predictions and classes above. R names the columns factor(cyl)6 and wt:factor(cyl)6, which are "(factor cyl)6" and "wt:(factor cyl)6" here. The example’s companion test checks each number against R’s to six significant digits.