1.6.12 Formula transforms
Source: "glmnet/examples/10-formula-polynomial.rkt"
Fuel economy falls with horsepower, but not along a straight line: an extra horsepower costs a small car more miles per gallon than a powerful one. R fits the bend with a quadratic, mpg ~ hp + I(hp^2); the formula front end writes (mpg . ~ . hp + (sqr hp)), where (sqr hp) is a transform, a Racket function of the column (see Transforms).
1.6.12.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"), horsepower ("hp") and weight in thousands of pounds ("wt"). sqr is racket/math’s:
> (require racket/math glmnet/datasets) > (table-column-names mtcars) '("mpg" "cyl" "disp" "hp" "drat" "wt" "qsec" "vs" "am" "gear" "carb")
1.6.12.2 A quadratic
(sqr hp) is a column of the design matrix, the square of each car’s horsepower, named by its source:
> (define X (formula-design-matrix (mpg . ~ . hp + (sqr hp)) mtcars)) > (design-matrix-column-names X) '("hp" "(sqr hp)")
> (design-matrix->rows (design-matrix-select-rows X '(0 1 2))) '((110.0 12100.0) (110.0 12100.0) (93.0 8649.0))
formula-path fits the lasso path on it:
> (define quadratic (formula-path (mpg . ~ . hp + (sqr hp)) mtcars)) > (coef quadratic #:lambda 1.0)
'(("(Intercept)" . 27.925167222547564)
("hp" . -0.05340974672380106)
("(sqr hp)" . 0.0))
> (coef quadratic #:lambda 0.1)
'(("(Intercept)" . 36.37237345928531)
("hp" . -0.15810816975843855)
("(sqr hp)" . 0.000265072072077869))
> (coef quadratic #:lambda 0.01)
'(("(Intercept)" . 39.974599695276595)
("hp" . -0.207377580883709)
("(sqr hp)" . 0.0004041137060060064))
Horsepower enters the path first, alone. The square follows only at λ ≈ 0.26, with a positive coefficient, which bends the line up: at λ = 0.01 the coefficients are close to those of R’s unpenalized lm(mpg ~ hp + I(hp^2), mtcars), 40.41 − 0.2133 hp + 0.000421 hp², whose lowest point is at about 253 horsepower. R’s I(hp^2) is the same column, and (I (expt hp 2)) writes it as R does:
> (define as-in-r (formula-path (mpg . ~ . hp + (I (expt hp 2))) mtcars)) > (formula-model-predictor-names as-in-r) '("hp" "(I (expt hp 2))")
> (equal? (formula-model-fit as-in-r) (formula-model-fit quadratic)) #t
1.6.12.3 R’s hp^2
In a formula ^ crosses terms, and hp crossed with itself is hp, so R’s mpg ~ hp + hp^2 is the linear model mpg ~ hp. ^ is crossing here too:
> (define trap (formula-path (mpg . ~ . hp + hp ^ 2) mtcars)) > (formula-model-predictor-names trap) '("hp")
> (equal? (formula-model-fit trap) (formula-model-fit (formula-path (~ mpg hp) mtcars))) #t
R’s glmnet cannot fit this one, since it needs at least two columns; formula-path fits the one predictor.
1.6.12.4 A logarithm
Any function is a transform. R’s mpg ~ log(hp) + wt is (mpg . ~ . (log hp) + wt):
> (define log-path (formula-path (mpg . ~ . (log hp) + wt) mtcars)) > (coef log-path #:lambda 0.1)
'(("(Intercept)" . 58.77272600056763)
("(log hp)" . -5.798449939727945)
("wt" . -3.2246578459578537))
At the end of the path the logarithm of horsepower and the weight explain 85.9% of the deviance, where horsepower itself and the weight, (mpg . ~ . hp + wt), explain 82.7%:
> (define (final-deviance-ratio model) (define ratios (deviance-ratio model)) (vector-ref ratios (sub1 (vector-length ratios)))) > (final-deviance-ratio log-path) 0.8591063852776801
> (final-deviance-ratio (formula-path (mpg . ~ . hp + wt) mtcars)) 0.826745845280224
1.6.12.5 New data
predict evaluates the transforms again, on the new table’s columns:
> (define new-cars (list (cons "hp" '(100 200 300)) (cons "wt" '(2.5 3.5 4.0)))) > (predict quadratic new-cars #:lambda 0.01) '(23.27797866696576 14.663631758775054 14.131558970704468)
> (predict log-path new-cars #:lambda 0.1) '(24.008232598293418 16.76439552499516 12.800997470344662)
The quadratic flattens out: from 100 to 200 horsepower it predicts 8.6 miles per gallon fewer, and from 200 to 300 only half a mile fewer. At λ = 0.01 its curve bottoms out near 257 horsepower and rises after it, which is the quadratic’s shape more than the cars’: only two of the 32 cars have more horsepower. The quadratic reads only "hp" from the new table, and the logarithm model both columns.
1.6.12.6 Matching R
The same fits in R are
x <- model.matrix(mpg ~ hp + I(hp^2), mtcars)[, -1] |
fit <- glmnet(x, mtcars$mpg) |
coef(fit, s = c(1, 0.1, 0.01)) |
new <- data.frame(hp = c(100, 200, 300), wt = c(2.5, 3.5, 4)) |
predict(fit, newx = model.matrix(~ hp + I(hp^2), new)[, -1], s = 0.01) |
and the same for mpg ~ log(hp) + wt. R 4.5.3 with glmnet 4.1.10 gives the same paths, 78 λ values from 4.604 down to 0.003565 for the quadratic and 55 from 5.147 down to 0.03386 for the logarithm, and the coefficients and predictions above. The example’s companion test checks each against R’s numbers to six significant digits.