Multiple Regression Analysis: Formulas and Solved Examples

Multiple regression analysis helps you study how an outcome relates to several predictors at once. For example, corn production may vary with fertilizer and insecticide use, while income may relate to education and agricultural employment. However, a useful analysis requires more than inserting numbers into an equation: you must also examine uncertainty, model assumptions, and the meaning of each coefficient.

This guide moves from the basic model to hand calculations, hypothesis tests, partial correlations, and matrix notation. Along the way, you will work through three small datasets and 22 exercises aligned with the topics in the supplied textbook pages. Finally, original numerical checks show why a model can fit extremely well while still leaving important questions unanswered.

Data provenance: The three datasets come from the user-supplied photographs of Chapter 7, “Multiple Regression Analysis,” pages 154–172. The visible pages do not establish the book’s complete bibliographic details. Therefore, this article identifies the tables directly, uses original explanations, and treats the observations as historical teaching data rather than current economic or agricultural evidence. All new statistics and charts represent calculations for this article.

What Is Multiple Regression Analysis?

Multiple regression estimates the relationship between one response variable and two or more predictors. Unlike a simple regression, it lets you describe the association with one predictor while accounting for the others in the model. Consequently, its slopes usually differ from the slopes in separate one-predictor regressions.

Here, Yᵢ is the response for observation i, β₀ is the intercept, and βⱼ is a population slope. Meanwhile, uᵢ captures influences outside the stated linear relationship. An estimated equation replaces the unknown coefficients with b₀, b₁, and the remaining fitted slopes.

Throughout this guide, p counts slopes and k counts all coefficients, including the intercept. However, some textbooks and software manuals use p for the total number of coefficients instead. Always check that convention before copying a degrees-of-freedom formula.

How to interpret a slope

In an additive model, b₁ describes the change in the fitted response for a one-unit increase in X₁ at a fixed value of X₂. For example, a coefficient of 0.65 means a fitted increase of 0.65 response units per additional unit of X₁. Nevertheless, that interpretation describes a conditional association unless the research design justifies a causal conclusion.

The phrase “holding other variables constant” refers to the model comparison. It does not mean the researcher literally held those variables constant in an experiment. Moreover, the comparison becomes less credible when the sample contains few observations with similar control variables but different values of the predictor of interest.

What “linear” actually means

A linear regression is linear in its coefficients. Therefore, it can include a squared predictor, a logarithm, or a product of predictors without leaving the linear-regression framework. For example, a model with X and X² still estimates coefficients that enter additively.

However, the marginal relationship now depends on X: the derivative is β₁ + 2β₂X. Similarly, an interaction term changes the meaning of the main-effect coefficients. As a result, the constant-slope interpretation applies to the additive specification, not automatically to every model fitted by OLS.

Assumptions: separate estimation from inference

Under full rank and a zero conditional error mean, OLS has the familiar unbiasedness property. With a spherical error covariance, the Gauss–Markov theorem also makes OLS best among linear unbiased estimators. Normality adds exact small-sample inference; it is not a requirement for calculating the least-squares coefficients.

For the corn data, observations run across years. Therefore, serial dependence and shared trends deserve attention before anyone treats the classical p-values as definitive agricultural evidence. With only ten observations, diagnostic tests also have limited power.

Reference: Penn State STAT 501: MLR Estimation, Prediction & Model Assumptions.

Essential Multiple Regression Formulas

OLS chooses the smallest sum of squared residuals

Squaring gives larger mistakes more weight and prevents positive and negative residuals from canceling. Consequently, the minimizing coefficients balance all observations simultaneously. The resulting residuals sum to zero when the model includes an intercept, apart from numerical precision.

Centered sums for two predictors

Define yᵢ = Yᵢ − Ȳ, x₁ᵢ = X₁ᵢ − X̄₁, and x₂ᵢ = X₂ᵢ − X̄₂. Next, write S₁₁ = ∑x₁ᵢ², S₂₂ = ∑x₂ᵢ², S₁₂ = ∑x₁ᵢx₂ᵢ, S₁y = ∑x₁ᵢyᵢ, and S₂y = ∑x₂ᵢyᵢ. These quantities compress the data into the information needed for the slopes.

For nonconstant predictors, D equals zero under perfect collinearity. In that case, the data cannot identify the two slopes separately. However, a small positive D can still allow a numerical solution while making separate slope estimates imprecise.

