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

## Goals for today

* Identify whether a model is linear in the predictors

> - Recognize that linear models can approximate nonlinear functions

> - Understand the difference between categorical and continuous models

> - Recognize the difference between written and coded factors


## Forms of linear models

| Errors | Single random process | Multiple random processes |
|:-------|:--------------:|:------------------:|
| Normal | Linear Model (LM) | Linear Mixed Model (LMM) |
| Non-normal | Generalized Linear Model (GLM) | Generalized Linear Mixed Model (GLMM) |


## Forms of linear models
 
```{r model_diagram_lm, fig.height = 5, fig.width = 7, fig.align = 'center'}
par(mai = rep(0, 4), omi = rep(0, 4), bg = NA)

## order: GLM, LM, GLMM, LMM
xx <- c(1, 2, 2, 3) * 10
yy <- c(2, 1, 3, 2) * 10

plot(c(7, 33), c(7, 33), type = "n", xlim = c(7, 33), ylim = c(7, 33),
     xaxt = "n", yaxt = "n", xlab = "", ylab = "",
     bty = "n")
symbols(xx, yy, circle = rep(2, 4), inches = 0.5, add = TRUE,
        lwd = 2, fg = c("black", "blue", "black", "black"), bg = "white")
text(xx, yy, c("GLM", "LM", "GLMM", "LMM"), cex = 1.5,
     col = c("black", "blue", "black", "black"))
## from LM to GLM
arrows(xx[2]-2, yy[2]+2, xx[1]+2, yy[1]-2, length = 0.2)
text(15, 14, "multiple forms of errors", pos = 2)
## from LM to LMM
arrows(xx[2]+2, yy[2]+2, xx[4]-2, yy[4]-2, length = 0.2)
text(25, 14, "multiple random processes", pos = 4)
## from GLM to GLMM
arrows(xx[1]+2, yy[1]+2, xx[3]-2, yy[3]-2, length = 0.2)
text(15, 26, "multiple random processes", pos = 2)
## from LMM to GLMM
arrows(xx[4]-2, yy[4]+2, xx[3]+2, yy[3]-2, length = 0.2)
text(25, 26, "multiple forms of errors", pos = 4)
```


## What is a linear model?

A relationship that defines a response variable as a linear function of one or more predictor variables


## Which of these are linear models?

<div class="columns-2">

$1) ~ y_i = \delta x_i$

$2) ~ y_i = \alpha + \beta x_i$

$3) ~ y_i = \alpha x_i^{\beta}$

$4) ~ y_i = \alpha + \beta x_i + \gamma z_i$

$5) ~ y_i = \alpha + \beta \frac{1}{x_i}$

$6) ~ y_t = \mu + \phi(y_{t-1} - \mu)$

$7) ~ y_i =   (\alpha + x_i) \beta x_i$

$8) ~ y_i =   \frac{\alpha x_i}{1 + \beta x_i}$

</div>


## Which of these are linear models?

<div class="columns-2">

$1) ~ y_i = \delta x_i$ &#10004;

$2) ~ y_i = \alpha + \beta x_i$ &#10004;

$3) ~ y_i = \alpha x_i^{\beta}$

$4) ~ y_i = \alpha + \beta x_i + \gamma z_i$ &#10004;

$5) ~ y_i = \alpha + \beta \frac{1}{x_i}$

$6) ~ y_t = \mu + \phi(y_{t-1} - \mu)$ &#10004;

$7) ~ y_i =   (\alpha + x_i) \beta x_i$

$8) ~ y_i =   \frac{\alpha x_i}{1 + \beta x_i}$

</div>


## What is a linear model?

*A relationship that defines a response variable as a linear function of one or more predictor variables*

> - characterized by a sum of terms, each of which is the product of a parameter and a single predictor


## Is this a linear model?

$$
y_i = \alpha (1 + \beta x_i)
$$


## Is this a linear model?

$$
y_i = \alpha (1 + \beta x_i)
$$

Yes, *if*

$$
\begin{align}
y_i &= \alpha (1 + \beta x_i) \\
    &= \alpha + \alpha \beta x_i \\
    &= \alpha + \gamma x_i ~~ \text{with} ~~ \gamma = \alpha \beta
\end{align}
$$


## What is a linear model?

*A relationship that defines a response variable as a linear function of one or more predictor variables*

* characterized by a sum of terms, each of which is the product of a parameter and a single predictor

* the predictor can be a transformed variable


## Is this a linear model?

$$
y_i = \alpha + \beta x_i^2 \\
$$


## Linear transformations

Yes, if we allow for a transformation

$$
y_i = \alpha + \beta x_i^2 \\
\Downarrow \\
y_i = \alpha + \beta z_i \\
z_i = x_i^2
$$


