Fitting a Multiple Linear Regression

Concepts

• A response often depends on more than one predictor, and using only one can leave out information that actually matters for explaining it.
  • Example: college GPA may depend on SAT score, credit hours, and study time together, not any one of those alone.
  • Leaving out a predictor that truly matters doesn’t just lose explanatory power — Residual and Sum-of-Squares Proofs’s omitted-term mechanics (and HW 2 - Regression Inference and Diagnostics Problem 2) show it leaves a systematic, non-random pattern in the residuals.
• The general multiple linear regression (MLR) model adds $k$ predictors instead of one: $Y_i={\beta }_0+{\beta }_1{X}_{i1}+\cdots +{\beta }_k{X}_{ik}+{ϵ}_i$.
  • $E(Y_i\mid X_{i1},\dots ,X_{ik})={\beta }_0+{\sum }_{j=1}^k{\beta }_j{X}_{ij}$ is the population response plane (a hyperplane once $k>2$).
  • The ${\beta }_j$‘s are unknown population parameters; ${\hat{\beta }}_j$ denotes their least-squares estimate from the sample, exactly as in SLR.
  • The same four error assumptions from Fitting a Simple Linear Regression carry over unchanged: $E({ϵ}_i\mid X)=0$, constant variance $\mathrm{Var}({ϵ}_i\mid X)={\sigma }^2$, independence across $i$, and (for exact small-sample inference) normality, ${ϵ}_i\stackrel{ind}{\sim }N(0,{\sigma }^2)$.
• “Linear” describes the coefficients, not the geometric shape in $X$ — this is why polynomial and interaction terms still count as MLR.
  • $E(Y\mid X)={\beta }_0+{\beta }_{1}X+{\beta }_2{X}^2$ is a linear regression model, because it’s linear in ${\beta }_0,{\beta }_1,{\beta }_2$ even though the curve itself is a parabola in $X$.
  • In R, wrap a derived term like $X^2$ in I() — lm(y ~ x + I(x^2)) — so it’s treated as its own new predictor column rather than literal exponentiation inside the formula syntax.
  • With two raw predictors $X_1,X_2$, a complete second-order model can include $X_1,X_2,X_1^2,X_2^2,X_1{X}_2$; the cross term $X_1{X}_2$ is called an interaction and is second order even though neither variable alone is squared.
• Matrix form packages the whole model as one equation, which is what makes a closed-form solution for $k$ predictors possible at all.
  • Stack the model as $Y=X\beta +ϵ$, where $Y$ is $n\times 1$, $X$ is the $n\times (k+1)$ design matrix (a leading column of $1$‘s for the intercept, then one column per predictor), $\beta$ is $(k+1)\times 1$, and $ϵ$ is $n\times 1$.
  • Minimizing $RSS={\sum }_i(Y_i-b_0-b_1{X}_{i1}-\cdots -b_k{X}_{ik}{)}^2=(Y-Xb{)}^{\prime }(Y-Xb)$ over $b$ and setting the matrix derivative to zero gives the normal equations $X^{\prime }Xb=X^{\prime }Y$.
  • As long as $X$ has full column rank (no predictor is an exact linear combination of the others), $X^{\prime }X$ is invertible and $\boxed{\hat{\beta }=(X^{\prime }X{)}^{-1}{X}^{\prime }Y}$ — the direct generalization of ${\hat{\beta }}_1=S_{xy}/S_{xx}$.
  • Fitted values are ${\hat{Y}}_i={\hat{\beta }}_0+{\sum }_j{\hat{\beta }}_j{X}_{ij}$ and residuals are ${\hat{e}}_i=Y_i-{\hat{Y}}_i$, same definitions as SLR, just with more terms inside ${\hat{Y}}_i$.
• The hat matrix writes the whole fitting step as one linear map of $Y$.
  • Define $H=X(X^{\prime }X{)}^{-1}{X}^{\prime }$; substituting $\hat{\beta }=(X^{\prime }X{)}^{-1}{X}^{\prime }Y$ into $\hat{Y}=X\hat{\beta }$ gives $\boxed{\hat{Y}=HY}$ — $H$ sends the observed responses directly to the fitted values.
  • The residual vector is then $\hat{e}=Y-\hat{Y}=Y-HY=(I-H)Y$ — also a linear map of $Y$, via $I-H$ instead of $H$.
  • $H$ depends only on $X$ (the predictor values), never on $Y$ — the same fact that made SLR’s leverage $h_{ii}$ depend only on the $x$‘s.
