Multiple Regression Exercises in R: 17 Real-World Practice Problems

Explanation: The formula mpg ~ wt + hp tells lm() to estimate one intercept and two partial slopes. coef() returns a named numeric vector, which is the cleanest way to grab a single slope by name later (e.g. coef(ex_1_1)["wt"]). Each slope is conditional: the wt slope is the effect of weight after holding horsepower fixed, not the marginal effect of weight on its own.

Explanation: A partial slope is the predicted change in the outcome per unit change in that predictor, holding the others fixed. Multiplying the slope by the size of the move (0.5) gives the predicted change directly. unname() strips the carried-over wt name so the result is a clean scalar, which matters when you feed numbers into reporting tables.

Explanation: confint() returns a matrix with one row per coefficient. Passing parm = "hp" limits the rows it computes, and [1, ] collapses the single row to a named vector. The interval excludes zero, which is the geometric statement of "p < 0.05" for that coefficient. Pair the slope with its interval whenever you report a regression result, because a slope without a width is incomplete information.

Explanation: Always read summary values programmatically, never from the printout. summary(fit) returns a list with named elements like adj.r.squared, r.squared, coefficients, fstatistic, and sigma. The adjusted version penalises model size, so it goes up only when an added predictor pulls its statistical weight. Compare adj-R² across nested models to decide whether a new variable is worth keeping.

Explanation: summary(fit)$coefficients is a 4-column matrix (Estimate, Std. Error, t value, Pr(>|t|)) indexed by predictor name. Indexing by name in both dimensions is robust to predictor reordering, which is why this pattern is safer than [3, 3]. The t-value here is borderline, which is exactly the kind of "should I keep this predictor?" signal you want to make a deliberate call on rather than eyeballing stars.

Explanation: anova() on two nested lm objects runs a partial F-test: does the extra block of predictors explain a statistically meaningful chunk of residual variance beyond the smaller model? The p-value is on the second row because the first row is the "smaller model" baseline. Use this test, not a stack of individual t-tests, whenever you want to add or drop a group of predictors together.

Explanation: step() greedily drops one predictor at a time, recomputes AIC, and stops when no single drop lowers AIC further. direction = "backward" starts from the full model and prunes; "both" lets it add variables back. AIC trades fit against complexity (AIC = -2 log L + 2 k), so unlike adjusted R² it imposes a hard penalty per parameter. Stepwise selection is a screening tool, not a final answer; always cross-check the chosen model on a holdout slice.

Explanation: fitted() returns the predicted values at each training row and resid() returns observed minus predicted. Stacking them into a tibble is the gateway to every regression diagnostic plot: residuals-vs-fitted (heteroscedasticity), QQ plot (normality), scale-location (variance trend), residuals-vs-leverage (influence). If the residuals fan out as the fitted value grows, your variance is not constant and a log transform or robust standard errors are on the table.

Explanation: Shapiro-Wilk is a formal test of "the data are drawn from a normal distribution." With n = 32 it has limited power, so failing to reject (p > 0.05) does not mean "the residuals are perfectly normal," it means you do not have strong evidence against normality. Always pair the test with a QQ plot; the eyeball test catches tail issues that small-sample tests miss.

Explanation: Breusch-Pagan regresses the squared residuals on the original predictors and tests whether any predictor has explanatory power over the variance. A p of 0.09 sits in the "smoke but not flame" zone: there is mild evidence the variance is not flat. Common fixes are to log-transform the outcome, use weighted least squares, or compute heteroscedasticity-consistent (HC3) standard errors via sandwich::vcovHC().

Explanation: Cook's distance measures how much every fitted value would shift if you deleted one row. The 4/n rule of thumb (0.125 here) flags candidates to inspect, not rows to delete blindly. Always look at flagged rows in context: the Maserati Bora has unusually high horsepower for its weight, so it is an honest outlier, not a data error. Removing it would bias the slope estimate toward the bulk of the sample.

Explanation: A correlation matrix is the first multicollinearity check: anything above 0.8 in absolute value is a warning sign. Here disp is strongly correlated with wt, cyl, and hp, which means putting all four in a model is likely to inflate standard errors and flip coefficient signs. The matrix is necessary but not sufficient; even with all pairwise correlations below 0.8, three predictors can still co-move enough to inflate VIF.

Explanation: VIF for predictor j is 1 / (1 - R_j^2), where R_j^2 is the R-squared of regressing predictor j on every other predictor. A VIF of 8.8 for cyl means its standard error is sqrt(8.8) = 2.97 times larger than it would be if cyl were orthogonal to the rest. The fix is usually to drop one of the redundant predictors, combine them into a single derived variable, or shrink coefficients with ridge or lasso.

Explanation: Dropping cyl removes the most redundant predictor; the remaining three are pairwise correlated but no longer flagged. This is the typical VIF workflow: identify the worst offender, drop it (or merge with a sibling), and confirm the rest fall under the threshold. If multiple predictors had VIF above 5, drop them one at a time and re-check, because dropping one often pulls the others below threshold on its own.

Explanation: wt * hp expands to wt + hp + wt:hp, so you get both main effects and the interaction in one stroke. The positive interaction coefficient means the negative effect of weight on mpg shrinks as horsepower grows, which makes physical sense: heavy cars with weak engines are penalised most. Always include main effects when including an interaction; dropping a main effect forces it through zero and silently biases the interaction.

Explanation: A log transform on the outcome turns additive effects into multiplicative ones; each coefficient becomes "approximate percentage change in mpg per unit predictor." The fit improves here because mpg has a long right tail and shrinks when logged, which evens out the residual variance. Watch out: R-squared on log(mpg) is not directly comparable to R-squared on mpg, because they live on different scales. Compare on the same metric (e.g. residual standard error on the original scale after back-transforming predictions).

Explanation: poly(hp, 2) builds two columns: a linear term and a quadratic term, both orthogonal to each other. That orthogonality matters because raw hp and hp^2 are correlated by construction (VIF in the dozens), which destabilises the coefficient estimates even when the fit is fine. The adj-R² climbs from 0.815 to 0.844, suggesting the mpg-vs-hp relationship is genuinely curved, not linear.

Explanation: Comparing nested models with anova() is the principled way to decide whether the interaction is worth its degree of freedom. The F-statistic of 14.8 with p well under 0.001 says the interaction does pay rent; the residual sum of squares drops sharply when you allow the wt slope to flex with hp. Without this test you would be guessing from the t-stat of the interaction alone, which works here but breaks for multi-degree-of-freedom factor interactions.

Explanation: interval = "confidence" returns the interval for the mean prediction at that x. interval = "prediction" would be wider because it also accounts for the residual noise of a single new observation. Use the confidence interval when reporting "what does the model say the expected mpg is for cars like this?" and the prediction interval when you need to bound an actual future observation, like the next car off the line.

Explanation: Holdout RMSE is the honest scorecard for a regression model; in-sample R-squared always overstates performance because the model has already seen those rows. With only 8 test rows the estimate is noisy, so in practice repeat the split many times or use k-fold cross-validation via caret::trainControl() or rsample::vfold_cv(). The takeaway here is the pattern: fit on train, predict on test, score with a metric that lives on the outcome's scale.

Explanation: broom::tidy() turns a model object into a flat tibble that plays nicely with dplyr, ggplot2, and report generation. conf.int = TRUE adds confidence bounds in two extra columns. This format is the right shape for forest plots, regression tables in Quarto reports, and bulk comparison across many models (purrr::map(fits, tidy) |> bind_rows(.id = "model")).