## Linear vs nonlinear models

There are only 2 forms of a linear model with 2 parameters

$$
y_i = \alpha + \beta x_i
$$

<center>or</center>

$$
y_i = \beta_1 x_{1,i} + \beta_2 x_{2,i}
$$


## Linear vs nonlinear models

There are *many* forms of nonlinear models with 2 parameters

$$
y_i = \alpha x_i^{\beta} \\
y_i = \alpha + x_i^{\beta} \\
y_i = \alpha^{\beta x_i} \\
y_i = \alpha + \beta \frac{1}{x} \\
\vdots
$$


## Linear vs nonlinear models

In linear models, effect sizes of different predictors are directly comparable

* intercept: units = response (eg, grams)

* slope: units = response per predictor (eg, grams per cm)


## Linear vs nonlinear models

In linear models, effect sizes of different predictors are directly comparable

* intercept: units = response (eg, grams)

* slope: units = response per predictor (eg, grams per cm)

In nonlinear models, common inference tools (*p*-values, confidence intervals) may not be available


## Locally linear models

If we reduce the interval of interest enough, we can approximate a nonlinear function with a linear model

$$
y = x^2 \\
\Downarrow \\
\frac{dy}{dx} = 2x
$$


## Locally linear models | Consider the quadratic $y = \alpha + \beta x + x^2$

```{r nonlinear_scaling, fig.align='center', fig.height=3.5, fig.width=7}
par(mfrow = c(1,2), mai = c(0.9,0.9,0.1,0.1))
alpha <- 2
beta <- 0.5
curve(alpha + beta*x + x^2, -2, 2,
      xlim = c(-2,2), ylim = c(0,8), col = "blue", lwd = 2,
      xlab = expression(italic(x)), ylab = expression(italic(y)))
curve(alpha + beta*x + x^2, 0, 1,
      xlim = c(-2,2), ylim = c(0,8), col = "orange", lwd = 2,
      xlab = expression(italic(x)), ylab = "")
```


## Locally linear models | A stochastic example with $y = 2 + \tfrac{1}{2} x + x^2 + \epsilon_i$

```{r sim_quadratic, echo = TRUE}
set.seed(514)
nn <- 30
alpha <- 2
beta <- 1/2
eps <- rnorm(nn, 0, 1) ## errors ~ N(0,1)
x_all <- runif(nn, -2, 2)
y_all <- alpha + beta*x_all + x_all^2 + eps
x_loc <- x_all[x_all >= 0 & x_all <= 1]
y_loc <- y_all[x_all >= 0 & x_all <= 1]
```


## Locally linear models | A stochastic example with $y = 2 + \tfrac{1}{2} x + x^2 + \epsilon_i$

```{r plot_sim_quadratic, fig.align='center', fig.height=3.5, fig.width=7}
par(mfrow = c(1,2), mai = c(0.9,0.9,0.1,0.1))
plot(x_all, y_all, pch = 16,
     xlim = c(-2,2), ylim = c(0,8), col = "blue",
     xlab = expression(italic(x)), ylab = expression(italic(y)))
plot(x_loc, y_loc, pch = 16,
     xlim = c(-2,2), ylim = c(0,8), col = "orange",
     xlab = expression(italic(x)), ylab = "")
```


## Linear model for size of fish

In **R**, we can use `lm()` to fit linear regression models

$y_i = \alpha + \beta x_i + e_i$

`lm(y ~ x)`

(notice that the intercept $\alpha$ is implicit here)


## Linear model for size of fish

In **R**, we use `summary()` to get info about a fitted model

`fitted_regr_model <- lm(L10_mass ~ L10_length)` 

`summary(fitted_regr_model)` 


## Locally linear models {.smaller}

```{r fit_quad_1, echo = TRUE}
## model 1: full dataset
fit_1 <- lm(y_all ~ x_all)
summary(fit_1)
```


## Locally linear models {.smaller}

```{r fit_quad_2, echo = TRUE}
## model 2: "local" data
fit_2 <- lm(y_loc ~ x_loc)
summary(fit_2)
```


## Linear model for size of fish

In **R**, we use `coef()` to extract the intercept(s) and slope(s)

`fitted_regr_model <- lm(y ~ x)` 

`coef(fitted_regr_model)` 


## Locally linear models

```{r coef_quad_2, echo = TRUE}
## intercept and slope for model 2
coef(fit_2)
```


## Locally linear models

```{r coef_quad_2b, echo = TRUE}
## intercept and slope for model 2
coef(fit_2)
```

True model: $y = \tfrac{1}{2} + 2 x + x^2$