Variance, standard errors, and confidence intervals

The symbol s² estimates the error variance, while SE(bⱼ) estimates uncertainty in a coefficient. Thus, these quantities answer different questions. A confidence interval describes uncertainty about a population coefficient, rather than the spread of individual responses.

R², adjusted R², and the overall F test

These identities assume ordinary least squares with an intercept and a nonconstant response. Here, SSR means explained sum of squares and SSE means residual sum of squares. However, notation varies across books, so inspect the definitions rather than relying on an abbreviation alone.

R² measures in-sample fit around the response mean. Adding predictors to a nested model on the same observations cannot lower ordinary R², because the larger model can reproduce the smaller model. By contrast, adjusted R² penalizes the loss of residual degrees of freedom.

Reference: Penn State STAT 501: Multiple Linear Regression.

Matrix notation

Each row of X contains an intercept entry of 1 followed by that observation’s predictor values. The matrix expression extends the same OLS calculation to any number of predictors. In actual numerical work, however, use a least-squares solver rather than explicitly forming an inverse to obtain b.

Reference: NumPy documentation: numpy.linalg.lstsq.

Worked Example: Corn, Fertilizer, and Insecticide

The first dataset contains ten observations labeled 1971–1980. In the supplied example, Y measures corn production in bushels per acre, while X₁ and X₂ measure fertilizer and insecticide in pounds per acre. Therefore, the fitted slopes use bushels per acre per additional pound per acre of the relevant input.

Step 1: calculate the centered sums

Step 2: write and interpret the fitted equation

At a fixed insecticide level, an additional pound per acre of fertilizer corresponds to about 0.6501 extra bushels per acre in the fitted model. Likewise, the insecticide slope corresponds to about 1.1099 extra bushels per acre at a fixed fertilizer level. However, neither estimate alone establishes what a real farm would gain from changing its inputs.

The intercept predicts about 31.98 bushels per acre when both inputs equal zero. Yet the sample contains no such combination. Therefore, treat the intercept primarily as part of the fitted equation, rather than as a validated no-input yield estimate.

Step 3: evaluate fit and coefficient uncertainty

Both slope p-values fall below 0.05 under the classical assumptions. Moreover, the joint test rejects the hypothesis that both slopes equal zero. Nevertheless, the exceptional in-sample fit does not validate the assumptions or guarantee accurate forecasts outside the observed input combinations.

Step 4: predict a response and distinguish two intervals

Suppose fertilizer equals 20 and insecticide equals 12 pounds per acre. Substituting those values gives a fitted yield of 58.3001 bushels per acre. Under the classical model, the interval for the mean response differs from the interval for a new individual response.

The prediction interval is wider because it includes both coefficient uncertainty and the new response’s unexplained variation. By comparison, the mean-response interval addresses uncertainty about the conditional mean alone. Both intervals remain conditional on the model assumptions and the relevance of the new input combination.

Reference: Penn State STAT 501: mean-response and prediction intervals.

Worked Example: Income, Agriculture, and Schooling

The second dataset contains 15 country observations in a historical teaching example labeled 1981. Here, Y represents real income per capita in thousands of U.S. dollars, X₁ is the agricultural share of the labor force in percent, and X₂ is average schooling in years for people over 25. Accordingly, one unit of X₁ means one percentage point, not a 1% proportional increase.

At a fixed schooling level, one extra percentage point in agricultural employment corresponds to about $376.16 lower fitted income per capita. Meanwhile, an extra year of average schooling corresponds to about $452.51 higher fitted income at a fixed agricultural share. These are country-level associations, so they do not directly measure an individual worker’s return to education.

The model explains about 69.32% of the sample variation in income around its mean. Furthermore, adjusted R² equals 0.642070 and the overall F statistic equals 13.556904 with 2 and 12 degrees of freedom. The joint p-value is about 0.000834, while both individual slope tests reject zero at the 5% level.

Worked Example: Historical Latin American Income Data

Our third example combines GDP per capita, agriculture’s share of the economy, and male literacy for 15 Latin American countries. According to the supplied table, its historical inputs trace back to World Bank World Development Indicators for 1999. However, the calculations here reproduce the photographed table, not a new download from the World Bank or a claim about present-day country conditions.