• $H$‘s algebraic properties are what make variances and expectations of $\hat{e}$ collapse to clean closed forms, so they’re worth having on hand, not just as trivia.
  • $H$ is symmetric, $H^{\prime }=H$: since $(ABC{)}^{\prime }=C^{\prime }{B}^{\prime }{A}^{\prime }$, $H^{\prime }=X[(X^{\prime }X{)}^{-1}{]}^{\prime }{X}^{\prime }$, and $X^{\prime }X$ is itself symmetric (so its inverse is too), giving back $H^{\prime }=X(X^{\prime }X{)}^{-1}{X}^{\prime }=H$.
  • $H$ is idempotent, $H^2=H$: $H^2=X(X^{\prime }X{)}^{-1}{X}^{\prime }X(X^{\prime }X{)}^{-1}{X}^{\prime }$, and the middle $X^{\prime }X\cdot (X^{\prime }X{)}^{-1}=I$ collapses the expression straight back to $H$ — applying the projection twice does nothing new.
  • $I-H$ is also symmetric and idempotent, by the same two facts: $(I-H{)}^{\prime }=I-H$ directly, and $(I-H{)}^2=I-2H+H^2=I-2H+H=I-H$.
  • $HX=X$: the design matrix’s own columns already sit in the space $H$ projects onto, so projecting them changes nothing; consequently $(I-H)X=0$.
  • $\mathrm{tr}(H)=p$, where $p=k+1$ is the number of estimated parameters (including the intercept) — the trace of the hat matrix always equals exactly the number of parameters fit, regardless of the data.
  • Every diagonal entry satisfies $0\le h_{ii}\le 1$.
• Why the fitted plane always passes through the centroid $({\bar{X}}_1,\dots ,{\bar{X}}_k,\bar{Y})$ — a direct consequence of fitting with an intercept, not a coincidence.
  • Minimizing $RSS$ with respect to $b_0$ gives the normal equation ${\sum }_i(Y_i-{\hat{Y}}_i)=0$.
  • Substituting ${\hat{Y}}_i={\hat{\beta }}_0+{\sum }_j{\hat{\beta }}_j{X}_{ij}$ and dividing by $n$: ${\hat{\beta }}_0=\bar{Y}-{\sum }_j{\hat{\beta }}_j{\bar{X}}_j$.
  • Plugging $X_1={\bar{X}}_1,\dots ,X_k={\bar{X}}_k$ into the fitted equation then gives exactly $\hat{Y}({\bar{X}}_1,\dots ,{\bar{X}}_k)=\bar{Y}$ — the plane contains the point of means.
• The variance of the residual vector has a clean matrix form, but it is not simply a diagonal matrix of per-observation variances.
  • Using $\mathrm{Var}(AY)=A\mathrm{Var}(Y)A^{\prime }$ for a fixed matrix $A$, with $A=I-H$ and $\mathrm{Var}(Y)={\sigma }^{2}I$: $\mathrm{Var}(\hat{e})=(I-H){\sigma }^{2}I(I-H{)}^{\prime }={\sigma }^2(I-H)(I-H)={\sigma }^2(I-H)$, using symmetry and idempotence of $I-H$.
  • Written out, the diagonal entries are ${\sigma }^2(1-h_{ii})$ — the familiar per-observation residual variance — but the off-diagonal entries are $-{\sigma }^2{h}_{ij}$, nonzero covariances: fitted residuals are correlated with each other even though the original $Y_i$‘s are independent.
  • $E(\hat{e})=(I-H)E(Y)=(I-H)X\beta =X\beta -HX\beta =X\beta -X\beta =0$, using $HX=X$ — residuals are unbiased for the true errors on average, same conclusion as SLR.
• Leverage generalizes directly: $h_{ii}$ is the $i$-th diagonal entry of $H$, and it only ever depends on the predictor values.
  • In simple linear regression this collapses to the familiar closed form $h_{ii}=\frac{1}{n}+\frac{(x_i-\bar{x}{)}^2}{S_{xx}}$, from Regression Diagnostics and Leverage — a special case of the general $H=X(X^{\prime }X{)}^{-1}{X}^{\prime }$ formula with only one predictor.
  • With $k$ predictors there’s no single clean algebraic formula like SLR’s — $h_{ii}$ has to come from the matrix computation (or hatvalues() in R) — but ${\sum }_i{h}_{ii}=\mathrm{tr}(H)=p$ still always holds, generalizing SLR’s “leverages sum to $2$.”
  • The high-leverage flag generalizes the same way: flag observation $i$ when $h_{ii}$ exceeds twice the average leverage, $\boxed{h_{ii}>\frac{2p}{n}}$ — SLR’s $4/n$ cutoff is exactly this with $p=2$.
• Polynomial regression is just MLR with derived predictor columns, not a separate technique.
  • $Y_i={\beta }_0+{\beta }_1{X}_i+{\beta }_2{X}_i^2+{ϵ}_i$ fits the matrix framework above directly, with the design matrix’s second column being $X_i$ and third column being $X_i^2$ — still linear in ${\beta }_0,{\beta }_1,{\beta }_2$, so $\hat{\beta }=(X^{\prime }X{)}^{-1}{X}^{\prime }Y$ applies unchanged.
  • Because it’s an ordinary MLR fit, every tool above — the hat matrix, its leverage, the $2p/n$ cutoff — applies to a polynomial fit exactly as written, just with $p$ counting however many polynomial terms (plus intercept) were fit.
