Answer

RMSE vs standard deviation: what does comparing them tell you?

Inspired by a question on Cross Validated ·

regressionstandard deviation

The short answer

Yes, the comparison makes sense, but it means the opposite of what is often assumed. The standard deviation of the observed values is (almost exactly) the RMSE you would get by predicting the mean for everyone. So a model whose RMSE is close to that SD has learned next to nothing, and one whose RMSE is larger is worse than the mean. The useful lower limit is the size of the unpredictable noise, not the SD of the outcome. In-sample, 1 − (RMSE/SD)² is the model's R².

The short answer

Root mean squared error (RMSE) is the typical size of a model's prediction errors: square each error (observed minus predicted), average the squares, and take the square root. The standard deviation of the outcome is built the same way, except that every "prediction" is the mean. So the two are directly comparable, in the same units, and the SD is simply the RMSE of the simplest possible model: one that ignores every predictor and guesses the average.

That turns the comparison into a benchmark. If your model's RMSE is about equal to the SD, it predicts no better than the mean. If it is smaller, the model explains some of the variation, and the further below the SD it falls, the more it explains. If it is larger, the model is doing worse than the mean, which happens with badly biased or overfitted predictions.

RMSE ≈ SD signals a useless model, not a perfect one

It is tempting to read "my error is the same size as the natural variation" as "my model has captured everything there is to capture". The logic runs the other way. The SD of the outcome contains two parts: variation that could in principle be predicted from the available information, and variation that is pure noise. A model that explains nothing leaves both parts in its errors, so its RMSE equals the SD.

The real floor is the noise alone. No model can predict below the standard deviation of the part of the outcome that is unrelated to the predictors, so a good model has an RMSE close to that noise SD, which is smaller than the SD of the outcome whenever the predictors matter. In practice the noise SD is unknown, which is why RMSE is judged against the SD (how much better than the mean?) and against competing models, ideally on data the model has not seen.

The link to R-squared

For a least-squares regression with an intercept, evaluated on the data it was fitted to, R² = 1 − SSE/SST, where SSE is the sum of squared residuals and SST the sum of squared deviations from the mean. Dividing both by n gives R² = 1 − (RMSE / SD)², with the SD computed using n in the denominator. So RMSE/SD and R² carry the same information: an RMSE at 70% of the SD means about half the variance is explained, since 1 − 0.7² = 0.51.

Two small details trip people up. Software usually reports the SD with n − 1 in the denominator, so the identity is only approximate with that version. And R's sigma() (the "residual standard error" in summary(lm)) divides the sum of squared residuals by n − p (p = number of coefficients), so it is slightly larger than the plain RMSE. The same idea appears in hydrology as the Nash-Sutcliffe efficiency, 1 − SSE/SST, and as the ratio of RMSE to the observations' SD.

This page compares the error of a model with the spread of the data. For why the standard deviation squares deviations in the first place, which is the reason both measures share a formula, see Why does standard deviation square the differences?.

See it in R and Python

The code simulates 200 observations in which y rises with x plus normal noise with a true SD of 3. It compares the outcome's SD with the RMSE of three sets of predictions: the mean, a fitted regression, and a line whose slope is twice too steep. It then checks the R² identity and repeats the comparison on a held-out half of the data.

R

set.seed(242787)

# Simulate: y depends on x, plus noise with a true SD of 3
n <- 200
x <- runif(n, 0, 10)
y <- 5 + 1.2 * x + rnorm(n, sd = 3)

rmse <- function(obs, pred) sqrt(mean((obs - pred)^2))
sd_n <- sqrt(mean((y - mean(y))^2))      # SD with denominator n

# 1. Three sets of predictions
fit <- lm(y ~ x)
round(c(sd_y        = sd(y),
        rmse_mean   = rmse(y, mean(y)),              # always predict the mean
        rmse_model  = rmse(y, fitted(fit)),          # the fitted regression
        rmse_bad    = rmse(y, 5 + 2.4 * x),          # slope twice too steep
        sigma_fit   = sigma(fit)), 3)                # residual SE, denominator n - 2