In this dataset, Y records GDP per capita in units of $100, X₁ records the agricultural share in percent, and X₂ records male literacy in percent. Thus, a slope of −1.9483 response units means approximately −$194.83 per percentage point. The units differ from the schooling example, even though both responses concern income.

Agriculture has a negative fitted slope, while literacy has a positive slope. However, their individual p-values are approximately 0.1170 and 0.4836, so neither test rejects zero at 5%. A failure to reject does not establish that either population relationship equals zero.

Meanwhile, the overall test gives F(2, 12) = 6.771170 and p ≈ 0.010753. Therefore, the model rejects the joint hypothesis that both slopes are zero, despite the two individually nonsignificant tests. This dataset provides a concrete example of why joint and individual tests answer different questions.

R² equals 0.530192 and adjusted R² equals 0.451890. In addition, the coefficient intervals remain wide, reflecting the limited information available for separate associations. Small samples, predictor relationships, and unexplained variation all matter; this pattern does not uniquely diagnose severe multicollinearity.

22 Solved Exercises in Multiple Regression Analysis

The following exercises follow the sequence of topics visible in the supplied chapter. Each solution uses original wording, explicit notation, and recalculated results where a dataset applies. Because the calculations retain full precision, some final numbers differ from the rounded printed answers.

Exercise 7.1. Specify the model and assumptions

The two-predictor specification has three coefficients. Therefore, its residual degrees of freedom equal n − 3. Full rank and E(u | X) = 0 support unbiased estimation; constant conditional variance and uncorrelated errors support the classical covariance calculation, while normality gives exact finite-sample t and F tests.

Exercise 7.2. Derive the normal equations

Start with SSE = ∑eᵢ² and differentiate with respect to each coefficient. Next, set all three derivatives equal to zero. These first-order conditions imply that residuals are orthogonal to the intercept and to each predictor.

With full column rank, the quadratic objective has a unique minimum. Otherwise, more than one coefficient vector may produce the same fitted values. Therefore, the normal equations alone do not guarantee unique coefficients.

Exercise 7.3. Obtain the centered slope formulas

Subtract the mean equation from each fitted equation to remove the intercept. As a result, the two slope equations become S₁y = b₁S₁₁ + b₂S₁₂ and S₂y = b₁S₁₂ + b₂S₂₂. Solve this pair simultaneously to obtain the formulas in the earlier formula section.

Finally, recover b₀ from the sample means. Notice that n = 3 can identify three coefficients if X has full rank, but it leaves no residual degrees of freedom. Consequently, you cannot estimate the error variance in the usual way with that saturated model.

Exercise 7.4. Interpret the coefficients and BLUE

The intercept gives the fitted response when all predictors equal zero, while each additive slope describes a conditional fitted change. Under the Gauss–Markov conditions, OLS is the best linear unbiased estimator, or BLUE. Here, “best” means minimum covariance in the relevant positive-semidefinite comparison among linear unbiased estimators, not superiority over every possible estimator or prediction method.

Exercise 7.5. Estimate the income–schooling model

Therefore, the fitted equation is ŷ = 6.202980 − 0.376164X₁ + 0.452514X₂. Its slopes describe lower fitted income with a higher agricultural labor share and higher fitted income with more schooling, conditional on the other predictor. However, the country-level data do not establish causal effects.

Exercise 7.6. Estimate the Latin American GDP model

Insert those centered sums into the two-slope formulas. The resulting estimates are b₀ = 5.030977, b₁ = −1.948288, and b₂ = 0.530123. For comparison, regressing GDP on agriculture alone gives −1,149 / 442 = −2.599548, so the conditional agricultural slope is less negative.

This change illustrates the effect of adding a correlated predictor. However, it does not prove that the longer model has identified a causal agricultural effect. Also, the intercept corresponds to zero agriculture and zero male literacy, a combination far outside these observations.

Exercise 7.7. Distinguish error variance, coefficient variance, and standard error

The intercept’s variance is the first diagonal entry of s²(XᵀX)⁻¹. Although an intercept test may have little substantive value when zero predictor values are unrealistic, its uncertainty still matters for some predictions and linear combinations. Therefore, you should not assume that the intercept is inherently unimportant.

Exercise 7.8. Calculate standard errors for the schooling model

First calculate the variances at full precision, and then take their square roots. Otherwise, rounding 0.01428 to 0.01 before taking the square root produces 0.10 instead of 0.11951. That change materially alters the t statistic.

