2
$\begingroup$

I have a burning question. First, in Python:

import os
import time
import numpy as np
from sklearn.linear_model import LassoCV
from sklearn.datasets import make_regression


X, y = make_regression(1000, 5000, noise = 4, random_state = 123)
if not os.path.exists("tmp"): os.mkdir("tmp")


# Save on disk to be read in R session.
np.save("tmp/X.npy", X)
np.save("tmp/y.npy", y)


# Let first 70% rows be training examples.
Nrow = int(round(len(X) * 0.7))


tik = time.time()
reg = LassoCV(cv = 5).fit(X[:Nrow, :], y[:Nrow])
time.time() - tik 
# 6.59s


# Compute mse on the validation examples:
np.mean((reg.predict(X[Nrow:, :]) - y[Nrow:]) ** 2)
# 19.036

Then in R:

X = RcppCNPy::npyLoad("tmp/X.npy") # Load the data saved from Python session.
y = RcppCNPy::npyLoad("tmp/y.npy")


# First 70% rows are training data.
trainInd = 1:as.integer(round(nrow(X) * 0.7))


system.time({
reg = glmnet::cv.glmnet(
  X[trainInd, ], y[trainInd], type.measure = "mse", nfolds = 5)
})
# 1.047s


mean((predict(reg, X[-trainInd, ], s = "lambda.min") - y[-trainInd]) ^ 2)
# 32.01

Why does R glmnet result in an MSE 68% higher than that from Python sklearn? I anticipated glmnet would be faster, but a 6.6x speedup over sklearn seems to suggest glmnet skipped something important, which might contribute to its significantly worse error metric? I couldn't figure out where I did wrong..

$\endgroup$
4
  • 1
    $\begingroup$ To understand what's going on, compute the variance of Y and then compare var(Y) to the MSE. And it always helps to look at plots. In your case plotting the predictions against the true Ys will also help you figure out things. $\endgroup$
    – dipetkov
    Commented Jan 25, 2023 at 8:40
  • $\begingroup$ @dipetkov Thanks! Please post ur answer if u'd like and i'll accept it $\endgroup$ Commented Jan 25, 2023 at 15:20
  • $\begingroup$ So what did actually happen here? Where did the discrepancy come from? $\endgroup$ Commented Jan 25, 2023 at 17:49
  • 1
    $\begingroup$ @RichardHardy 32.01 and 19.036 are trivial compared to y's magnitude.. $\endgroup$ Commented Jan 25, 2023 at 17:56

1 Answer 1

3
$\begingroup$

As always, it helps to make plots. Here I plot the cv.glmnet predictions on the left and the LassoCV predictions on the right, against the true values in the hold out set.

We see immediately that both implementations do a great job: they explain most of the variability in $y$ and the mean squared error (MSE) is small compared to the total variability in $y$.

We can quantify the difference between the implementations in terms of the proportion of variance explained (PVE) = 1 - MSE / Var$(y_{\text{holdout}})$.

  • PVEcv.glmnet = .9989
  • PVELassoCV = .9993

Note: PVE is not the same as the R-squared, which is computed on the training data. And it's not equivalent to comparing the Lasso to the "naive" model that predicts the average $y$ observed in the training data.

So there is a difference but it is small relative to the variability of the outcome variable. I wouldn't use this example to conclude that LassoCV is a better implementation than cv.glmnet in general.

PS. An increase of .0004 could be a meaningful improvement and worthwhile to put in practice. However, to conclude that a new model is better than the production model by such a small increment, we would want to compare the models on more than 300 examples.

$\endgroup$
4
  • $\begingroup$ Do you think the discrepancy is because of some kind of difference in the numerical optimization of the LASSO parameters? $\endgroup$
    – Dave
    Commented Jan 25, 2023 at 20:53
  • 2
    $\begingroup$ @Dave Or the default settings perhaps? One could repeatedly simulate data with a different seed and see how often LassoCV does better. If it was easier to jump between R and python.... But I'm not sure what we'll learn from the exercise. Visually, on this dataset, LassoCV does better at predicting the extremes. $\endgroup$
    – dipetkov
    Commented Jan 25, 2023 at 20:57
  • 1
    $\begingroup$ There’s always the R package reticulate! $\endgroup$
    – Dave
    Commented Jan 25, 2023 at 20:59
  • 2
    $\begingroup$ @Dave as well as the python package rpy2 :) $\endgroup$ Commented Jan 25, 2023 at 22:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.