During a regression task, I am getting low $R^2$ values, but elementwise difference between test set and prediction values is huge

Question

I am doing a random forest regression on my dataset (which has abut 15 input features and 1 target feature). I am getting a decently low $R^2$ of <1 for both the train and test sets (please do let me know if <1 is not a good-enough $R^2$ score).

import pandas as pd
import numpy as np
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split

# load dataset
df = pd.read_csv('Dataset.csv')

# split into input (X) and output (Y) variables
X = df.drop(['ID_COLUMN', 'TARGET_COLUMN'], axis=1)
Y = df.TARGET_COLUMN

# Split the data into 67% for training and 33% for testing
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.33)

# Fitting the regression model to the dataset
regressor = RandomForestRegressor(n_estimators = 100, random_state = 50)
regressor.fit(X_train, Y_train.ravel()) # Using ravel() to avoid getting 'DataConversionWarning' warning message


print("Predicting Values:")
y_pred = regressor.predict(X_test)

print("Getting Model Performance...")

# Get regression scores
print("R^2 train = ", regressor.score(X_train, Y_train))
print("R^2 test = ", regressor.score(X_test, Y_test))

This outputs the following:

Predicting Values:
Getting Model Performance...
R^2 train =  0.9791000275450427
R^2 test = 0.8577464692386905

Then, I checked the difference between the actual target column values in the test dataset versus the predicted values, like so:

diff = []
for i in range(len(y_pred)):
    if Y_test.values[i]!=0: # a few values were 0 which was causing the corresponding diff value to become inf
        diff.append(100*np.abs(y_pred[i]-Y_test.values[i])/Y_test.values[i]) # element-wise percentage error

I found that the majority of the element-wise differences were between 40-60% and their mean was almost 50%!

np.mean(diff)
>>> 49.07580695857447

So, which one is correct? Is the regression score correct and my model is good for this data, or is the element-wise error I calculated correct and the model didn't do well for this data? If its the latter, please advise on how to increase the prediction accuracy.

---

I also checked the rmse score:

import math
rmse = math.sqrt(np.mean((np.array(Y_test) - y_pred)**2))
rmse
>>> 3.67328471827293

This seems quite high for the model to have done a good job, but please correct me if I'm wrong.

And I also checked the $R^2$ scores for different number of estimators:

import matplotlib.pyplot as plt
model = RandomForestRegressor(n_jobs=-1)
# Try different numbers of n_estimators
estimators = np.arange(10, 200, 10)
scores = []
for n in estimators:
    model.set_params(n_estimators=n)
    model.fit(X_train, Y_train)
    scores.append(model.score(X_test, Y_test))
plt.title("Effect of n_estimators")
plt.xlabel("n_estimator")
plt.ylabel("score")
plt.plot(estimators, scores)

Please advise.

---

I tried using linear regression first, and got a very high MSE (which is why I was trying out random forest):

from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score

lr = LinearRegression()
lr.fit(X_train, y_train)
y_pred = lr.predict(X_test)

# The coefficients
print('Coefficients: \n', lr.coef_)
# The mean squared error
print("Mean squared error: %.2f" % mean_squared_error(y_test, y_pred))
# Explained variance score: 1 is perfect prediction
print('Variance score: %.2f' % r2_score(y_test, y_pred))


Coefficients: 
 [ 1.93829229e-01 -4.68738825e-01  2.01635420e-01  6.35902010e-01
  6.57354434e-03  5.13180293e-03  2.84015810e-01 -1.31469084e-06
  1.95335035e+00]
Mean squared error: 86.92
Variance score: 0.08

Answer 1

For starters, does the split_train_test shuffle the data first? I'd also say that 100 trees is not enough as you can see in your plot, it keeps getting better as the number of trees is growing.

It also seems that you first started by using a linear regression without regularisation (try Ridge instead, or Lasso or elasticNet).

In general, I would say that error measures such as the MSE, RMSE for a model are only possible to interpret once you compare them to the same scores for other models (unless it is either incredibly small or large, with respect to the scale/magnitude of your targets). So, I woul advise you compute MSE using Ridge, and then compare it to the RMSE of RFs.

And to answer your question, if you only fit one model, I would trust the R^2 score more than the MSE for the reason stated above. The R^2 is more easily interpretable on its own and tells you whether your algorithm "solved"/did well for the given problem. Your R^2 scores are pretty good, but you may be overfitting a bit.

Hope this helps

Answer 2

Remember the definition of $R^2$ as it is defined by that sklearn function.

$$ R^2=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 \right)^2 }\right) $$

The numerator of that fraction is proportional to the $MSE$, and the denominator is proportional to the variance of the $y$ observations.

Your objection is to getting an $R^2$ that looks good yet also getting an $MSE$ that looks big.

If that denominator term is extremely gigantic compared to the $MSE$, then it is entirely possible to have an $R^2$ that looks like a strong letter grade in school yet corresponds to a high $MSE$.

Remember that, under certain conditions, $R^2$ describes the proportion of variance explained. If you start out with a gargantuan variance, then you can leave a huge amount of variance unexplained while explaining a large proportion of the variance.

Let's do a simulation. In R:

set.seed(2021)
N <- 1000
a_vals <- seq(1, 10, 0.01)   
err <- rnorm(N)
r2s <- mses <- vars <- rep(NA, length(a_vals))

for (i in 1:length(a_vals)){
        
    x <- runif(N, -a_vals[i], a_vals[i])
    y <- x + err 
    L <- lm(y ~ x)
    r2s[i] <- summary(L)$r.squared
    mses[i] <- mean((predict(L) - y)^2)    
    vars[i] <- var(y)
}
plot(vars, r2s)
plot(vars, mses)
plot(ecdf(vars)(vars), ecdf(mses)(mses)) # Copula is more-or-less flat, so independence

If you plot these (I can't save and post pictures right now), you will see that, as you increase the variance of $y$, you increase the $R^2$. However, variance of $y$ is independent of the $MSE$.

(There's this other issue where you seem to be looking at mean absolute error, rather than mean square error. That is more complicated, because there is not a perfect relationship between $MAE$ and $R^2$. However, if you're noticing that the deviations between truth and prediction are large, in some sense, that corresponds to a large $MSE$.)