Exercise 7.9. Test both schooling-model slopes

Both absolute t statistics exceed 2.178813, so both tests reject their zero-slope null hypotheses at 5%. Moreover, the exact two-sided p-values are 0.015058 and 0.002593. Thus, schooling also passes a 1% test, while agriculture does not.

Exercise 7.10. Construct 95% coefficient intervals

Neither interval includes zero, which agrees with the two-sided 5% tests. In dollar units, the schooling interval corresponds to approximately $192.12 to $712.91 per additional year of average schooling. However, its interpretation remains conditional and country-level.

Exercise 7.11. Calculate uncertainty for the Latin American GDP model

These standard errors are large relative to the slope estimates. Consequently, the model supplies limited precision for separating the agricultural and literacy associations. The fitted response itself also contains substantial unexplained variation.

Exercise 7.12. Test the GDP-model slopes

Neither absolute statistic exceeds the 5% critical value of 2.178813. Therefore, fail to reject both zero-slope null hypotheses individually. This wording preserves the distinction between insufficient evidence against zero and evidence that the true association must equal zero.

Exercise 7.13. Construct GDP-model confidence intervals

Both intervals contain zero, consistent with the individual tests. Since Y uses $100 units, the agricultural interval corresponds to approximately −$446.15 to $56.50 per percentage point. Meanwhile, the literacy interval corresponds to approximately −$106.76 to $212.78 per percentage point.

Exercise 7.14. Derive R² from the centered products

With an intercept, fitted deviations and residuals are orthogonal. Therefore, SST equals SSR plus SSE. Also, the centered normal equations imply that SSR = b₁S₁y + b₂S₂y.

Use unrounded OLS coefficients when checking this identity. Otherwise, a coefficient vector that no longer exactly satisfies the normal equations can produce slightly inconsistent totals. That numerical mismatch does not invalidate the algebra.

Exercise 7.15. Calculate schooling-model R² in three ways

All three methods agree when the calculation retains precision. Therefore, about 69.32% of the observed income variation around the sample mean appears in the fitted values. This percentage describes fit, not the proportion of income that schooling and agriculture causally create.

Exercise 7.16. Derive adjusted R² and discuss its range

For a fixed n and k with an intercept, nonconstant Y, and n > k, adjusted R² cannot exceed 1. Its minimum is 1 − (n − 1)/(n − k), which is negative when k > 1. Unlike ordinary in-sample R², it can fall below zero because it compares variance estimates after the degrees-of-freedom adjustment.

Exercise 7.17. Calculate adjusted R² for the schooling model

The agriculture-only model has R² = 0.326666667 and adjusted R² = 0.274871795. After adding schooling, both measures increase substantially. Nevertheless, model selection should also consider the research question, diagnostics, and performance on new observations.

Exercise 7.18. Recover SSE without calculating every fitted response

This shortcut avoids a separate fitted-value calculation for each observation. However, it does not remove the need to inspect individual residuals when evaluating the model. A correct aggregate SSE can coexist with an influential observation or systematic residual pattern.

Exercise 7.19. State and perform the overall significance test

The numerator measures explained variation per slope, and the denominator estimates unexplained error variance. Under the classical normal-error null, their ratio follows F with p and n − k degrees of freedom. Importantly, the null excludes the intercept.

Exercise 7.20. Compare a simple-regression F test with individual tests

With one slope, the overall F test and the corresponding two-sided t test give the same decision. However, with several slopes, the overall test can reject while every individual test fails to reject. The historical Latin American dataset demonstrates that pattern numerically.

Exercise 7.21. Express F using R²

Substitute the two identities into the F ratio and cancel SST. As a result, the overall test can use R², n, and the number of slopes without reconstructing the full ANOVA table. Under the stated assumptions, it tests whether the population slopes are jointly zero, rather than asking whether sample R² is merely numerically positive.

Exercise 7.22. Test the schooling model jointly

Since 13.556904 exceeds 3.885294, reject the hypothesis that both slopes equal zero. Therefore, the two predictors jointly improve fit relative to the intercept-only model under the classical test. The result still does not establish causation or guarantee useful forecasts.

Original Insights Beyond the Basic Formulas

1. A rounding audit explains why software may disagree with a worked answer

