3
$\begingroup$

I am performing a simple linear regression and have started to examine the performance of my model. One action I've taken is to stratify the dependent variable into deciles and summarize the performance of the model across those deciles:

library(dplyr)
library(ggplot2)

set.seed(100)

r2 <- function(y_true, y_pred){
    return(cor(y_true, y_pred) ** 2)
}

rmse <- function(y_true, y_pred){
    return(sqrt(mean((y_true - y_pred) ** 2)))
}

data(cars)

model <- lm(speed ~ dist, data = cars)

y_true <- cars$speed
y_pred <- predict(model, data = cars)

r2(y_true, y_pred)
# 0.6510794

results_df <- data.frame(
    this_y_true = y_true,
    this_y_pred = y_pred,
    decile = as.factor(ntile(y_true, 10))
)

gb <- as.data.frame(
    results_df %>%
        group_by(decile) %>%
        mutate(
            mean_y_true = mean(this_y_true),
            sd_y_pred = sd(this_y_true),
            mean_y_pred = mean(this_y_pred),
            sd_y_pred = sd(this_y_pred),
            n = sum(!is.na(this_y_true)),
            r2 = r2(this_y_true, this_y_pred),
            rmse = rmse(this_y_true, this_y_pred)
        ) %>%
        select(-one_of(c('this_y_true', 'this_y_pred'))) %>%
        arrange(decile) %>%
        unique()
)
gb
#    decile mean_y_true sd_y_pred mean_y_pred n          r2     rmse
# 1       1         6.0 1.3772993    10.07204 5 0.298307184 4.318254
# 2       2        10.0 1.5264576    11.76082 5 0.072058824 2.213814
# 3       3        11.8 0.9823056    12.05885 5 0.240056818 1.159017
# 4       4        13.2 1.3572502    13.78075 5 0.241071429 1.564547
# 5       5        14.4 4.1517534    15.63511 5 0.606923240 4.288315
# 6       6        16.2 1.4882674    14.84038 5 0.515558699 2.368759
# 7       7        17.8 2.9375999    18.48287 5 0.133576874 2.600321
# 8       8        19.4 2.3179460    15.90001 5 0.153061224 4.193200
# 9       9        21.0 1.0313146    17.95305 5 0.003221649 3.406219
# 10     10        24.2 3.0054049    23.51612 5 0.046471927 2.883963

I have summarized the R-squared and RMSE values across each decile. One thing I noticed was the substantial variation in R-squared values across deciles – while the overall regression has an R-squared of 0.65, no decile yields an R-squared higher than that. I don't understand how none of the individual deciles have at least as high an R-squared as the overall regression. How can this be? Is my model useless?

I don't think the model is useless. Here is my line of thinking:

ggplot(data = results_df, aes(x = this_y_pred, y = this_y_true)) +
    geom_point() +
    geom_smooth(method = lm, se = FALSE) +
    ylab('Actual') +
    xlab('Predicted')

enter image description here

The above graph is a simple scatterplot of predicted vs. actual y values, with a line of least squares added. There is definitely a positive linear association between the model and the outcome (this is where the 0.65 R-squared comes from).

Now let's look at the same plot, but this time we'll color the points and add a best-fit line by decile:

ggplot(data = results_df, aes(x = this_y_pred, y = this_y_true, color = decile)) +
    geom_point() +
    geom_smooth(method = lm, se = FALSE) +
    ylab('Actual') +
    xlab('Predicted')

enter image description here

Whoa – at a by-decile level, any positive association between the actual and predicted values seems to have disappeared, as evidence by the nearly horizontal lines of best fit. This is where those low intra-decile R-squareds are coming from.

I notice, though, that some of the deciles with the highest R-squareds also have the highest RMSEs. Looking at the by-decile graph, we see that the top two deciles by R-squared (5 and 6) look dramatically different visually. We also note decile 5's high RMSE vs. 6's relatively low RMSE. The horizontal dispersion of the two deciles' points gives us an idea as to why 6 has a lower RMSE.

I'm having trouble taking the last step, though, and reconciling the high R-squareds of deciles 5 and 6, for example, and their dramatically different RMSEs. I'm beginning to think R-square should not be used in this way – looked at by-decile – but perhaps there's more nuance.

I know this is a small sample size (only 5 points per decile), but this mirrors a real-world problem I'm dealing with.

Additional resources I've found on this topic include:

$\endgroup$
2
+50
$\begingroup$

I like the way you have investigated this. You'll want also to look into the concept of restriction of range. RSQ within any one of such deciles is bound to be attenuated due to restriction of range. You need full variability in X to see full correlation with Y.

What you've found doesn't invalidate your overall model. And if you simulate further examples that mimic these conditions you are going to find the same sort of results. With sample sizes per decile much greater than 5, you'd find within-decile RSQ to be more consistent, but still low. If you divided your data into fewer and fewer groups instead of 10, you'd see within-group RSQs climb to more and more closely match the 0.65.

$\endgroup$
1
  • $\begingroup$ Thank you for formalizing this idea – I had in my mind that it had something to do with the decile groupings automatically limiting how much variation was possible, but didn't know the correct term, "restriction of range." $\endgroup$ – blacksite Apr 27 '20 at 14:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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