Estimate: $\hat{y} \approx 1.1 + 2.7 x + 0 x^2$


## {.flexbox .vcenter .bigger}

<font size="10">
<center>Linear models can be *good approximations* to nonlinear functions</center>
</font>


# QUESTIONS?


## {.flexbox .vcenter .bigger}

<font size="10">
<center>Common forms for linear models</center>
</font>


## A simple starting point

### Data = (Deterministic part) + (Stochastic part)


## Types of linear models

We classify linear models by the form of their deterministic part

Discrete predictor $\rightarrow$ **AN**alysis **O**f **VA**riance (ANOVA)

Continuous predictor $\rightarrow$ Regression

Both $\rightarrow$ **AN**alysis of **COVA**riance (ANCOVA)


## Possible models for growth of fish

| Model | Description |
|:------|:-----------:|
| $\text{growth}_i = \alpha + \beta \text{species}_i + \epsilon_i$ | 1-way ANOVA |
| $\text{growth}_i = \alpha + \beta_{1,\text{species}} + \beta_{2,\text{tank}} + \epsilon_i$ | 2-way ANOVA |
| $\text{growth}_i = \alpha + \beta \text{ration}_i + \epsilon_i$ | simple linear regression |
| $\text{growth}_i = \alpha + \beta_1 \text{ration}_i + \beta_2 \text{temperature}_i + \epsilon_i ~ ~$ | multiple regression |
| $\text{growth}_i = \alpha + \beta_{1,\text{species}} + \beta_2 \text{ration}_i + \epsilon_i$ | ANCOVA |


## Example | Fish growth during an experiment

* A biologist at the WA Dept of Fish & Wildlife contacts you for help with an experiment

* She wants to know how growth of hatchery salmon is affected by ration size

* She sends you a spreadsheet with 2 cols:  

    1. fish growth (mm)  
    2. ration size (2g, 4g, 6g)


## Fish growth data

```{r fish_data, fig.align='center', fig.height=4, fig.width=5}
## sample size
nn <- 30
## groups
pp <- 3
## global intercept
alpha <- 5
## offsets
beta_1 <- c(1,2,3)*5
## slope
beta_2 <- 2
## vector of lienar parameters
BETA <- matrix(c(alpha, beta_1, beta_2), ncol = 1)
## linear predictors
x_avg <- rep(1, nn*pp)
x_int <- matrix(c(rep(c(rep(1,nn),
                        rep(0,nn*pp)),2),
                  c(rep(1,nn))),
                nn*pp, pp)
x_cov <- runif(nn*pp, 0, 12)
## matrix of predictors
xx <- cbind(x_avg, x_int, x_cov)
## Gaussian errors
ee <- rnorm(nn*pp, 0, 2)
## simulated data
yy <- xx %*% BETA + ee
## bin data for ANOVA
vv <- cbind(x_cov, yy)
v1 <- vv[x_cov <=4, ]
v2 <- vv[x_cov > 4 & x_cov <=8, ]
v3 <- vv[x_cov > 8 & x_cov <=12, ]
## plot all data
par(mai = c(0.9,0.9,0.1,0.1))
## low
plot(rep(1, nrow(v1)), v1[,2], pch = 16, col = "red",
     xlim = c(0.5,3.5), ylim = range(yy),
     xaxt = "n",
     cex.lab = 1.5, cex.axis = 1.2,
     xlab = "Ration size (g)", ylab = "Growth (mm)")
## med
points(rep(2, nrow(v2)), v2[,2], pch = 16, col = "blue")
# high
points(rep(3, nrow(v3)), v3[,2], pch = 16, col = "orange")
axis(1, at = seq(3), cex.axis = 1.2,
     labels = c("Low (2)", "Med (6)", "High (10)"))
```


## ANOVA model

```{r sim_anova, fig.align='center', fig.height=4, fig.width=5}
## plot all data
par(mai = c(0.9,0.9,0.1,0.1))
## low
plot(rep(1, nrow(v1)), v1[,2], pch = 16, col = "red",
     xlim = c(0.5,3.5), ylim = range(yy),
     xaxt = "n",
     cex.lab = 1.5, cex.axis = 1.2,
     xlab = "Ration size (g)", ylab = "Growth (mm)")
points(1, mean(v1[,2]), col = "red", pch = "-", cex = 5)
## med
points(rep(2, nrow(v2)), v2[,2], pch = 16, col = "blue")
points(2, mean(v2[,2]), col = "blue", pch = "-", cex = 5)
# high
points(rep(3, nrow(v3)), v3[,2], pch = 16, col = "orange")
points(3, mean(v3[,2]), col = "orange", pch = "-", cex = 5)
axis(1, at = seq(3), cex.axis = 1.2,
     labels = c("Low (2)", "Med (6)", "High (10)"))
```

