---
title: Inference in multiple linear regression
author: STATS 191
date: 'April 1, 2024'
engine: knitr
format:
  revealjs:
    html-math-method: mathjax
webr:
  show-startup-message: false
filters:
  - webr
jupyter:
  jupytext:
    formats: 'qmd:quarto,ipynb'
    text_representation:
      extension: .qmd
      format_name: quarto
      format_version: '1.0'
      jupytext_version: 1.16.1
  kernelspec:
    display_name: R
    language: R
    name: ir
---


## Outline

- Case studies:

    A. Galileo's falling bodies

    B.  Energy costs of echolocation

## Energy costs of echolocation

```{webr-r}
url = 'https://raw.githubusercontent.com/StanfordStatistics/stats191-data/main/Sleuth3/bats_energy.csv'
bats = read.table(url, sep=',', header=TRUE)
bats
```

## Fitting a model

```{webr-r}
bats.lm = lm(Energy ~ Type + Mass, data=bats)
summary(bats.lm)
```

## Other models

```{webr-r}
simpler.lm = lm(Energy ~ Type, data=bats)
coef(simpler.lm)
```

```{webr-r}
complex.lm = lm(Energy ~ Type + Mass + Type:Mass, data=bats)
coef(complex.lm)
```

## Case study B: Galileo's data

```{webr-r}
url = 'https://raw.githubusercontent.com/StanfordStatistics/stats191-data/main/Sleuth3/galileo.csv'
galileo = read.table(url, sep=',', header=TRUE)
with(galileo, plot(Height, Distance))
```

---

- Galileo fit a quadratic model to his data

- Note the notation `I(Height^2)` -- without `I` a quadratic term will not be added...

```{webr-r}
galileo1.lm = lm(Distance ~ Height + I(Height^2), data=galileo)
summary(galileo1.lm)
```

---

- A different way to fit the model

```{webr-r}
galileo2.lm = lm(Distance ~ poly(Height, 2), data=galileo)
summary(galileo2.lm)
```

---

- Predictions / CIs are the same

```{webr-r}
predict(galileo1.lm, list(Height=250), interval='confidence')
predict(galileo2.lm, list(Height=250), interval='confidence')
```

---

- Tests of quadratic effect are the same

```{webr-r}
galileo0.lm = lm(Distance ~ Height, data=galileo)
anova(galileo0.lm, galileo1.lm)
anova(galileo0.lm, galileo2.lm)
```

## Confidence intervals 

-   Suppose we want a $(1-\alpha)\cdot 100\%$ CI for
    $\sum_{j=0}^p a_j\beta_j$.

-   Just as in simple linear regression:

    $$\sum_{j=0}^p a_j \widehat{\beta}_j \pm t_{1-\alpha/2, n-p-1} \cdot SE\left(\sum_{j=0}^p a_j\widehat{\beta}_j\right).$$


```{webr-r}
confint(bats.lm, level=0.90)
predict(bats.lm, list(Mass=300, Type="non-echolocating bats"), interval='confidence')
```




## $T$-statistics revisited

Of course, these confidence intervals are based on the standard ingredients of a
$T$-statistic.

-   Suppose we want to test

\begin{equation}
H_0:\sum_{j=0}^p a_j\beta_j= h.
\end{equation}

---

- As in simple linear regression, it is based on

\begin{equation}
T = \frac{\sum_{j=0}^p a_j \widehat{\beta}_j - h}{SE(\sum_{j=0}^p a_j \widehat{\beta}_j)}.
\end{equation}

-   If $H_0$ is true, then $T \sim t_{n-p-1}$, so we reject $H_0$ at
    level $\alpha$ if

\begin{equation}
\begin{aligned}
|T| &\geq t_{1-\alpha/2,n-p-1}, \qquad \text{ OR} \\
p-\text{value} &= {\tt 2*(1-pt(|T|, n-p-1))} \leq \alpha.
\end{aligned}
\end{equation}


---

Let's do a quick calculation to remind ourselves the relationships of the variables
in the table above.


```{webr-r}
T_Mass = 0.057495 / 0.007557
P_Mass = 2 * (1 - pt(abs(T_Mass), 16))
c(T_Mass, P_Mass)
```