• If the model isn’t linear in the parameters (e.g. $Y={\beta }_0{e}^{{\beta }_{1}X}+ϵ$), there’s no closed-form $\hat{\beta }=(X^{\prime }X{)}^{-1}{X}^{\prime }Y$ at all — fitting falls back to an iterative numerical search for where $\partial (RSS)/\partial \beta =0$, done by computer since there’s no clean algebraic solution.

Practice Problems

• Problem 1 — polynomial regression as MLR (HW2’s profsalary.txt, $n=143$, fitting $Salary={\beta }_0+{\beta }_{1}Experience+{\beta }_{2}Experienc{e}^2+ϵ$). Fit the quadratic model, verify $\mathrm{tr}(H)=p$, find the leverage cutoff, and compare to the plain linear fit from HW 2 - Regression Inference and Diagnostics Problem 1.

  • Solution.
    • m2 <- lm(Salary ~ Experience + I(Experience^2), data = prof) gives $\widehat{Salary}=34.720+2.872Experience-0.0533Experienc{e}^2$
    • $p=3$ parameters (${\beta }_0,{\beta }_1,{\beta }_2$); sum(hatvalues(m2)) $=3.000$, confirming $\mathrm{tr}(H)=p$
    • cutoff $2p/n=6/143\approx 0.0420$; 18 of the 143 observations exceed it — far more flagged than the linear fit’s handful, since adding a squared term changes which experience values count as “extreme”
    • $R^2$ jumps from $0.675$ (linear, HW2 Problem 1) to $0.925$ (quadratic) — the curvature HW2 Problem 1(g)/(h) found by eye (an upside-down-U residual pattern) is exactly what the squared term captures
    • adjusted $R^2=0.924$, barely below plain $R^2=0.925$ — the one added parameter cost almost nothing against the huge gain in explained variation, unlike adding a genuinely useless predictor (see MLR Inference, ANOVA, and Adjusted R-Squared Problem 2)

• Problem 2 — hat matrix trace and leverage cutoff on a real 4-predictor model (mtcars, $Y=mpg$, $X=$disp, hp, drat, wt, $n=32$). Confirm $\mathrm{tr}(H)=p$, compute the leverage cutoff, and identify which car is flagged.

  • Solution.
    • model1 <- lm(mpg ~ disp + hp + drat + wt, data = mtcars); $p=5$ (intercept plus 4 predictors)
    • sum(hatvalues(model1)) $=5.000$, confirming $\mathrm{tr}(H)=p$
    • cutoff $2p/n=10/32=0.3125$
    • only the Maserati Bora exceeds it, $h_{ii}\approx 0.508$ — its combination of disp, hp, drat, wt sits far enough from the other 31 cars’ predictor values to pull extra weight in the fit, regardless of its actual mpg

• Problem 3 — derive the hat matrix’s symmetry and idempotence from scratch, cold. Show $H^{\prime }=H$ and $H^2=H$ for $H=X(X^{\prime }X{)}^{-1}{X}^{\prime }$.

  • Solution.
    • $H^{\prime }=H$:
      • $H^{\prime }=[X(X^{\prime }X{)}^{-1}{X}^{\prime }{]}^{\prime }=(X^{\prime }{)}^{\prime }[(X^{\prime }X{)}^{-1}{]}^{\prime }{X}^{\prime }$, reversing the order of the product’s transpose
      • $(X^{\prime }{)}^{\prime }=X$
      • $X^{\prime }X$ is symmetric, since $(X^{\prime }X{)}^{\prime }=X^{\prime }(X^{\prime }{)}^{\prime }=X^{\prime }X$, so its inverse is symmetric too: $[(X^{\prime }X{)}^{-1}{]}^{\prime }=(X^{\prime }X{)}^{-1}$
      • substituting both back in: $H^{\prime }=X(X^{\prime }X{)}^{-1}{X}^{\prime }=H$
    • $H^2=H$:
      • $H^2=X(X^{\prime }X{)}^{-1}{X}^{\prime }\cdot X(X^{\prime }X{)}^{-1}{X}^{\prime }$
      • the middle factors multiply to $X^{\prime }X$: $H^2=X(X^{\prime }X{)}^{-1}[X^{\prime }X](X^{\prime }X{)}^{-1}{X}^{\prime }$
      • $(X^{\prime }X{)}^{-1}(X^{\prime }X)=I$, collapsing the bracketed piece
      • $H^2=X\cdot I\cdot (X^{\prime }X{)}^{-1}{X}^{\prime }=X(X^{\prime }X{)}^{-1}{X}^{\prime }=H$