R ile Çoklu Doğrusal Regresyon: Modelden Tahmine

Fitting a regression model in R takes one line. The real work starts after that line: what do the coefficients say, do the assumptions hold, and how well does the model perform on data it has never seen?

In this second post of the series I fit a multiple linear regression on the Carseats data we prepared in the first post — from the model to coefficient interpretation, and from assumption checks to real performance on the test set.

“A high R² does not mean a good model; it only means the model fits the training data well.”

1. Fitting the Model

The general form of multiple linear regression is:

In the Carseats data the target is Sales, unit sales per store. We split 80% training / 20% test and fit the model on the training set only.

library(ISLR)
data(Carseats)

set.seed(5364)
idx <- sample(seq_len(nrow(Carseats)), size = floor(0.8 * nrow(Carseats)))
train <- Carseats[idx, ]
test  <- Carseats[-idx, ]

lm_fit <- lm(Sales ~ CompPrice + Income + Advertising + Price + ShelveLoc + Age,
             data = train)

summary(lm_fit)

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)
(Intercept)      5.461708   0.581852   9.387  < 2e-16 ***
CompPrice        0.093870   0.004771  19.674  < 2e-16 ***
Income           0.016768   0.002122   7.902 4.76e-14 ***
Advertising      0.124613   0.008737  14.263  < 2e-16 ***
Price           -0.097200   0.003116 -31.192  < 2e-16 ***
ShelveLocGood    4.813149   0.172728  27.866  < 2e-16 ***
ShelveLocMedium  1.978825   0.145542  13.596  < 2e-16 ***
Age             -0.047426   0.003690 -12.852  < 2e-16 ***

Residual standard error: 1.044 on 312 degrees of freedom
Multiple R-squared:  0.8719,	Adjusted R-squared:  0.869
F-statistic: 303.4 on 7 and 312 DF,  p-value: < 2.2e-16

The model as a whole is significant: F = 303.4 with a p-value below 2.2e−16. All seven variables are individually significant too. Adjusted R² is 0.869, so roughly 87% of the variation in sales in the training set is explained by these seven variables.

2. Interpreting the Coefficients

Written out, the model is:

Each coefficient shows the effect of a one-unit increase in that variable on sales, holding the other variables constant.

The striking finding here is shelving location: according to this model, putting the product on a good shelf makes a far bigger difference than raising the advertising budget. Matching a 4.81 thousand-unit effect through advertising alone would take roughly 39 thousand dollars of extra spend.

Which variable matters most?

Reading the table and saying "shelf location has the largest effect" is fair. Comparing Advertising (0.125) with Income (0.017) is not. Each coefficient is measured in its own unit: one is a thousand dollars of advertising budget, another a thousand dollars of income, a third is years. When the units differ, the size of the numbers cannot be compared.

The fix is to put every numeric variable on a common scale and refit. A standardised coefficient answers a different question: when this variable rises by one standard deviation, how many standard deviations does sales move?

scaled <- train
numeric_vars <- c("Sales", "CompPrice", "Income",
                   "Advertising", "Price", "Age")
scaled[numeric_vars] <- lapply(scaled[numeric_vars], scale)

round(coef(lm(Sales ~ CompPrice + Income + Advertising +
                 Price + ShelveLoc + Age, data = scaled))[-1], 3)

The ordering is not the one the raw table suggested. There CompPrice (0.094) and Price (−0.097) looked almost equal; on a common scale Price turns out to be roughly 1.6 times stronger. The reason is simple: the competitor's price varies over a narrower range, so in practice it moves sales less.

Advertising lands in its proper place too. Its raw coefficient is seven times that of Income, but its standardised effect is only 1.8 times larger — because income varies over a far wider range. "Large coefficient" and "large effect" are not the same claim.

3. Checking the Assumptions

Linear regression has four core assumptions, and none of them is about whether the model "works" — they are all about whether you can trust the coefficients and p-values.

par(mfrow = c(2, 2))
plot(lm_fit)
par(mfrow = c(1, 1))

# Normality of residuals
shapiro.test(residuals(lm_fit))

# Autocorrelation
library(lmtest)
dwtest(lm_fit)

# Heteroscedasticity
bptest(lm_fit)

plot(lm_fit) produces four plots, each addressing one assumption:

From plots to tests: checking assumptions with numbers

The four plots are a good first look, but they are a judgement call. Whether a band is "flat enough" depends on who is looking. Each assumption also has a formal test, and writing a number in the report is easier to defend than "nothing looked wrong".

e <- residuals(lm_fit)

# 1) Normality of errors
shapiro.test(e)                 # W = 0.9963, p = 0.6546

# 2) Constant variance - Breusch-Pagan (LM = n * R2)
aux <- lm(I(e^2) ~ CompPrice + Income + Advertising +
            Price + ShelveLoc + Age, data = train)
LM <- length(e) * summary(aux)$r.squared          # 3.3946
pchisq(LM, df = 7, lower.tail = FALSE)          # 0.8463

# 3) Independence - Durbin-Watson
sum(diff(e)^2) / sum(e^2)                     # 2.0266

