2 Linear regression

The full rank Gauss-Markov model

\[\underset{n\times 1}{\b{y}} = \underset{n\times p}{\b{X}}\quad \underset{p\times 1}{\bs{\beta}}+\underset{n\times 1}{\b{e}}\]

with \(\b{e}\sim N_n(\b{0}, \sigma^2\b{I})\) and \(r(\b{X})=p\)

From hereafter:

• \(y\) may represent any transformation of the original response variable

• \(\b{X}\) may be the design matrix for any transformation of the covariates \((\mathbb{R}^k\rightarrow \mathbb{R}^{p-1})\)

Usually, the entry point to regression models is the 1^st order model that simply includes all the covariates in the design matrix. If required, the inclusion of higher order terms can be considered:

Note

The usual interpretation of \(\beta_i\), \(i=1,\ldots, p-1\), as the variation in \(E[y\mid \b{x}]\) due to an unit increase in \(x_i\), keeping all other covariates fixed, that is \[E[y\mid x_1,\ldots, x_i=x+1,\ldots, x_{p-1}]-E[y\mid x_1,\ldots, x_i=x,\ldots, x_{p-1}]\] is only valid for the 1^st order model.

The matrix \(\b{X}'\b{X}\) is:

1. symmetric

2. full rank
3. invertible
4. positive definite

Therefore \[\bs{\beta}=(\b{X}'\b{X})^{-1}\b{X}'\b{y}\] is the unique solution.

Properties of the estimator

\[\hat{\bs{\beta}}=(\b{X}'\b{X})^{-1}\b{X}'\b{y}\]

1. Linear estimator

2. Unbiased estimator of \(\bs{\beta}\)

3. \(Var[\hat{\bs{\beta}}] = \sigma^2(\b{X}'\b{X})^{-1}\)

ML estimation

Maximizing \[\log \mathcal{L}(\bs{\beta},\sigma^2)=-\frac{n}{2}\log \sigma^2-\frac{1}{2\sigma^2}(\b{y}-\b{X}\bs{\beta})'(\b{y}-\b{X}\bs{\beta}) + C\] we get:

1. the same estimator for \(\bs{\beta}\)

2. \(\hat{\sigma}^2_{ML}=\frac{(\b{y}-\b{X}\hat{\bs{\beta}})'(\b{y}-\b{X}\hat{\bs{\beta}})}{n}=\frac{SSE}{n}\)

Residuals

\(\hat{\b{y}}=\hat{E}[\b{y}\mid \b{x}]=\b{X}\hat{\bs{\beta}}\) – the fitted values

A similar thing for the residuals:

\[\b{e}=\b{y}-\hat{\b{y}}=(\b{I}-\b{H})\b{y}\]

Some properties of the residuals

1. \(E[\b{e}]=\b{0}\)

2. \(Var[\b{e}]=\sigma^2(\b{I}-\b{H})\)

3. \(\hat{\b{y}}'\b{e}=0\)

4. \(Cov[\hat{\b{y}},\b{e}] = \b{0}\)

The SSE in matrix form:

\[SSE=\b{e}'\b{e}=\b{y}'(\b{I}-\b{H})\b{y}\]

Therefore, \(E[SSE]=\sigma^2(n-p)\) which shows that \(E[\hat{\sigma}^2_{ML}]=\frac{n-p}{n}\sigma^2\).

Distributions of the estimators

1. \(\hat{\bs{\beta}} \sim N_p\left(\bs{\beta}, \sigma^2(\b{X}'\b{X})^{-1}\right)\)

2. \(\frac{(\hat{\bs{\beta}}-\bs{\beta})'\b{X}'\b{X}(\hat{\bs{\beta}}-\bs{\beta})}{\sigma^2}\sim \chi^2_{(p)}\)

3. \(\frac{SSE}{\sigma^2} = \frac{(n-p)\hat{\sigma}^2}{\sigma^2}\sim \chi^2_{(n-p)}\)

4. \(\hat{\bs{\beta}}\) and \(\hat{\sigma}^2\) are independent

The previous distributions provide a few more useful results.

\(\forall \b{c}\in\mathbb{R}^p :\)

\[\b{c}'\hat{\bs{\beta}}\sim N\left(\b{c}'\bs{\beta}, \sigma^2\b{c}'(\b{X}'\b{X})^{-1}\b{c}\right)\]

Particular case 1

For \(\b{c}' = (0,\ldots, 0, 1, 0, \ldots, 0)\):

\[\frac{\hat{\beta}_j-\beta_j}{se(\hat{\beta}_j)}\sim t_{(n-p)}\]

with \(se(\hat{\beta}_j)=\hat{\sigma}\sqrt{(\b{X}'\b{X})^{-1}_{jj}}\)

Particular case 2

For \(y_0=E[y\mid \b{x}=\b{x}_0] = \b{x}'_0\bs{\beta}\):

\[\frac{\b{x}'_0\hat{\bs{\beta}}-\b{x}'_0\bs{\beta}}{\hat{\sigma}\sqrt{\b{x}'_0 (\b{X}'\b{X})^{-1} \b{x}_0}}\sim t_{(n-p)}\]

Prediction

The goal is to predict a new observation under the model, that is

\[y_u\sim N(\b{x}'_u\bs{\beta}, \sigma^2)\]

independent from observed data.

Considering the prediction error, \(\hat{y}_u-y_u\), we have:

1. \(E[\hat{y}_u-y_u]=0\)

2. \(Var[\hat{y}_u-y_u]=\sigma^2\left(1+\b{x}'_u(\b{X}'\b{X})^{-1}\b{x}_u\right)\)

3. \(\dfrac{\hat{y}_u-y_u}{se(\hat{y}_u-y_u)}\sim t_{(n-p)}\)

   where \(se(\hat{y}_u-y_u)=\hat{\sigma}\sqrt{1+\b{x}'_u(\b{X}'\b{X})^{-1}\b{x}_u}\)

Another result

From (2) to (4) we have:

\[\frac{n-p}{p}\frac{(\hat{\bs{\beta}}-\bs{\beta})'\b{X}'\b{X}(\hat{\bs{\beta}}-\bs{\beta})}{SSE}\sim F_{(p, n-p)}\]

How to interpret the following results?

ci <- confint(model, level = 0.96)
ci

               2 %  98 %
(Intercept) -81.05 138.7
education     2.69  23.3

These results where obtained from

\[\frac{\hat{\bs{\beta}}_j-\bs{\beta}_j}{se(\hat{\bs{\beta}}_j)}\sim t_{(n-p)}\] with \(se(\hat{\bs{\beta}}_j)=\hat{\sigma}\sqrt{(\b{X}'\b{X})^{-1}_{jj}}\).

The problem of simultaneous estimation

Using the distributions of the estimators, we have

\[T(\b{y})=\frac{n-p}{p}\frac{(\hat{\bs{\beta}}-\bs{\beta})'\b{X}'\b{X}(\hat{\bs{\beta}}-\bs{\beta})}{SSE}\sim F_{(p, n-p)}\]

A possible correction

Let \(C_i(\b{y})\) be a confidence interval for \(\beta_i\) with a \(\gamma_i\) confidence level.

Then, \[\begin{align*} P\left(\bs{\beta}\in \bigcap_{i=0}^{p-1} C_i(\b{y})\right) & \geq 1-\sum_{i=0}^{p-1} P(\beta_i\not\in C_i(\b{y}) )=\\ & = 1-\sum_{i=0}^{p-1} (1-\gamma_i)=\\ & = 1-p(1-\gamma), \text{ if } \gamma_i=\gamma, \forall i \end{align*}\]

Application

Inferences for the mean response

From

\[\frac{\b{x}'_0\hat{\bs{\beta}}-\b{x}'_0\bs{\beta}}{se(\b{x}'_0\hat{\bs{\beta}})}\sim t_{(n-p)}\]

a \(100\gamma\%\) confidence interval is given by

\[\b{x}'_0\hat{\bs{\beta}}\pm F_{t_{(n-p)}}^{-1}\left(\frac{1+\gamma}{2}\right)\times se(\b{x}'_0\hat{\bs{\beta}})\]

96% confidence bands for the mean response

As before, corrections are also required when considering more than one prediction interval:

1. the Bonferroni method

2. the Scheffé method for \(k\) intervals

   \[\hat{y}_u\pm S \times se(\hat{y}_u-y_u)\]

   with \(S^2 = k F_{F_{(k, n-p)}}^{-1}\left(\gamma\right)\)

96% prediction bands

R package: api2lm

Notes

1. \(H_0\) means that the set of terms is not capable of explaining \(y\) under the model or

   \[H_0: E[y\mid \b{X}] = \beta_0\]

The test procedure is based on the following decomposition of deviation:

\[y_i-\bar{y} = (y_i-\hat{y}_i) + (\hat{y}_i-\bar{y})\]

We have already seen that

\[SSE=\b{y}'(\b{I}-\b{H})\b{y}\]

In matrix notation, it is clear that:

1. \(SST = SSE + SSR\)

And so, to test the hypotheses we can use the overall F-test with the statistic:

\[F=\frac{\frac{SSR}{p-1}}{\frac{SSE}{n-p}} = \frac{MSR}{MSE}\stackrel{H_0}{\sim}F_{(p-1, n-p)}\]

Related to the previous test we already know

\[R^2=\frac{SSR}{SST}=1-\frac{SSE}{SST}\] that is now called the multiple coefficient of determination.

Because of this \(R^2\) inflation, it should be preferred the

The general linear hypotheses

Many hypotheses on the regression coefficients can be written in the generic form

\[H_0: \b{C}\bs{\beta}= \b{m}\]

versus

\[H_1: \b{C}\bs{\beta}\neq \b{m}\]

Under \(H_0\):

\[\b{C}\hat{\bs{\beta}} - \b{m} \sim N_c\left(\b{0}, \sigma^2\b{C}(\b{X}'\b{X})^{-1}\b{C}'\right)\]

This framework provides a way to compare a larger model (L) with a reduced one (R) defined by the restritions.

For nested models:

\[Q = SSE(R) - SSE(L)\]

and the F-statistic can be written in a more practical way:

\[F = \frac{SSE(R) - SSE(L)}{SSE(L)}\frac{df_L}{df_R-df_L} \stackrel{H_0}{\sim}F_{(df_R-df_F, df_F)}\]

Note

Application of this general framework to the initial examples of hypotheses would produce exactly the previous tests. For the univariate t-tests we would have \(F_o =T_o^2\).

The tests to compare nested models are usually called partial F-tests.

Extra sums of squares

Each regression model \(M\) leads to a particular decomposition of the SST: \[SST = SSR(M) + SSE(M)\]

To do that we define partial or extra sums of squares as

\[SSR(L) - SSR(R) = SSE(R) - SSE(L)\]

Another way of using the extra sums of squares is

\(SSR(x_1,\ldots,x_{p-1})=\)

\(\quad\quad = SSR(x_1,\ldots,x_{j}) + SSR(x_{j+1}, \ldots, x_{p-1}\mid x_1,\ldots,x_{j})\)

The extra sums of squares are also useful to define the coefficients of partial determination.

For example,

\[R^2_{3\mid 1,2}=\frac{SSR(x_3\mid x_1, x_2)}{SSE(x_1, x_2)}\]

is the coefficient of partial determination between \(y\) and \(x_3\), given that \(x_1\) and \(x_2\) are in the model, and it is interpreted as the proportion of variation not explained by the reduced model that can be explained by the larger model.

Finally, we can consider a sequence of partial F-tests.

anova(model)

Analysis of Variance Table

Response: demand
             Df Sum Sq Mean Sq F value Pr(>F)   
education     1    754     754   19.06 0.0072 **
income        1    325     325    8.21 0.0352 * 
urbanization  1      2       2    0.06 0.8221   
Residuals     5    198      40                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Another view of the previous ANOVA table with added information:

m0 <- lm(demand ~ 1, data= product)
m1 <- update(m0, . ~ . + education)
m2 <- update(m1, . ~ . + income)
m3 <- update(m2, . ~ . + urbanization)
anova(m0, m1, m2, m3)

Analysis of Variance Table

Model 1: demand ~ 1
Model 2: demand ~ education
Model 3: demand ~ education + income
Model 4: demand ~ education + income + urbanization
  Res.Df  RSS Df Sum of Sq     F Pr(>F)   
1      8 1279                             
2      7  525  1       754 19.06 0.0072 **
3      6  200  1       325  8.21 0.0352 * 
4      5  198  1         2  0.06 0.8221   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Binary covariates (2 levels)

Suppose that \(p=3\) and \(x_2\) is binary (0/1). The response surface for the first-order model is:

\[\begin{align*} E[Y\mid \b{x}] & = \beta_0 + \beta_1 x_1 + \beta_2 x_2 =\\[0.25cm] & = \begin{cases}\beta_0+\beta_1 x_1, & x_2=0\\[0.25cm](\beta_0 + \beta_2)+\beta_1 x_1, & x_2=1\end{cases} \end{align*}\]

Note

• With \(q\) binary covariates we would have \(2^q\) parallel response surfaces.

• \(E[y\mid x_1, x_2=1] - E[y\mid x_1, x_2=0] = \beta_{2}\)

To allow non-parallel surfaces we need to include interactions:

\[\begin{align*} E[Y\mid \b{x}] & = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_1 x_2 =\\[0.25cm] & = \begin{cases}\beta_0+\beta_1 x_1, & x_2=0\\[0.25cm](\beta_0 + \beta_2)+(\beta_1 + \beta_3) x_1, & x_2=1\end{cases} \end{align*}\]

More than 2 levels

Suppose now that \(x_2\) is a categorical variable or factor with 3 levels. If we code the levels as 0, 1 and 2 we have:

\[\begin{align*} E[Y\mid \b{x}] & = \beta_0 + \beta_1 x_1 + \beta_2 x_2 =\\[0.25cm] & = \begin{cases}\beta_0+\beta_1 x_1, & x_2=0\\[0.25cm](\beta_0 + \beta_2)+\beta_1 x_1, & x_2=1\\[0.25cm](\beta_0 + 2\beta_2)+\beta_1 x_1, & x_2=2\end{cases} \end{align*}\]

We could associate a binary variable to each level of the factor \(x_2\rightsquigarrow (x_{2,0}, x_{2,1}, x_{2,2})\) but this would lead to an incomplete rank linear model.

\[\begin{align*} E[Y\mid \b{x}] & = \beta_0 + \beta_1 x_1 + \beta_{2,1} x_{2,1} + \beta_{2,2} x_{2,2} =\\[0.25cm] & = \begin{cases}\beta_0+\beta_1 x_1, & x_2=0\\[0.25cm](\beta_0 + \beta_{2,1})+\beta_1 x_1, & x_2=1\\[0.25cm](\beta_0 + \beta_{2,2})+\beta_1 x_1, & x_2=2\end{cases} \end{align*}\]

Again, the inclusion of interactions will remove the parallelism.

\[\begin{align*} E[Y\mid \b{x}] & = \beta_0 + \beta_1 x_1 + \beta_{2,1} x_{2,1} + \beta_{2,2} x_{2,2} + \\[0.25cm] & + \beta_{3} x_1 x_{2,1} + \beta_4 x_1 x_{2,2} \\[0.25cm] & = \begin{cases}\beta_0+\beta_1 x_1, & x_2=0\\[0.25cm](\beta_0 + \beta_{2,1})+(\beta_1 + \beta_3) x_1, & x_2=1\\[0.25cm](\beta_0 + \beta_{2,2}) + (\beta_1 + \beta_4) x_1, & x_2=2\end{cases} \end{align*}\]

Notes

1. Instead of the binary encoding in the dummy variables, other values (contrasts) are also used that provide different interpretations for the regression parameters.

Another application of categorical variables

Piecewise linear regression

Suppose \(x_1\) be a quantitative covariate and let \(x_2=I_{[c, +\infty[}(x_1)\).

2. Consider

   \[\begin{align*} E[Y\mid \b{x}] & = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 (x_1-c)x_2 =\\[0.25cm] & = \begin{cases}\beta_0+\beta_1 x_1, & x_1<c\\[0.25cm](\beta_0 + \beta_2 - c \beta_3)+(\beta_1 + \beta_3) x_1, & x_1\geq c\end{cases} \end{align*}\]

   \(\leadsto\) discontinuous piecewise regression model