Почему кросс-валидация в scikit возвращает отрицательное R^2 для моих сильно коррелированных данных?

Question 1

У меня есть следующий предварительно обработанный набор данных в небольшом датафрейме Pandas:

df = pd.DataFrame({'protesters_estimate': [-0.025277219233618427,  -1.1209981799804989, -0.7584382556131943, 0.3357505651218889, -0.4730177055285195, -1.0359399607957092, 0.4222519531704431, 0.3686864558273824, -0.3912614033608771, 3.4038543116306723, -0.7557908523751271, 0.8585956431868133, -0.36278846361162215, -0.7752270840612929, 1.434618254960961, -0.5714029624055086, 0.42736811052416984, -0.13314848352841915, -0.7043146521442091, 0.9568348193244159, -0.37506384503413215, -0.7252910460740197],

'climate_mentions': [0.014620553379808637, 0.014523096436234, 0.009729460025372921, 0.019436345966958212, 0.017640028535268398, 0.011277283001333295, 0.019327228087934106, 0.015160273972011584, 0.021741381095351484, 0.03670954832162733, 0.005625206385871155, 0.02724524233334451, 0.008399949825246915, 0.012832028743744384, 0.015122959766868112, 0.009519105767012719, 0.012243145968000543, 0.02004639571430332, 0.016822301175172555, 0.012241593155590624, 0.004016740988627743, 0.01021722632961321]})

Когда я запускаю df.corr(), я получаю корреляцию 0.7 между двумя переменными. Но когда я запускаю следующий код (с использованием sklearn):

X = df[['protesters_estimate']]
y = df['climate_mentions']
print(f"Cross val score: {cross_val_score(LinearRegression(), X, y).mean()}")

Я получаю вывод ‘Cross val score: -0.4412345093041985’.

Что здесь происходит? Как я понимаю, R^2 должен быть буквально (0.7)^2 для такой линейной регрессии, и если шум вносится разбиением cross-val, то его можно ожидать в пределах +- десятичного разряда или двух, а не отрицательной оценки.

Далее запутывая меня, если я фактически обучу модель и проверю её:

model = LinearRegression()
model.fit(X, y)
print(f"R-squared: {model.score(X, y)}")

Я получаю вывод, который я ожидал: ‘R-squared: 0.48415257262102906’

Question 2

Я закодировал эксперимент на R, который использует ваши данные и вручную выполняет некоторое перекрестное тестирование.

(Ладно, я не совсем выполнял перекрестное тестирование, но это повторное использование различных наборов для удержания, чтобы увидеть, как себя ведет вневыборочная $R^2$.)

library(ggplot2)
set.seed(2024)
x <- c(-0.025277219233618427,  -1.1209981799804989, -0.7584382556131943, 0.3357505651218889, -0.4730177055285195, -1.0359399607957092, 0.4222519531704431, 0.3686864558273824, -0.3912614033608771, 3.4038543116306723, -0.7557908523751271, 0.8585956431868133, -0.36278846361162215, -0.7752270840612929, 1.434618254960961, -0.5714029624055086, 0.42736811052416984, -0.13314848352841915, -0.7043146521442091, 0.9568348193244159, -0.37506384503413215, -0.7252910460740197)
y <- c(0.014620553379808637, 0.014523096436234, 0.009729460025372921, 0.019436345966958212, 0.017640028535268398, 0.011277283001333295, 0.019327228087934106, 0.015160273972011584, 0.021741381095351484, 0.03670954832162733, 0.005625206385871155, 0.02724524233334451, 0.008399949825246915, 0.012832028743744384, 0.015122959766868112, 0.009519105767012719, 0.012243145968000543, 0.02004639571430332, 0.016822301175172555, 0.012241593155590624, 0.004016740988627743, 0.01021722632961321)
n_train <- 17
ncv <- 75
r2_sklearn <- function(y, yhat){
  num <- sum((y - yhat)^2)
  den <- sum((y - mean(y))^2)
  return(1 - (num/den))
}
r2_dave <- function(y, yhat, ybar){
  num <- sum((y - yhat)^2)
  den <- sum((y - ybar)^2)
  return(1 - (num/den))
}
r2s_sklearn <- r2s_dave <- rep(NA, ncv)
for (i in 1:ncv){
  
  idx <- sample(
    seq(1, length(x), 1),
    n_train,
    replace = F
  )
  
  x_train <- x[idx]
  y_train <- y[idx]
  x_test <- y[-idx]
  y_test <- y[-idx]
  
  L <- lm(y_train ~ x_train)
  yhat <- (coef(L) %*% t(cbind(1, x_test)))[1, ]
  
  r2s_sklearn[i] <- r2_sklearn(y_test, yhat)
  r2s_dave[i]    <- r2_dave(   y_test, yhat, mean(y_train))

  
}
d_sklearn <- data.frame(
  R2 = r2s_sklearn,
  CDF = ecdf(r2s_sklearn)(r2s_sklearn),
  Calculation = "Scikit-learn"
)
d_dave <- data.frame(
  R2 = r2s_dave,
  CDF = ecdf(r2s_dave)(r2s_dave),
  Calculation = "Dave"
)
d <- rbind(d_dave, d_sklearn)
ggplot(d, aes(x = R2, y = CDF, col = Calculation)) +
  geom_line() +
  theme(legend.position = "bottom")