$\text{growth}_i = \alpha + \beta_{\text{ration}} + \epsilon_i$


## More info arrives

It turns out that targeting the exact ration is hard, but they know how much each fish ate during the trial


## Continuous predictor

```{r sim_regr, fig.align='center', fig.height=4, fig.width=5}
par(mai = c(0.9,0.9,0.1,0.1))
plot(v1[,1], v1[,2], pch = 16, col = "red",
     cex.lab = 1.5, cex.axis = 1.2,
     xlim = range(x_cov), ylim = range(yy),
     xlab = "Ration size (g)", ylab = "Growth (mm)")
lines(c(0,4), rep(mean(v1[,2]), 2), col = "red", lwd = 3)
## med
points(v2[,1], v2[,2], pch = 16, col = "blue")
lines(c(4,8), rep(mean(v2[,2]), 2), col = "blue", lwd = 3)
# high
points(v3[,1], v3[,2], pch = 16, col = "orange")
lines(c(8,12), rep(mean(v3[,2]), 2), col = "orange", lwd = 3)
```

$\text{growth}_i = \alpha + \beta_{\text{ration}} + \epsilon_i$


## Continuous predictor

```{r sim_regr_all, fig.align='center', fig.height=4, fig.width=5}
par(mai = c(0.9,0.9,0.1,0.1))
plot(x_cov, yy, pch = 16, 
     cex.lab = 1.5, cex.axis = 1.2,
     xlab = "Ration size (g)", ylab = "Growth (mm)")
```


## Linear regression

```{r sim_regr_all_2, fig.align='center', fig.height=4, fig.width=5}
par(mai = c(0.9,0.9,0.1,0.1))
plot(x_cov, yy, pch = 16, 
     cex.lab = 1.5, cex.axis = 1.2,
     xlab = "Ration size (g)", ylab = "Growth (mm)")
abline(a = 15, b = 2)
```

$\text{growth}_i = \alpha + \beta \text{ration}_i + \epsilon_i$


## More info arrives

It also turns out that there are 3 lineages of fish in the trials


## Continuous & discrete predictors

```{r sim_ancova, fig.align='center', fig.height=4, fig.width=5}

## all data
par(mai = c(0.9,0.9,0.1,0.1))
plot(x_cov[1:nn], yy[1:nn], pch = 16, col = "red", ylim = range(yy),
     cex.lab = 1.5, cex.axis = 1.2,
     xlab = "Ration size (g)", ylab = "Growth (mm)")
points(x_cov[1:nn+nn], yy[1:nn+nn], pch = 16, col = "blue")
points(x_cov[1:nn+nn*2], yy[1:nn+nn*2], pch = 16, col = "orange")
```


## ANCOVA

```{r sim_ancova_fit, fig.align='center', fig.height=4, fig.width=5}

## all data
par(mai = c(0.9,0.9,0.1,0.1))
plot(x_cov[1:nn], yy[1:nn], pch = 16, col = "red", ylim = range(yy),
     cex.lab = 1.5, cex.axis = 1.2,
     xlab = "Ration size (g)", ylab = "Growth (mm)")
points(x_cov[1:nn+nn], yy[1:nn+nn], pch = 16, col = "blue")
points(x_cov[1:nn+nn*2], yy[1:nn+nn*2], pch = 16, col = "orange")

abline(a = 10, b = 2, col = "red")
abline(a = 15, b = 2, col = "blue")
abline(a = 20, b = 2, col = "orange")
```

$\text{growth}_i = \alpha + \beta_{1,\text{lineage}} + \beta_2 \text{ration}_i + \epsilon_i$


## Notation for categorical effects

Here we have specified categorical effects in AN(C)OVA models as discrete parameters

For example, for a one-way ANOVA with 3 factors

$$
y_i = \alpha + \beta_j + \epsilon_i
$$

the definition of $\beta_j$ is

$$
\beta_j = 
\left\{
\begin{matrix}
\beta_1 ~ \text{if factor 1} \\
\beta_2 ~ \text{if factor 2} \\
\beta_3 ~ \text{if factor 3}
\end{matrix}
\right.
$$


## Notation for categorical effects

In practice, we will use a combination of -1/0/1 predictors, so our model becomes

$$
y_i = \alpha + \beta_1 x_{1,i} + \beta_2 x_{2,i} + \beta_3 x_{3,i} + \epsilon_i
$$

and each of the $x_{j,i}$ indicates whether the $i^{th}$ observation was assigned factor $j$

<br>

(We'll visit this again when we discuss design matrices)