# 4) Functional form - Ramsey RESET
yhat <- fitted(lm_fit)
wider <- lm(Sales ~ CompPrice + Income + Advertising + Price +
              ShelveLoc + Age + I(yhat^2) + I(yhat^3), data = train)
anova(lm_fit, wider)            # F = 0.060, p = 0.9416

All four are clean. The residuals have a skewness of 0.124 and a kurtosis of 3.013 — essentially the 0 and 3 of a normal distribution. Nothing stands in the way of trusting the coefficients and their p-values here.

And if one of the assumptions had failed? Each has a standard answer:

Multicollinearity

When predictors are strongly related to one another the coefficients become unstable: signs can flip and standard errors inflate. The measure is VIF:

Here R²j is the R² from regressing the j-th predictor on all the others. The common thresholds are 5 or 10.

library(car)
vif(lm_fit)

# Result:
# CompPrice   1.546
# Income      1.022
# Advertising 1.018
# Price       1.568
# Age         1.033

Every value sits between 1.0 and 1.6, far below the thresholds. There is no multicollinearity problem in this data set — the coefficients can be interpreted with confidence.

4. Predicting on the Test Set

The model looks good on training data. The real question: what does it do on 80 stores it has never seen?

pred <- predict(lm_fit, newdata = test)

rmse <- sqrt(mean((test$Sales - pred)^2))
mae  <- mean(abs(test$Sales - pred))
r2   <- 1 - sum((test$Sales - pred)^2) /
             sum((test$Sales - mean(test$Sales))^2)

cat(sprintf("RMSE=%.4f MAE=%.4f R2=%.4f\n", rmse, mae, r2))
# RMSE=0.9340 MAE=0.7694 R2=0.8639

A test error below the training error looks odd at first but is perfectly normal: in a test set of 80 observations, "easy" stores can land together by chance. What matters is that the two numbers are close — a large gap would point to overfitting.

How uncertain is a single store's prediction?

RMSE gives the average error across the whole test set. But the question on the table is usually narrower: "what do you predict for this store?" A single number answers it incompletely; the uncertainty around that number belongs in the answer too. R produces two different intervals, and they are not the same thing.

new_store <- data.frame(CompPrice = 125, Income = 70, Advertising = 10,
                         Price = 115, ShelveLoc = "Good", Age = 50)

predict(lm_fit, new_store, interval = "confidence")
#      fit      lwr      upr
# 1 10.879   10.634   11.125

predict(lm_fit, new_store, interval = "prediction")
#      fit      lwr      upr
# 1 10.879    8.811   12.948

The gap is more than eightfold, and the reason is mathematical: the confidence interval carries only the uncertainty in the estimated coefficients. The prediction interval adds the irreducible error on top — the store-specific randomness the model can never explain. Here the residual standard error is σ = 1.044, and the prediction interval's width of roughly 2 × 1.96 × 1.044 ≈ 4.09 comes straight from it.

5. Improving the Model

The first model has seven variables, while Carseats offers ten predictors. The missing ones — Population, Education, Urban and US — were left out deliberately: none of them was significant in the full model.

library(caret)

ctrl <- trainControl(method = "cv", number = 10)

cv_fit <- train(Sales ~ CompPrice + Income + Advertising + Price + ShelveLoc + Age,
                data = train,
                method = "lm",
                trControl = ctrl)

cv_fit$results

Cross-validation stops performance from hanging on a single random split. With 400 observations, the number a single test set reports can move noticeably with the chosen seed.

Running 10-fold cross-validation on this dataset gives the following picture:

One detail stands out: the RMSE we measured on the single test set was 0.934 — at the optimistic end of the cross-validated distribution. The best of the ten folds comes in at 0.833, the worst at 1.168, and the 0.34 between them is purely a matter of which observations landed in the test fold. "Test RMSE is 0.934" is not wrong, but it is the product of a lucky split.

Were those four excluded variables really useless?

"None of them was significant" is not a sufficient argument; variables that are individually insignificant can matter jointly. The right tool is a nested F test: since the smaller model is a restricted version of the larger one, do the four extra variables add explanatory power as a group?

full <- lm(Sales ~ ., data = train)
anova(lm_fit, full)

# Res.Df    RSS Df Sum of Sq      F Pr(>F)
#    312 340.02
#    308 334.82  4    5.2008 1.1961 0.3125

AIC(lm_fit, full)  # 945.54  vs  948.61

p = 0.313: the four variables add nothing as a group either. Adjusted R² rises by 0.0004 while AIC actually gets worse. Preferring the smaller model here is not a stylistic choice — it is what the data says.

6. Common Mistakes

Conclusion

Multiple linear regression is one of the simplest models in machine learning, but its interpretability keeps it valuable. Here the model explained 87% of the variation on the training set and 86% on the test set, with every assumption satisfied.

The real business output is in the coefficients: improving shelving location is a far stronger lever in this data than increasing the advertising budget. A model's worth lies not in the R² it reports but in the decision it supports.

“Fitting the model is one line; deciding whether to trust it is a whole chapter.”