```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
library(faraway)
```

## Goals for today

* Understand the concept and practice of partitioning sums-of-squares

> - Understand the uses of R<sup>2</sup> and adjusted-R<sup>2</sup> for linear models


## General approach

### Question $\rightarrow$ Data $\rightarrow$ Model $\rightarrow$ <u>Inference</u> $\rightarrow$ Prediction


## Inference | Relationship(s) between predictor(s) and response

* How big is an effect? Is an effect "important"? 

> - How certain are we about the sign and size of the effect?  


## General approach

### Question $\rightarrow$ Data $\rightarrow$ Model $\rightarrow$ Inference $\rightarrow$ <u>Prediction</u>


## Prediction 

* What do we predict for a novel case? 
      
> - How good will our predictions be? 


## Partitioning variance

In general, we have something like

$$
DATA = MODEL + ERRORS
$$

and hence

$$
\text{Var}(DATA) = \text{Var}(MODEL) + \text{Var}(ERRORS)
$$


## Partitioning total deviations

The total deviations in the data equal the sum of those for the model and errors

$$
\underbrace{y_i - \bar{y}}_{\text{Total}} = \underbrace{\widehat{y}_i - \bar{y}}_{\text{Model}} + \underbrace{y_i - \widehat{y}_i}_{\text{Error}}
$$


## Partitioning variance | An example

Let's consider a model for the relationship between

* soil surface temperature (C) measured across 15 small plots on a hill, and

* the aspect of each plot (direction in degrees; N to S)


## Partitioning variance | An example

```{r ex_SS, echo = FALSE, fig.width = 5, fig.height = 4, fig.align = "center"}
nn <- 15

aspect <- round(runif(nn, 0, 180))

b1 <- 5
b2 <- 0.05

ee <- rnorm(nn, 0, 1.5)

temp <- b1 + b2 * aspect + ee

par(mai = c(0.9,0.9,0.1,0.1), omi = c(0, 0, 0, 0))

## plot the data
plot(aspect, temp, type = "p", pch = 19, las = 1,
     ylim = c(floor(min(temp)) - 1, ceiling(max(temp))),
     bty = "l", cex.axis = 1.2, cex.lab = 1.3,
     xlab = "Aspect (Degrees)",
     ylab = "Temperature (C)")
```


## Partitioning variance | An example

Let's fit this model:

$$temperature_i = \beta_0 + \beta_1 \times aspect_i + e_i$$

and find $\hat{\boldsymbol{\beta}}$ and $\hat{\boldsymbol{y}}$

<br>

<center>_Let's use matrix algebra like we learned last time!_</center>


## Partitioning variance | An example 

Create our design matrix $\mathbf{X}$ 

```{r ex_SS_2, echo=TRUE}
## construct our "Design matrix"
XX <- cbind(intercept = rep(1, nn), x = aspect)
head(XX)
``` 

## Partitioning variance | An example  

We can calculate our betas using