A useful solution should allow readers to reconcile the hand calculation with software output. Therefore, the table below compares selected numbers printed in the supplied pages with calculations from the underlying observations. It distinguishes harmless coefficient display rounding from larger changes caused by rounding intermediate variances or means.

These differences do not justify dismissing the entire example. Instead, they show why you should preserve precision until the final reporting step. In particular, calculate a square root before rounding the variance beneath it.

2. Excellent fit can coexist with strong multicollinearity

Here, Rⱼ² comes from regressing predictor j on the other predictors, not from the main response regression. A VIF of about 18.46 means that the classical slope variance is 18.46 times its orthogonal-predictor benchmark at the same error variance and centered predictor sum of squares. Consequently, the corresponding standard-error factor is √18.46, or about 4.30.

Reference: Penn State: Detecting Multicollinearity Using Variance Inflation Factors.

The corn model can still estimate the joint direction of increasing inputs very accurately within this sample. However, it has less information about changing one input while fixing the other, because the two inputs move together. That distinction explains why fit and coefficient precision can tell different stories.

3. Partial correlation measures what remains after adjustment

To isolate the linear relationship between Y and X₁ after accounting for X₂, remove the fitted X₂ component from both variables. Then correlate the two residual series. This gives a partial correlation, which differs from the ordinary correlation between the original variables.

For the corn model, the fertilizer partial correlation is approximately 0.700711 and the insecticide partial correlation is approximately 0.843220. Their squared values describe the share of remaining variation associated with the additional predictor after the other predictor enters. However, a larger partial correlation does not establish greater causal importance or economic value.

4. Partial R² and the partial F test answer a practical model-comparison question

Suppose you already use agricultural employment to predict income and want to add schooling. The restricted model contains agriculture and an intercept, while the full model contains both predictors. Therefore, the comparison asks whether schooling adds information conditional on agriculture.

For the schooling example, the restricted SSE equals 26.933333 and the full SSE equals 12.271881. Thus, the partial R² is approximately 0.544361: schooling accounts for about 54.44% of the variation left unexplained by the agriculture-only model. With one added coefficient, the partial F equals the schooling t statistic squared, approximately 14.336631.

Reference: Penn State STAT 501: general linear F tests and partial R-squared.

5. Adjusted R² can increase without 5% significance

Compare nested models on identical observations, with the full model adding one predictor. Let ν denote the full model’s residual degrees of freedom. The restricted model then has ν + 1 degrees of freedom.

A 5% two-sided t test usually requires a threshold well above 1. Consequently, adjusted R² can favor a variable that does not meet the conventional 5% significance threshold. This is an algebraic consequence of the adjustment, not a software inconsistency.

However, do not treat either threshold as an automatic rule for keeping or deleting controls. A variable may matter for the research design even when its individual p-value is large. Likewise, repeated selection based on the same sample changes the interpretation of later confidence intervals and tests.

6. Leave-one-out prediction tells a less flattering story

The same observations used to estimate a model also make its training fit look favorable. To examine that optimism, leave out one observation, fit on the others, and predict the omitted response. Repeat this process for every observation and sum the squared prediction errors.

The displayed shortcut applies to these full-rank OLS fits with an intercept. It reproduces explicit leave-one-out refitting while avoiding a new fit for each omitted row. Nevertheless, its interpretation depends on whether leaving out isolated rows resembles the prediction task.

The Latin American example falls from training R² of about 0.53 to predicted R² of about 0.17 under this definition. Therefore, the fitted relationship predicts omitted countries much less well than it describes the training sample. With only 15 countries, however, this diagnostic also remains sensitive to individual observations.

For corn, leave-one-out prediction mixes earlier and later years. As a result, it should not substitute for a chronological forecasting evaluation. A genuine forecast study would train on past data and evaluate later observations without using future information.

7. Centering changes the intercept, not the underlying fit

Replace fertilizer by X₁ − 18 and insecticide by X₂ − 12. The corn model then has an intercept of exactly 57, because the centered predictors both equal zero at their sample means. Meanwhile, its slopes, fitted values, residuals, and R² remain unchanged.

This makes the intercept easier to interpret as fitted yield at the average inputs. However, centering these two additive predictors does not eliminate their correlation or reduce their VIF. It can help with some numerical or polynomial-term issues, but it cannot create independent information that the sample lacks.

Diagnostics Before You Trust a Regression