```{webr-r}
summary(bats.lm)$coef
```

## One-sided tests

-   Suppose, instead, we wanted to test the one-sided hypothesis

\begin{equation}
H_0:\sum_{j=0}^p a_j\beta_j \leq  h, \  \text{vs.} \ H_a: \sum_{j=0}^p a_j\beta_j > 
 h
\end{equation}

-   If $H_0$ is true, then $T$ is no longer exactly $t_{n-p-1}$ but we still have

\begin{equation}
\mathbb{P}\left(T > t_{1-\alpha, n-p-1}\right) \leq 1 - \alpha
\end{equation}

---

- We reject $H_0$ at level $\alpha$ if

\begin{equation}
\begin{aligned}
T &\geq t_{1-\alpha,n-p-1}, \qquad \text{ OR} \\
p-\text{value} &= {\tt (1-pt(T, n-p-1))} \leq \alpha.
\end{aligned}
\end{equation}
       
## Standard error of $\sum_{j=0}^p a_j \widehat{\beta}_j$

- In order to form these $T$ statistics, we need the $SE$ of our estimate $\sum_{j=0}^p a_j \widehat{\beta}_j$.

-   Based on matrix approach to regression

\begin{equation}
SE\left(\sum_{j=0}^p a_j\widehat{\beta}_j \right) = SE\left(a^T\widehat{\beta} \right) =
\sqrt{\widehat{\sigma}^2 \cdot a^T (X^TX)^{-1} a}.
\end{equation}

---

-   Don’t worry too much about specific implementation – for much of the effects
    we want `R` will do this for you in
    general.



```{webr-r}
X = model.matrix(bats.lm)
sigma.hat = sqrt(sum(resid(bats.lm)^2 / 16))
cov.beta = sigma.hat^2 * solve(t(X) %*% X)
cov.beta
```

```{webr-r}
vcov(bats.lm)
```

---

The standard errors of each coefficient estimate are the square root of the diagonal entries. They appear as the
`Std. Error` column in the `coef` table.


```{webr-r}
sqrt(diag(vcov(bats.lm)))
summary(bats.lm)$coef[,2]
```

## Prediction / forecasting interval

-   Basically identical to simple linear regression.

-   Prediction interval at $X_{1,new}, \dots, X_{p,new}$:

\begin{equation}
\begin{aligned}
\widehat{\beta}_0 + \sum_{j=1}^p X_{j,new} \widehat{\beta}_j\pm t_{1-\alpha/2, n-p-1} \times & \\
\qquad \sqrt{\widehat{\sigma}^2 + SE\left(\widehat{\beta}_0 + \sum_{j=1}^p X_{j,new}\widehat{\beta}_j\right)^2} &
\end{aligned}
\end{equation}

## Questions about many (combinations) of $\beta_j$’s

-   In multiple regression we can ask more complicated questions than in
    simple regression.

-   For instance, in `bats.lm` we could ask whether `Type` is important at all?

-   These questions can be answered answered by $F$-statistics.


## Dropping one or more variables

-   Suppose we wanted to test the above hypothesis
    Formally, the null hypothesis is:

$$
H_0: \beta_1=\beta_2=0
$$

-    The alternative is

$$
H_a = \text{one of $ \beta_1,\beta_2$ is not 0}. 
$$

---

-   This test is an $F$-test based on two models

\begin{equation}
\begin{aligned}
\text{Full:}  & \qquad
\texttt{Energy $\widetilde{}$ Type + Mass} \\
\text{Reduced:} & \qquad \texttt{Energy $\widetilde{}$ Mass} \\
\end{aligned}
\end{equation}


## $SSE$ of a model

-   In the graphic, a "model", ${\cal M}$ is a subspace of
    $\mathbb{R}^n$ (e.g. column space of ${X}$).

-   Least squares fit = projection onto the subspace of ${\cal M}$,
    yielding predicted values $\widehat{Y}_{{\cal M}}$

-   Error sum of squares:

$$SSE({\cal M}) = \|Y - \widehat{Y}_{{\cal M}}\|^2.$$

## Least squares for $F$ statistic