$\hat{\boldsymbol{\beta}} = (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbf{y}$

```{r ex_SS_3, echo = TRUE}
## calculate betas  
beta_hat <- solve(t(XX) %*% XX) %*% t(XX) %*% temp
``` 

<br>

Recall that `t()` transforms (rotates) a matrix and `solve()` inverts it


## Partitioning variance | An example  

We can calculate our model estimates ($\hat{\mathbf{y}}$) using

$\hat{\mathbf{y}} = \mathbf{X}\hat{\boldsymbol{\beta}}$

```{r ex_SS_4, echo = TRUE}
## calculate model estimates  
temp_hat <- XX %*% beta_hat
``` 


## Partitioning variance | An example  

Let's plot our model estimates (ie, the fitted regression line)

```{r ex_S_5, echo = FALSE, fig.width = 5, fig.height = 4, fig.align = "center"}
par(mai = c(0.9,0.9,0.1,0.1), omi = c(0, 0, 0, 0))

## plot the data
plot(aspect, temp, type = "p", pch = 19, las = 1,
     ylim = c(floor(min(temp)) - 1, ceiling(max(temp))),
     bty = "l", cex.axis = 1.2, cex.lab = 1.3,
     xlab = "Aspect (Degrees)",
     ylab = "Temperature (C)")
abline(a = beta_hat[1], b = beta_hat[2], col = "blue", lwd = 2)
```


## Partitioning total deviations

```{r partition_SS, echo = FALSE, fig.width = 8, fig.height = 3.5, fig.align="center"}
## set plot area
par(mfrow = c(1,3),
    mai = c(0.5,0.3,0.5,0.2),
    omi = c(0, 0, 0, 0),
    cex = 0.7)

## SSTO - plot y and mean(y)
plot(aspect, temp, type = "n",
     ylim = c(floor(min(temp)) - 1, ceiling(max(temp))),
     bty = "l", xaxt = "n", yaxt = "n",
     xlab = "",
     ylab = "")
segments(x0 = aspect, y0 = temp, y1 = mean(temp), col = "gray")
abline(h = mean(temp), lty = "dashed")
points(aspect, temp, pch = 19)
mtext(expression(Total: ~ italic(y[i])-bar(italic(y))), side = 3, line = 1)
mtext(expression(italic(x)), 1, line=0.2, cex=1.1, at=max(aspect))
mtext(expression(italic(y)), 2, line=0.3, cex=1.1, at=max(temp), las=1)
mtext(expression(italic(bar(y))), 2, line=0.3, cex=1.1, at=mean(temp), las=1)

## SSR - plot temp and mean(temp)
plot(aspect, temp, type = "n",
     ylim = c(floor(min(temp)) - 1, ceiling(max(temp))),
     bty = "l", xaxt = "n", yaxt = "n",
     xlab = "",
     ylab = "")
segments(x0 = aspect, y0 = temp_hat, y1 = mean(temp_hat), col = "gray")
points(aspect, temp_hat, pch = 21, bg = "white")
abline(a = beta_hat[1], b = beta_hat[2])
abline(h = mean(temp), lty = "dashed")
points(aspect, temp_hat, pch = 21, bg = "white")
mtext(expression(Model: ~ italic(hat(y)[i])-bar(italic(y))), side = 3, line = 1)
mtext(expression(italic(x)), 1, line=0.2, cex=1.1, at=max(aspect))
mtext(expression(italic(y)), 2, line=0.3, cex=1.1, at=max(temp), las=1)
mtext(expression(italic(bar(y))), 2, line=0.3, cex=1.1, at=mean(temp), las=1)
mtext(expression(italic(hat(y))), 4, line=0.2, cex=1.1, at=max(temp), las=1)

## SSE - plot observed and temp
plot(aspect, temp, type = "n",
     ylim = c(floor(min(temp)) - 1, ceiling(max(temp))),
     bty = "l", xaxt = "n", yaxt = "n",
     xlab = "",
     ylab = "")
segments(x0 = aspect, y0 = temp_hat, y1 = temp, col = "gray")
points(aspect, temp, pch = 19)
abline(a = beta_hat[1], b = beta_hat[2])
mtext(expression(Error: ~ italic(y[i])-italic(hat(y)[i])), side = 3, line = 1)
mtext(expression(italic(x)), 1, line=0.3, cex=1.1, at=max(aspect))
mtext(expression(italic(y)), 2, line=0.3, cex=1.1, at=max(temp), las=1)
mtext(expression(italic(hat(y))), 4, line=0.2, cex=1.1, at=max(temp), las=1)
```


## Partitioning sums-of-squares

The sums-of-squares have the same additive property as the deviations

$$
\underbrace{\sum (y_i - \bar{y})^2}_{SSTO} = \underbrace{\sum (\widehat{y}_i - \bar{y})^2}_{SSR} + \underbrace{\sum (y_i - \widehat{y}_i)^2}_{SSE}
$$


## Sum-of-squares: Total

The total sum-of-squares $(SSTO)$ measures the total variation in the data as the differences between the data and their mean

$$
SSTO = \sum \left( y_i - \bar{y} \right)^2
$$


## Sum-of-squares: Model

The model (regression) sum-of-squares $(SSR)$ measures the variation between the model fits and the mean of the data

$$
SSR = \sum \left( \widehat{y}_i - \bar{y} \right)^2
$$


## Sum-of-squares: Error

The error sum-of-squares $(SSE)$ measures the variation between the data and the model fits

$$
SSE = \sum \left( y_i - \widehat{y}_i \right)^2
$$


## Sum-of-squares | An example  

Let's calculate the total sum-of-squares (SSTO) using

$SSTO = \sum \left( y_i - \bar{y} \right)^2$

```{r ex_SS_6, echo = TRUE}
## mean of the response
y_bar <- mean(temp)

## total sum-of-squares
SSTO <- t(temp - y_bar) %*% (temp - y_bar)
``` 

<br>

Recall that $\mathbf{x}^{\top} \mathbf{x}$ will give the sum of the squared elements in $\mathbf{x}$


## Sum-of-squares | An example  

Let's calculate the model sum-of-squares (SSR) using

$SSR = \sum \left( \widehat{y}_i - \bar{y} \right)^2$


```{r ex_SS_7, echo = TRUE}
## model sum-of-squares
SSR <- t(temp_hat - y_bar) %*% (temp_hat - y_bar)
``` 


## Sum-of-squares | An example  

Let's calculate the error sum-of-squares (SSE) using

$SSE = \sum \left( y_i - \widehat{y}_i \right)^2$


```{r ex_SS_8, echo = TRUE}
## error sum-of-squares
SSE <- t(temp - temp_hat) %*% (temp - temp_hat)
``` 


# QUESTIONS?


## Goodness-of-fit

How about a measure of how well a model fits the data?

> - $SSTO$ measures the variation in $y$ *without* considering $X$

> - $SSE$ measures the reduced variation in $y$ *after* considering $X$

> - Let's consider this reduction in variance as a proportion of the total


## Goodness-of-fit

A common option is the *coefficient of determination* or $(R^2)$ 

$$
R^2 = \frac{SSR}{SSTO} = 1 - \frac{SSE}{SSTO} \\
~ \\
0 < R^2 < 1
$$


## Goodness-of-fit | An example

Let's calculate the $(R^2)$ for our model

```{r ex_SS_9, echo = TRUE}
## coefficient of determination
SSR / SSTO
## or via
1 - SSE / SSTO
``` 


## Degrees of freedom

The number of *independent* elements that are free to vary when estimating quantities of interest


## Degrees of freedom | An example

* Imagine you have 7 hats and you want to wear a different one on each day of the week.

> - On day 1 you can choose any of the 7, on day 2 any of the remaining 6, and so forth

> - When day 7 rolls around, however, you are out of choices: there is only one unworn hat

> - Thus, you had 7 - 1 = 6 days of freedom to choose your hat


## Model in geometric space

```{r model_geom, fig.height=4, fig.width=6, fig.align="center"}
pgram <- data.frame(x = c(0, 3, 10, 7, 0), y = c(0, 3, 3, 0, 0))

## empty plot space
par(mai = rep(0,4), omi = rep(0, 4))
plot(c(0,10), c(0,5), type = "n", xaxt = "n", yaxt = "n",
     bty = "n", xlab = "", ylab = "")
## parallelogram
lines(pgram[,1], pgram[,2])
text(1, 0.5, expression(paste("space spanned by ",bold(X))),
     pos = 4, cex = 1.2, col = "darkgray")
## perpendicular
lines(c(6.6,6.6),c(2,2.3), col = "gray")
lines(c(6.6,7),c(2.3,2.3), col = "gray")
## vectors
## data
arrows(3, 2, 7, 4, length = 0.15, angle = 10, lwd = 2, col ="purple")
text(6.5, 4, expression(bold(y)), cex = 1.5)
## yhat
arrows(3, 2, 7, 2, length = 0.15, angle = 10, lwd = 2, col ="blue")
text(6, 1.6, expression(paste(bold(hat(y))," = ",bold(X),bold(Beta))), cex = 1.5)
## errors
arrows(7, 2, 7, 4, length = 0.15, angle = 10, lwd = 2, col ="red")
text(7.3, 3.5, expression(bold(e)), cex = 1.5)
```

<center>$\mathbf{y}$ is $n$-dim; $\widehat{\mathbf{y}}$ is $k$-dim; $\mathbf{e}$ is $(n-k)$-dim</center>


## Degrees of freedom | Linear models

Beginning with $SSTO$, we have

$$
SSTO = \sum \left( y_i - \bar{y} \right)^2
$$

The data are unconstrained and lie in an $n$-dimensional space, but estimating the mean $(\bar{y})$ from the data costs 1 degree of freedom $(df)$, so

$$
df_{SSTO} = n - 1
$$


## Degrees of freedom | Linear models

For the $SSR$ we have

$$
SSR = \sum \left( \widehat{y}_i - \bar{y} \right)^2
$$

We estimate the data $(\widehat{y})$ with a $k$-dimensional model, but we lose 1 $df$ when estimating the mean, so

$$
df_{SSR} = k - 1
$$


## Degrees of freedom | Linear models

The $SSE$ is analogous

$$
SSE = \sum \left( y_i - \widehat{y}_i \right)^2
$$

The data lie in an $n$-dimensional space and we represent them in a $k$-dimensional subspace, so

$$
df_{SSE} = n - k
$$


## Mean squares

The expectation of the sum-of-squares or "mean square" gives an indication of the variance for the model and errors

A mean square is a sum-of-squares divided by its degrees of freedom

$$
MS = \frac{SS}{df} \\
\Downarrow \\
MSR = \frac{SSR}{k - 1} ~~~ \& ~~~ MSE = \frac{SSE}{n - k} 
$$


## Variance estimates

We are typically interested in two variance estimates:

1. The variance of the residuals $\mathbf{e}$

2. The variance of the model parameters $\mathbf{B}$


## Variance estimates | Residuals

In a least squares context, we assume that the model errors (residuals) are independent and identically distributed with mean 0 and variance $\sigma^2$

The problem is that we don't know $\sigma^2$ and therefore we must estimate it


## Variance estimates | Residuals

If $z_i \sim \text{N}(0, 1)$ then

$$
\sum_{i = 1}^{n} z_i^2 = \mathbf{z}^{\top}\mathbf{z} \sim \chi^2_{n}
$$


## Variance estimates | Residuals

If $z_i \sim \text{N}(0, 1)$ then

$$
\sum_{i = 1}^{n} z_i^2 = \mathbf{z}^{\top}\mathbf{z} \sim \chi^2_{n}
$$

In our linear model, $e_i \sim \text{N}(0, \sigma^2)$ so

$$
\sum_{i = 1}^{n} e_i^2 = \mathbf{e}^{\top}\mathbf{e} \sim \sigma^2 \cdot \chi^2_{n - k} 
$$


## Variance estimates | Residuals

Thus, given 

$$
\mathbf{e}^{\top}\mathbf{e} \sim \sigma^2 \cdot \chi^2_{n - k} \\
\text{E}(\chi^2_{n - k}) = n - k \\
\text{E}(\mathbf{e}^{\top}\mathbf{e}) = SSE
$$

then

$$
SSE = \sigma^2 (n - k) ~ \Rightarrow ~ \sigma^2 = \frac{SSE}{n - k} = MSE
$$


## Variance estimates | Parameters

Recall that our estimate of the model parameters is 

$$
\widehat{\boldsymbol{\beta}} = (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbf{y}
$$


## Variance estimates | Parameters

Estimating the variance of the model parameters $\boldsymbol{\beta}$ requires some linear algebra

For a scalar $z$, if $\text{Var}(z) = \sigma^2$ then $\text{Var}(az) = a^2 \sigma^2$

For a vector $\mathbf{z}$, if $\text{Var}(\mathbf{z}) = \mathbf{\Sigma}$ then $\text{Var}(\mathbf{A z}) = \mathbf{A} \mathbf{\Sigma} \mathbf{A}^{\top}$


## Variance estimates | Parameters

The variance of the parameters is therefore 

$$
\begin{aligned}
\widehat{\boldsymbol{\beta}} &= (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbf{y} \\
   &= \left[ (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \right] \mathbf{y} \\
\end{aligned} \\
\Downarrow \\
\text{Var}(\widehat{\boldsymbol{\beta}}) = \left[ (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \right] \text{Var}(\mathbf{y}) \left[ (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \right]^{\top}
$$


## Variance estimates | Parameters

Recall that we can write our model in matrix form as

$$
\mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \mathbf{e} \\
\mathbf{e} \sim \text{MVN}(\mathbf{0}, \sigma^2 \mathbf{I})
$$


## Variance estimates | Parameters

We can rewrite our model more compactly as

$$
\mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \mathbf{e} \\
\mathbf{e} \sim \text{MVN}(\mathbf{0}, \sigma^2 \mathbf{I}) \\
\Downarrow \\
\mathbf{y} \sim \text{MVN}(\mathbf{X} \boldsymbol{\beta}, \underbrace{\sigma^2 \mathbf{I}}_{\text{Var}(\mathbf{y} | \mathbf{X} \boldsymbol{\beta})}) \\
$$


## Variance estimates | Parameters

Our estimate of $\text{Var}(\widehat{\boldsymbol{\beta}})$ is then  

$$
\begin{aligned}
\text{Var}(\widehat{\boldsymbol{\beta}}) &= \left[ (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \right] \text{Var}(\mathbf{y}) \left[ (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \right]^{\top} \\
  &= \left[ (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \right] \sigma^2 \mathbf{I} \left[ (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \right]^{\top} \\
  &= \sigma^2 (\mathbf{X}^{\top} \mathbf{X})^{-1} (\mathbf{X}^{\top} \mathbf{X}) \left[ (\mathbf{X}^{\top} \mathbf{X})^{-1} \right]^{\top} \\
  &= \sigma^2 (\mathbf{X}^{\top} \mathbf{X})^{-1}
\end{aligned}
$$


## Partitioning total deviations

Here is a plot of some data $y$ and a predictor $x$

```{r intro_SS, echo = FALSE, fig.width = 3.5, fig.height = 3.5, fig.align="center"}
## set random seed for reproducibility
set.seed(514)

## sample size
nn <- 30
## some values for the predictor x
xx <- sort(runif(nn, 0, 20))
## intercept
beta_0 <- 1
## slope
beta_1 <- 0.5
## random errors (mean = 0, SD = 1)
epsilon <- rnorm(nn)
## observed data
yy <- beta_0 + beta_1*xx + epsilon

## fit a regression model
mm <- lm(yy ~ xx)
## get predicted values
yhat <- predict(mm)

par(mai = c(0.3,0.3,0.1,0.1), omi = c(0, 0, 0, 0), cex = 0.7)

## plot the data
plot(xx, yy, type = "p", pch = 19,
     bty = "l", xaxt = "n", yaxt = "n",
     xlab = "",
     ylab = "")
mtext(expression(italic(x)), 1, line=0.2, cex=1.1, at=max(xx))
mtext(expression(italic(y)), 2, line=0.3, cex=1.1, at=max(yy), las=1)

```


## Partitioning total deviations

And let's consider this model: $y_i = \alpha + \beta x_i + e_i$

```{r intro_SS_2, echo = FALSE, fig.width = 3.5, fig.height = 3.5, fig.align="center"}
## plot the data
par(mai = c(0.3,0.3,0.1,0.1),
    omi = c(0, 0, 0, 0),
    cex = 0.7)
plot(xx, yy, type = "p", pch = 19,
     bty = "l", xaxt = "n", yaxt = "n",
     xlab = "",
     ylab = "")
mtext(expression(italic(x)), 1, line=0.2, cex=1.1, at=max(xx))
mtext(expression(italic(y)), 2, line=0.3, cex=1.1, at=max(yy), las=1)
```


## Partitioning total deviations

Here is our model fit to the data

```{r intro_SS_3, echo = FALSE, fig.width = 3.5, fig.height = 3.5, fig.align="center"}
par(mai = c(0.3,0.3,0.1,0.1), omi = c(0, 0, 0, 0), cex = 0.7)

## plot the data
plot(xx, yy, type = "p", pch = 19,
     bty = "l", xaxt = "n", yaxt = "n",
     xlab = "",
     ylab = "")
mtext(expression(italic(x)), 1, line=0.2, cex=1.1, at=max(xx))
mtext(expression(italic(y)), 2, line=0.3, cex=1.1, at=max(yy), las=1)
abline(a = beta_hat[1], b = beta_hat[2], col = "blue", lwd = 2)
```


## Variance estimates | Parameters

Let's think about the variance of $\widehat{\boldsymbol{\beta}}$

$$
\text{Var}(\widehat{\boldsymbol{\beta}}) = \sigma^2 (\mathbf{X}^{\top} \mathbf{X})^{-1}
$$

This suggests that our confidence in our estimate increases with the spread in $\mathbf{X}$


## Effect of $\mathbf{X}$ on parameter precision

Consider these two scenarios where the slope of the relationship is identical

```{r slope_comp, fig.height=3.5, fig.width=7, fig.align="center"}
## new x
xs <- sort(runif(nn, 5, 15))
## random errors (mean = 0, SD = 1)
epsilon <- rnorm(nn)
## observed data
ys <- beta_0 + beta_1*xs + epsilon

m1 <- lm(yy ~ xx)
yy_ci <- predict(m1, interval = "confidence")
m2 <- lm(ys ~ xs)
ys_ci <- predict(m2, interval = "confidence")

XX <- cbind(rep(1, nn), xx)
se_b1 <- round(sqrt(anova(m1)[["Mean Sq"]][2] * solve(t(XX) %*% XX)[2,2]), 3)
XS <- cbind(rep(1, nn), xs)
se_b2 <- round(sqrt(anova(m2)[["Mean Sq"]][2] * solve(t(XS) %*% XS)[2,2]), 3)


par(mfrow = c(1,2), mai = c(0.5, 0.5, 0.1, 0.1), omi = c(0,0,0,0))
## large range in X
plot(xx, yy, pch = 16, xlim = range(xx), ylim = range(yy),
     xlab = expression(italic(y)), ylab = expression(italic(y)))
lines(xx, yy_ci[,"fit"])
text(1, 10.5, substitute(paste(SE(hat(beta))," = ", seb1), list(seb1 = se_b1)), pos = 4)

## small range in X
plot(xs, ys, pch = 16, xlim = range(xx), ylim = range(yy),
     xlab = "", ylab = "",
     yaxt = "n")
lines(xs, ys_ci[,"fit"])
text(1, 10.5, substitute(paste(SE(hat(beta))," = ", seb2), list(seb2 = se_b2)), pos = 4)
```