# 2. In-sample R-squared is one minus the squared ratio of RMSE to SD
round(c(r2_from_ratio = 1 - (rmse(y, fitted(fit)) / sd_n)^2,
        r2_lm         = summary(fit)$r.squared), 4)

# 3. Out of sample: fit on the first 100 rows, score the other 100
train <- 1:100; test <- 101:200
fit_tr <- lm(y ~ x, subset = train)
pred_te <- predict(fit_tr, newdata = data.frame(x = x[test]))
round(c(sd_test        = sd(y[test]),
        rmse_mean_test = rmse(y[test], mean(y[train])),
        rmse_model_test = rmse(y[test], pred_te)), 3)

Python

import numpy as np

rng = np.random.default_rng(242787)

# Simulate: y depends on x, plus noise with a true SD of 3
n = 200
x = rng.uniform(0, 10, n)
y = 5 + 1.2 * x + rng.normal(0, 3, n)

def rmse(obs, pred):
    return np.sqrt(np.mean((obs - pred) ** 2))

def fit_line(xs, ys):                     # least-squares intercept and slope
    X = np.column_stack([np.ones_like(xs), xs])
    return np.linalg.lstsq(X, ys, rcond=None)[0]

sd_n = np.sqrt(np.mean((y - y.mean()) ** 2))   # SD with denominator n

# 1. Three sets of predictions
b = fit_line(x, y)
fitted = b[0] + b[1] * x
print("sd_y", round(np.std(y, ddof=1), 3),
      "rmse_mean", round(rmse(y, y.mean()), 3),
      "rmse_model", round(rmse(y, fitted), 3),
      "rmse_bad", round(rmse(y, 5 + 2.4 * x), 3),
      "sigma_fit", round(np.sqrt(np.sum((y - fitted) ** 2) / (n - 2)), 3))

# 2. In-sample R-squared is one minus the squared ratio of RMSE to SD
r2_lm = 1 - np.sum((y - fitted) ** 2) / np.sum((y - y.mean()) ** 2)
print("r2_from_ratio", round(1 - (rmse(y, fitted) / sd_n) ** 2, 4),
      "r2_lm", round(r2_lm, 4))

# 3. Out of sample: fit on the first 100 rows, score the other 100
tr, te = slice(0, 100), slice(100, 200)
b_tr = fit_line(x[tr], y[tr])
pred_te = b_tr[0] + b_tr[1] * x[te]
print("sd_test", round(np.std(y[te], ddof=1), 3),
      "rmse_mean_test", round(rmse(y[te], y[tr].mean()), 3),
      "rmse_model_test", round(rmse(y[te], pred_te), 3))

The figures below come from one seeded run of the R code. Python's random numbers differ from R's, so its figures differ slightly (for example, an SD of 4.68 and a model RMSE of 3.15), but they show the same pattern, and its two R² values also match each other exactly.

For more on why held-out data matters, see Training, validation and test sets explained, and for the limits of R² itself, Is R-squared useful or misleading?. More plain-language guides are on the DASS blog.

How to report RMSE in APA style (7th edition)

Give RMSE in the outcome's units, say whether it was computed on the fitting data or on held-out data, and set it beside a benchmark such as the outcome's SD or the RMSE of predicting the mean. With the example above:

"On the held-out test set (n = 100), the regression model's predictions had an RMSE of 2.86, compared with 4.29 for a baseline that predicted the training-set mean (outcome SD = 4.23)."

For in-sample fit you can add "R² = .58". RMSE can exceed 1, so it keeps its leading zero when below 1 (e.g., 0.42), whereas R² cannot exceed 1 and takes none. Define RMSE at first use (root mean square error) and, in the method section, state how the data were split or cross-validated.

Related tools and guides

More answered questions

Working with your own data?

General answers only go so far. Send us your situation and we'll reply by email within two business days.

Ask your question

Written by AskStats with AI assistance. This is general information, not advice for your specific data or study. When the results matter, check your approach with a qualified statistician.