Когда вы смотрите на функции распределения для $75$ итераций тренировки и затем тестирования на новых данных, действительно, расчет $R^2$ из sklearn обычно отрицателен. Это соответствует тому, что вы наблюдаете в вашем коде на Python.

Причина, по которой вы наблюдаете отрицательные значения $R^2$, заключается в том, что sklearn не вычисляет $R^2$ путем возведения в квадрат. Как я обсуждаю здесь, sklearn вычисляет следующее.

$$
R^2_{\text{out-of-sample, scikit-learn}}=
1-\left(\dfrac{
\overset{N}{\underset{i=1}{\sum}}\left(
y_i-\hat y_i
\right)^2
}{
\overset{N}{\underset{i=1}{\sum}}\left(
y_i-\bar y_{\text{input-}y}
\right)^2
}\right)
$$

Если числитель дроби превышает знаменатель, весь расчет отрицателен, показывая, что модель имеет худшую производительность (с точки зрения мер потерь, таких как MSE и RMSE), чем простая эталонная модель, которая всегда предсказывает среднее от входных данных.

Такой расчет делает API простым, но, на мой взгляд, лишает его полезной интерпретации для данных вне выборки, как я обсуждаю в той же ссылке, и я бы выбрал расчет $R^2$ для данных вне выборки следующим образом.

$$
R^2_{\text{out-of-sample, Dave}}=
1-\left(\dfrac{
\overset{N}{\underset{i=1}{\sum}}\left(
y_i-\hat y_i
\right)^2
}{
\overset{N}{\underset{i=1}{\sum}}\left(
y_i-\bar y_{\text{in-sample}}
\right)^2
}\right)
$$

Это дает красную линию “Dave” на графике, которая в основном положительна. Действительно, можно строго показать, что расчет sklearn является нижней границей моего расчета. (Это упражнение полезно.)

Возможно объяснение, почему вы получаете какие-либо значения $R^2$ ниже нуля, заключается в том, что происходит, когда вы исключаете определенную точку, находящуюся в верхнем правом углу.

d_data <- data.frame(x, y)
ggplot(d_data, aes(x = x, y = y)) +
  geom_point()

На всем наборе данных корреляция довольно высокая, почти $0.7$ для регрессии $R^2$ около $0.49$. Однако, исключение этой одной точки с использованием cor(x[y < max(y)], y[y < max(y)]) снижает корреляцию до гораздо более низкой величины около $0.41$ для регрессионной оценки $R^2$ равной $0.16$. Тренд становится намного слабее, когда эта одна точка отсутствует, что будет происходить четыре из пяти раз в вашем перекрестном тестировании. Если другие точки в тренировочном наборе данных приводят к довольно плоской или отрицательной наклонной линии регрессии, то эта верхняя точка будет предсказана очень плохо, с соответствующей высокой потерей квадратов и низким, возможно даже отрицательным, значением $R^2$, особенно когда вы используете расчет sklearn, который является нижней границей предпочитаемой мной эксплуатации.

(Опять же, эта ссылка объясняет, почему мне больше нравится мой расчет, включая ссылку на литературу, которая выступает в поддержку моей точки зрения. (Эта ссылка не на мою статью, хотя я периодически (в основном не) пишу статью о едином представлении $R^2$, в которой буду касаться вопросов, связанных с работой вне выборки.))

Почему кросс-валидация в scikit возвращает отрицательное R^2 для моих сильно коррелированных данных?

Вопрос или проблема

Ответ или решение