Tutorial ridge regression

This tutorial gives a brief introduction to ridge regression. The tutorial makes use of the acetylene dataset from "Marquardt, D. W., and Snee, R. D. (1975). “Ridge Regression in Practice.” American Statistician 29:3–20.", and the follow the same outline as the SAS documentation

First, creating the dataset.

We create the dataset with the help of the DataFrames.jl package.

using DataFrames
x1 = [1300, 1300, 1300, 1300, 1300, 1300, 1200, 1200, 1200, 1200, 1200, 1200, 1100, 1100, 1100, 1100]
x2 = [7.5, 9.0, 11.0, 13.5, 17.0, 23.0, 5.3, 7.5, 11.0, 13.5, 17.0, 23.0, 5.3, 7.5, 11.0, 17.0]
x3 = [0.012, 0.012, 0.0115, 0.013, 0.0135, 0.012, 0.04, 0.038, 0.032, 0.026, 0.034, 0.041, 0.084, 0.098, 0.092, 0.086]
y = [49.0, 50.2, 50.5, 48.5, 47.5, 44.5, 28.0, 31.5, 34.5, 35.0, 38.0, 38.5, 15.0, 17.0, 20.5, 29.5]
df = DataFrame(x1= x1, x2= x2, x3= x3, y= y)

Second, make a least square regression

We make a ordinary least squares (OLS) regression for comparison.

using LinearRegressionKit, StatsModels
using VegaLite

f = @formula(y ~ x1 + x2 + x3 + x1 & x2 + x1^2)
lm, ps = regress(f, df, "all", req_stats=["default", "vif"])
lm

Model definition:	y ~ 1 + x1 + x2 + x3 + :(x1 ^ 2) + x1 & x2
Used observations:	16
Model statistics:
  R²: 0.993708			Adjusted R²: 0.990562
  MSE: 1.33625			RMSE: 1.15596
  σ̂²: 1.33625
  F Value: 315.86 with degrees of freedom 5 and 10, Pr > F (p-value): 1.14778e-10
Confidence interval: 95%

Coefficients statistics:
Terms ╲ Stats │       Coefs      Std err            t     Pr(>|t|)         code       low ci      high ci          VIF
──────────────┼───────────────────────────────────────────────────────────────────────────────────────────────────────
(Intercept)   │     390.538      211.523      1.84632    0.0946146           .      -80.7641      861.841          0.0
x1            │   -0.776761     0.324481     -2.39385    0.0377095           *      -1.49975   -0.0537716      7682.37
x2            │     10.1735     0.943009      10.7883   7.89586e-7         ***       8.07235      12.2747      320.022
x3            │    -121.626      69.0175     -1.76225     0.108507                  -275.407      32.1545      53.5246
x1 ^ 2        │ 0.000398306  0.000125283      3.17925   0.00983192          **   0.000119158  0.000677454      6643.32
x1 & x2       │  -0.0080548  0.000772088     -10.4325   1.07686e-6         ***   -0.00977512  -0.00633448      344.545

	Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

We observe that the VIF for the coefficients are highs, and hence indicate likely multicollinearity.

Ridge regression

Ridge regression requires a parameter (k), while there are methods to numerically identify k. It is also possible to trace the coefficients and VIFs values to let the analyst choose a k. Here we are going to trace for the k between 0.0 and 0.1 by increment of 0.0005. We display only the results for the first 5 k.

rdf, ps = ridge(f, df, 0.0:0.0005:0.1, traceplots=true)
rdf[1:5 , :]

Here is the default trace plot for the coefficients:

ps["coefs traceplot"]

And here is the default trace plot for the VIFs:

ps["vifs traceplot"]

As it is difficult to see the VIFs traces, it is also possible to request the version of hte plot with the y-axis log scaled.

ps["vifs traceplot log"]

Once a k has been selected (in this case 0.004) a regular ridge regression with this value can executed.

rlm = ridge(f, df, 0.004)

Ridge Regression
Constant k:	0.004
Model definition:	y ~ 1 + x1 + x2 + x3 + :(x1 ^ 2) + x1 & x2
Used observations:	16
Model statistics:
  R²: 0.951075			Adjusted R²: 0.926612
  MSE: 10.3903			RMSE: 3.2234

Coefficients statistics:
Terms ╲ Stats │       Coefs          VIF
──────────────┼─────────────────────────
(Intercept)   │     -93.797          0.0
x1            │   0.0578795      4.35597
x2            │     2.83932      23.9391
x3            │    -21.3375      9.99599
x1 ^ 2        │  3.82194e-5       5.1494
x1 & x2       │ -0.00205277      25.7445

From there the regular predict_* functions can be used, although only the predicted and potentially the residuals statistics will be calculated.

res = predict_in_sample(rlm, df)