A strong model summary begins a diagnostic review rather than ending it. NIST emphasizes residual analysis because a high R² alone does not guarantee an adequate model. Therefore, inspect several aspects of the fit before interpreting coefficients or making predictions.

Reference: NIST/SEMATECH e-Handbook: How can I tell if a model fits my data?.

An influential observation is not automatically an error. Therefore, investigate its provenance before considering exclusion. Report how reasonable analysis choices change the result, rather than silently deleting inconvenient observations.

Reference: Penn State STAT 501: Influential Points.

What robust standard errors can and cannot repair

Heteroskedasticity-robust covariance estimators change the uncertainty calculation while leaving OLS coefficients unchanged. For example, HC3 adjusts residual contributions using leverage. However, it does not repair omitted-variable bias, reverse causality, an inappropriate mean function, or arbitrary serial dependence.

Reference: statsmodels documentation: HC3 covariance.

With small datasets, robust methods do not create missing information. Consequently, a careful report should describe both the estimator and the limitations of the sample. Clustered or time-dependent data require covariance choices that match their dependence structure.

A practical sequence for an applied project

1. Define the question. Decide whether you want description, prediction, or a causal estimate.
2. Audit the variables. Record units, dates, missing values, transformations, and population coverage.
3. Specify the model. Choose predictors and controls using subject knowledge and the target question.
4. Fit at full precision. Use a stable least-squares routine and confirm the design matrix has full rank.
5. Read the coefficients. Translate changes into meaningful units before interpreting p-values.
6. Inspect uncertainty and diagnostics. Evaluate intervals, residual patterns, leverage, and predictor relationships.
7. Validate the intended use. Match the validation split to future observations, new groups, or new time periods.
8. Report transparently. Present the specification, sample size, uncertainty method, and limitations.

Reproduce the Calculations in Python

The following script contains all three datasets and reproduces the main results using NumPy and SciPy. First, it fits coefficients with a least-squares solver. Next, it calculates classical standard errors, tests, confidence intervals, and leave-one-out diagnostics.

import numpy as np
from scipy.stats import t, f

data = {
    "corn": [
        [40,44,46,48,52,58,60,68,74,80],
        [6,10,12,14,16,18,22,24,26,32],
        [4,4,5,7,9,12,14,20,21,24]],
    "schooling": [
        [6,8,8,7,7,12,9,8,9,10,10,11,9,10,11],
        [9,10,8,7,10,4,5,5,6,8,7,4,9,5,8],
        [8,13,11,10,12,16,10,10,12,14,12,16,14,10,12]],
    "latin_gdp": [
        [76,10,44,47,23,19,13,19,8,44,4,31,24,59,37],
        [6,16,9,8,14,11,12,10,18,5,26,8,8,9,5],
        [97,92,85,96,91,83,93,81,74,93,67,92,94,97,93]]
}

for name, values in data.items():
    y, x1, x2 = np.asarray(values, dtype=float)
    X = np.column_stack([np.ones(len(y)), x1, x2])
    n, k = X.shape
    assert np.linalg.matrix_rank(X) == k
    b = np.linalg.lstsq(X, y, rcond=None)[0]
    fitted = X @ b
    e = y - fitted
    sse = e @ e
    sst = (y - y.mean()) @ (y - y.mean())
    df = n - k
    mse = sse / df
    # QR supplies the covariance geometry without inverting X.T @ X.
    Q, upper = np.linalg.qr(X, mode="reduced")
    upper_inv = np.linalg.solve(upper, np.eye(k))
    geometry = upper_inv @ upper_inv.T
    se = np.sqrt(np.diag(mse * geometry))
    t_stat = b / se
    p_value = 2 * t.sf(np.abs(t_stat), df)
    critical = t.ppf(0.975, df)
    ci = np.column_stack([b-critical*se, b+critical*se])
    r2 = 1 - sse/sst
    adj_r2 = 1 - (sse/df)/(sst/(n-1))
    F = ((sst-sse)/(k-1))/mse
    leverage = np.sum(Q**2, axis=1)
    press = np.sum((e/(1-leverage))**2)
    vif = 1/(1-np.corrcoef(x1, x2)[0,1]**2)
    print(name)
    print("b, SE, t, p, CI:", b, se, t_stat, p_value, ci)
    print("SSE, R2, adjusted R2:", sse, r2, adj_r2)
    print("F, F p-value:", F, f.sf(F, k-1, df))
    print("VIF, PRESS, predicted R2:", vif, press, 1-press/sst)
    assert np.allclose(X.T @ e, 0, atol=1e-8)