![](http://www.stanford.edu/class/stats191/Chapter9/figs/axes_SSR_multiple.png){height=400 fig-align="center"}

- Fits of a full and reduced model $\hat{Y}_F$ and $\hat{Y}_R$

- The difference $\hat{Y}_F-\hat{Y}_R$.

## Right triangle for $F$ statistic

![](http://www.stanford.edu/class/stats191/Chapter8/figs/right_triangle_reduced.png){height=300 fig-align="center"}

- Sides of the triangle: $SSE_R-SSE_F$, $SSE_F$

- Hypotenuse: $SSE_R$

---

### Right triangle with full and reduced model: degrees of freedom

![](http://www.stanford.edu/class/stats191/Chapter8/figs/right_triangle_reduced_df.png){height=300 fig-align="center"}

- Sides of the triangle: $df_R-df_F$, $df_F$

- Hypotenuse: $df_R$



## $F$-statistic for $H_0:\beta_{1}=\beta_{2}=0$

-   We compute the $F$ statistic the same to compare any two (nested) models

\begin{equation}
\begin{aligned}
F &=\frac{\frac{SSE(R) - SSE(F)}{2}}{\frac{SSE(F)}{n-1-p}} \\
& \sim F_{2, 16}       \qquad   (\text{if $H_0$ is true})
\end{aligned}
\end{equation}

---

-   Reject $H_0$ at level $\alpha$ if $F > F_{1-\alpha, 2, 16}$.


```{webr-r}
reduced.lm = lm(Energy ~ Mass, data=bats)
anova(reduced.lm, bats.lm)
```

## General $F$-tests

-   Given two models $R \subset F$ (i.e. $R$ is a subspace of $F$), we
    can consider testing

\begin{equation}
H_0: \text{$R$ is adequate (i.e. $\mathbb{E}(Y) \in R$)}
\end{equation}

\begin{equation}
H_a: \text{$F$ is adequate (i.e. $\mathbb{E}(Y) \in F$)}
\end{equation}
    
- The least squares picture has models $X_R$ and $X_F=X_R+(X_F \perp X_R) \dots$

---

- The test statistic is

\begin{equation}
F = \frac{(SSE(R) - SSE(F)) / (df_R - df_F)}{SSE(F)/df_F}
\end{equation}

-   If $H_0$ is true, $F \sim F_{df_R-df_F, df_F}$ so we reject $H_0$ at
    level $\alpha$ if $F > F_{1-\alpha, df_R-df_F, df_F}$.


## Constraining coefficients

- Suppose we wanted to test $H_0$: the line for non-echolocating
bats has the same intercept as the line for non-echolocating birds.

- Can be expressed as $H_0:\beta_1=\beta_2$ in `bats.lm`.

---

### Strategy 1: fit a model in which this is *forced* to be true

```{webr-r}
null_group = rep(1, nrow(bats))
null_group[bats$Type=="echolocating bats"] = 2
bats$null_group = factor(null_group)
null_bats.lm = lm(Energy ~ Mass + null_group, data=bats)
anova(null_bats.lm, bats.lm)
```

---

### Strategy 2: a $T$-statistic

- Hypothesis is $H_0:\beta_1-\beta_2=0$

- This method doesn't require fitting the special model `null_bats.lm`!

- Can be generalized to $F$ tests (hypotheses involving multiple contrasts of $\beta$)

```{webr-r}
A = c(0,1,-1,0)
var_diff = sum(A * vcov(bats.lm) %*% A)
se_diff = sqrt(var_diff)
T = sum(coef(bats.lm) * A) / se_diff
T^2
```

## Math aside: general linear hypothesis

- Suppose we want to test the null hypothesis

\begin{equation}
H_0:C_{q \times (p+1)}\beta_{(p+1) \times 1} = h
\end{equation}

- Alternative is

\begin{equation}
H_a :C_{q \times (p+1)}\beta_{(p+1) \times 1} \neq h.
\end{equation}

## Math aside: $F$ statistic in general linear hypothesis

- Numerator

$$
(C\hat{\beta}-h)^T \left(C(X^TX)^{-1}C^T \right)^{-1} (C\hat{\beta}-h) / q
$$

- Denominator: the usual MSE

- We just used special case $q=1$ above...
