Multiple Regression in Psychology: Tutorial with SPSS and R

Multiple regression is probably the most widely used multivariate analysis in psychology, and also one that requires the most methodological decisions before, during, and after its execution. One of the most controversial decisions is the variable selection method: should you include all predictors at once (enter method), add them one by one based on their statistical contribution (forward or stepwise method), or follow some other strategy? The answer depends fundamentally on whether your goal is to explain a phenomenon or to predict an outcome, and confusing these two purposes is a common source of problems.

The stepwise method was enormously popular for decades, but it has fallen out of favor among methodologists for well-founded reasons. In its classic version, the algorithm adds predictors to the model one at a time, selecting at each step the one that makes the largest statistically significant contribution, and removing those that cease to be significant when new predictors are included. The problem is that this process capitalizes on chance in a massive way. With 20 potential predictors and 100 participants, stepwise will find a model that appears to work well with your data, but that will probably not replicate in a new sample. Type I error rates are dramatically inflated, coefficients are biased toward more extreme values, and the R-squared of the final model overestimates the true explained variance.

Enter methods and theory-guided selection

The enter method (entering all predictors simultaneously) avoids the problems of stepwise because it does not rely on the data to decide which variables to include. The decision about which predictors to include is based on theory and prior evidence, not on the p-values of the current sample. This produces less biased and more replicable estimates, although it requires that you have a theoretical justification for each predictor you include. The drawback is that with many predictors and a moderate sample, standard errors can be large and multicollinearity can make individual coefficients difficult to interpret.

An intermediate approach is hierarchical regression, where you enter predictors in blocks according to a predetermined theoretical order. For example, in a study predicting academic performance, you might enter sociodemographic variables first (block 1), then cognitive variables (block 2), and finally motivational variables (block 3). At each step you evaluate whether the new block explains significant additional variance beyond the previous blocks. This approach is informative because it allows you to evaluate the incremental contribution of each set of variables, and the order of entry is determined by theory, not by the data.

Interpreting the coefficients

Interpreting multiple regression coefficients requires more caution than many researchers exercise. The partial regression coefficient (b or beta) of a predictor X represents the expected change in Y for each unit change in X, holding constant all other predictors in the model. This "holding constant the other variables" interpretation is crucial and is often forgotten. If self-esteem and optimism correlate .60 with each other and both predict well-being, the partial coefficients for each will be substantially attenuated relative to their bivariate correlations with well-being, because each one is controlling for the other.

Multicollinearity: that is, high correlation among predictors: does not strictly violate any assumption of regression, but it makes individual coefficients unstable and difficult to interpret. Variance inflation factors (VIF) are the most commonly used measure for detecting it: values above 5 or 10 (depending on the reference you consult) indicate problematic multicollinearity. Solutions may include removing one of the redundant predictors, combining them into a composite score, or using techniques such as ridge regression that penalize coefficients and produce more stable estimates at the cost of slight bias.

Beyond classic linear regression

Classic linear regression assumes a linear relationship between the predictors and the dependent variable, normally distributed residuals, homoscedasticity, and independence of observations. When the dependent variable is not continuous (it is dichotomous, ordinal, a count, or a proportion), you need generalized regression models. Logistic regression for binary dependent variables and Poisson regression for count variables are the most commonly used in psychology, but the general principles of predictor selection and coefficient interpretation are analogous.

If your data have a hierarchical or nested structure (students within classrooms, repeated measures within participants), classic linear regression is not appropriate because it violates the independence assumption. In these cases, multilevel or mixed models are the correct alternative. And if your primary interest is prediction (for example, classifying patients into diagnostic categories), machine learning techniques such as regularized regression (LASSO, elastic net) offer important advantages over traditional stepwise regression, including built-in cross-validation and better control of overfitting.

Diagnostics: residuals, outliers, and influential cases

A regression model is only as trustworthy as the diagnostics that accompany it. In my consulting work I find that researchers often report R-squared and the table of coefficients but skip the residual analysis entirely, which is precisely where most of the problems hide. Plotting residuals against fitted values is the single most informative diagnostic: a random cloud around zero suggests linearity and homoscedasticity are reasonable; a funnel shape signals heteroscedasticity; a curve signals non-linearity that the model is failing to capture. In R, plot(model) on a base lm() object gives you four diagnostic plots almost for free, and the performance package (from the easystats ecosystem) provides check_model(), which generates a complete visual assessment of all assumptions in a single call.

Influential cases deserve special attention because a single observation can dramatically alter your coefficients. Cook's distance, leverage values, and DFBETAS are the standard metrics. A common rule of thumb is to flag cases with Cook's distance above 4/n or with standardized residuals beyond ±3 for individual inspection. The right response is rarely "delete the outlier": the right response is to understand why that case is influential. Is it a data entry error? A genuinely unusual participant who falls outside the population of interest? Or simply tail behavior that your model fails to accommodate? Report what you find. Sensitivity analyses with and without influential cases give reviewers (and readers) confidence that your conclusions are not driven by a handful of observations. This is closely related to the broader topic of verifying statistical assumptions, which I cover in another article.

Interactions and centering: the part most papers get wrong

Interaction terms are where multiple regression earns its analytical power and where reporting most often goes off the rails. When you include X1 * X2 in a model, the main effects of X1 and X2 are no longer "average effects": they become conditional effects at the value zero of the other predictor. If zero is not a meaningful value (because your variables are scales from 1 to 7, for example), these "main effects" are uninterpretable, even though they are statistically valid. The standard fix is mean-centering the continuous predictors before computing the interaction, which makes the lower-order coefficients interpretable as effects at the mean of the moderator and also reduces non-essential multicollinearity between the product term and its components.

Once a significant interaction is detected, do not stop at the omnibus test. Probe it. Simple slopes analysis (computing the effect of X1 at ±1 SD of X2, for example) and the Johnson-Neyman technique (identifying the regions of the moderator where the effect is significant) are now standard in psychology. In R, the interactions package automates both. Reporting only the interaction p-value without showing the conditional effects leaves the reader unable to understand the substantive pattern. For a deeper treatment of when an interaction becomes a moderation or a mediation, see my article on mediation, moderation, and mechanisms.

Reporting in APA 7 format

APA 7 expects regression results to include unstandardized coefficients (B) with standard errors or 95% confidence intervals, standardized coefficients (beta) when the predictors are on different scales, the t or F statistic, exact p-values (not just "p < .05"), and an overall R-squared along with an effect size for the model (typically Cohen's f-squared for the model and semi-partial correlations or delta R-squared for individual predictors). Confidence intervals are no longer optional; reviewers at Q1 journals will ask for them if missing. A clean reporting sentence looks like: "Self-efficacy predicted academic performance, B = 0.34, 95% CI [0.18, 0.50], beta = .29, t(247) = 4.16, p < .001, controlling for prior GPA and study hours." Tables should follow the APA 7 template, and packages such as apaTables, gtsummary, or sjPlot::tab_model() generate them with minimal post-editing.

If the regression in your thesis involves multiple predictors and you need help checking assumptions or interpreting coefficients, my statistical consulting service reviews the model with you and delivers results ready to paste into your manuscript. More details on the doctoral thesis page.