The classical p-values in this script assume the conventional normal, homoskedastic, uncorrelated-error model. Therefore, matching the published arithmetic does not by itself validate those assumptions. The code intentionally keeps estimation and interpretation separate.

Frequently Asked Questions

What is the difference between multiple and multivariate regression?

Multiple regression has several predictors for one response. In contrast, multivariate regression models multiple response variables together. The number of response variables, rather than the number of predictors, distinguishes the two terms.

How many observations does multiple regression require?

A unique OLS solution requires a full-rank design matrix, and classical residual-variance estimation requires n greater than k. However, that algebraic minimum is not a practical sample-size recommendation. The needed sample depends on the effect sizes, noise, predictor relationships, design, and desired precision.

Do all predictors need normal distributions?

No. Normality in the classical finite-sample inference model concerns errors conditional on the predictors. Therefore, skewed predictors do not automatically invalidate multiple regression.

Can R² be negative?

In-sample OLS with an intercept and a nonconstant response produces R² between zero and one. However, adjusted R² and prediction-based R² can be negative. Some no-intercept definitions also produce values outside the familiar range.

Does a significant F test mean every predictor is significant?

No. The overall test evaluates the slopes jointly. Consequently, it can reject the all-zero null even when individual slope tests do not reject, as the Latin American example shows.

Should I remove every predictor with p greater than 0.05?

No. A predictor may serve as an important control or form part of a prespecified research model. Moreover, repeated selection based on p-values can make later uncertainty statements misleading. Consider the question, design, diagnostics, and validation evidence together.

Is a coefficient with a larger absolute value more important?

Not necessarily. Coefficients depend on measurement units, so comparing a dollar slope with a percentage-point slope can mislead. Furthermore, statistical precision, predictive contribution, and causal importance are distinct concepts.

Why do my software results differ from the printed solution?

Check the data transcription, intercept setting, missing-data sample, and variance estimator first. Next, preserve precision throughout intermediate calculations. The rounding audit in this article shows how premature rounding changes several textbook-style answers.

Can multiple regression establish causation?

A regression equation alone cannot establish causation. However, a credible experimental or identification design may support a causal interpretation under explicit assumptions. Adding controls does not automatically provide such a design.

Conclusion: Read the Whole Model

Multiple regression analysis becomes useful when you connect its equations to the data and the question. The three examples show how to estimate slopes, calculate standard errors, test individual and joint hypotheses, and interpret fitted relationships in the correct units. Moreover, the recalculations demonstrate why numerical precision matters even in small teaching datasets.

The most revealing results come from comparing different measures. A 99.16% R² can coexist with a VIF above 18; a significant overall F test can coexist with nonsignificant individual slopes. Therefore, use fit, uncertainty, diagnostics, and research design together before making a substantive claim.

References and Data Notes

All dataset-specific values, confidence intervals, charts, and diagnostic comparisons are original recalculations from the observations reproduced above. The source photographs supply the historical examples and exercise topics. For additional explanations and implementation details, consult the following institutional and official references.

1. Penn State STAT 501: Multiple Linear Regression
   https://online.stat.psu.edu/stat501/Lesson05
2. Penn State STAT 501: MLR Model Evaluation
   https://online.stat.psu.edu/stat501/Lesson06
3. Penn State STAT 501: MLR Estimation, Prediction & Model Assumptions
   https://online.stat.psu.edu/stat501/Lesson07
4. Penn State: Detecting Multicollinearity Using Variance Inflation Factors
   https://online.stat.psu.edu/stat462/node/180/
5. Penn State STAT 501: Influential Points
   https://online.stat.psu.edu/stat501/Lesson11
6. NIST/SEMATECH: How can I tell if a model fits my data?
   https://www.itl.nist.gov/div898/handbook/pmd/section4/pmd44.htm
7. NumPy: numpy.linalg.lstsq
   https://numpy.org/doc/stable/reference/generated/numpy.linalg.lstsq.html
8. statsmodels: HC3 covariance
   https://www.statsmodels.org/stable/generated/statsmodels.stats.sandwich_covariance.cov_hc3.html