Multiple Linear Regression

Models a continuous outcome from several predictors at once, giving each predictor's unique contribution — the workhorse of survey-based management research.

Data requirements

Assumptions

• Linearity of each predictor–outcome relation (partial regression plots).
• Independence of residuals (Durbin–Watson ≈ 2 for ordered data).
• Homoscedasticity of residuals.
• Normality of residuals.
• No severe multicollinearity (VIF < 5–10; tolerance > .1–.2).
• No unduly influential cases (Cook's distance, leverage).

Hypotheses

The concept

OLS finds the weights that minimize squared prediction errors using all predictors jointly. Each coefficient bj is the expected change in Y for a one-unit change in Xj holding the other predictors constant — a 'unique contribution' logic that makes multiple regression the standard tool for statistical control.

Standardized betas allow rough comparison of predictor importance on a common scale; squared semi-partial correlations give each predictor's unique slice of R². Adjusted R² corrects the optimism of adding predictors. Two different projects — explanation (theory-driven, enter variables by design) and prediction (maximize out-of-sample accuracy, use validation) — should not be mixed casually; avoid mechanical stepwise selection for explanatory work.

Worked example

Predicting work engagement (1–7) from autonomy, supervisor feedback, and tenure among 214 employees.

Result: R² = .31, F(3, 210) = 31.4, p < .001. Autonomy β = .38 (p < .001) and feedback β = .21 (p = .002) are unique predictors; tenure β = .06 (p = .34) adds nothing beyond them.

How to run it

model <- lm(engagement ~ autonomy + feedback + tenure, data = df)
summary(model)
confint(model)

car::vif(model)                     # multicollinearity
par(mfrow = c(2,2)); plot(model)    # diagnostics

library(effectsize)
standardize_parameters(model)       # standardized betas

import statsmodels.formula.api as smf
from statsmodels.stats.outliers_influence import variance_inflation_factor

model = smf.ols("engagement ~ autonomy + feedback + tenure", data=df).fit()
print(model.summary())

X = df[["autonomy", "feedback", "tenure"]].assign(const=1)
for i, c in enumerate(X.columns):
    print(c, variance_inflation_factor(X.values, i))

Interpreting the output

• Model level: R², adjusted R², F, p — does the set predict at all?
• Each b in raw units ('one more feedback point ≈ .21 engagement points, holding others constant'); β for relative weight.
• CIs and p per predictor; non-significance means no unique contribution — the variable may still correlate bivariately.
• Diagnostics (VIF, residual plots) are part of the result, not an optional extra.

APA-style reporting

Common mistakes

Related methods