lm.Rmd

---
title: "Linear models"
author: "Francisco Rodríguez-Sánchez"
institute: "https://frodriguezsanchez.net"
aspectratio: 43  # use 169 for wide format
fontsize: 10pt
output: 
  binb::metropolis:
    keep_tex: no
    incremental: yes
    fig_caption: no
    pandoc_args: ["--lua-filter=hideslide.lua"]
urlcolor: blue
linkcolor: blue
header-includes:
  - \definecolor{shadecolor}{RGB}{230,230,230}
  # - \setbeamercolor{frametitle}{bg=black}
---


```{r knitr_setup, include=FALSE, cache=FALSE}

library("knitr")

### Chunk options ###

## Text results
opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE, size = 'tiny')

## Code decoration
opts_chunk$set(tidy = FALSE, comment = NA, highlight = TRUE, prompt = FALSE, crop = TRUE)

# ## Cache
# opts_chunk$set(cache = TRUE, cache.path = "knitr_output/cache/")

# ## Plots
# opts_chunk$set(fig.path = "knitr_output/figures/")
opts_chunk$set(fig.align = 'center', out.width = '90%')

### Hooks ###
## Crop plot margins
knit_hooks$set(crop = hook_pdfcrop)

## Reduce font size
## use tinycode = TRUE as chunk option to reduce code font size
# see http://stackoverflow.com/a/39961605
knit_hooks$set(tinycode = function(before, options, envir) {
  if (before) return(paste0("\n \\", options$size, "\n\n"))
  else return("\n\n \\normalsize \n")
  })

```




## Example dataset: forest trees

- Download [this dataset](https://raw.githubusercontent.com/Pakillo/LM-GLM-GLMM-intro/trees/data/trees.csv) (or the entire [zip file](https://github.com/Pakillo/LM-GLM-GLMM-intro/raw/trees/datasets.zip))

- Import:

```{r}
trees <- read.csv("data/trees.csv")
head(trees)
```


## Questions

- What is the relationship between DBH and height?

\vspace{10mm}

- Do taller trees have bigger trunks?

\vspace{10mm}

- Can we predict height from DBH? How well?



# Always plot your data first!


## Always plot your data first!

```{r echo=FALSE}
include_graphics("images/anscombe.png")
```



## Exploratory Data Analysis (EDA)

Outliers

```{r indexplot}
plot(trees$height)
```



## Outliers impact on regression

```{r out.width="70%", echo=FALSE}
include_graphics("images/reg_outliers.png")
```

See http://rpsychologist.com/d3/correlation/


## Histogram of response variable

```{r histog}
hist(trees$height)
```


## Histogram of predictor variable

```{r}
hist(trees$dbh)
```

## Scatterplot

```{r scatterplot}
plot(height ~ dbh, data = trees, las = 1)
```

## Scatterplot

```{r echo=3:4, tinycode = TRUE}
library(ggplot2)
theme_set(theme_minimal(base_size = 18))
ggplot(trees) +
  geom_point(aes(dbh, height)) 
```




# Model fitting

## Now fit model

Hint: `lm`


## Now fit model

Hint: `lm`

\vspace{5mm}

```{r lm_trees}
m1 <- lm(height ~ dbh, data = trees)
```

\vspace{5mm}

which corresponds to

$$
  \begin{aligned}  
  Height_{i} = a + b \cdot DBH_{i} + \varepsilon _{i} \\  
  \varepsilon _{i}\sim N\left( 0,\sigma^2 \right) \\  
  \end{aligned}  
$$


## Package `equatiomatic` returns model structure

```{r}
library("equatiomatic")
m1 <- lm(height ~ dbh, data = trees)
equatiomatic::extract_eq(m1)
```


# Model interpretation


## What does this mean?

\footnotesize

```{r summary_lm, echo=TRUE}
summary(m1)
```


## Remember that in a Normal distribution

```{r echo=FALSE}
include_graphics("images/gaussian.png")
```


## Estimated distribution of the **intercept** parameter

```{r echo=FALSE}
library(arm)
library(ggplot2)

coefs <- as.data.frame(coef(sim(m1)))
names(coefs) <- c("intercept", "slope")

ggplot(coefs) +
  geom_density(aes(intercept), fill = "grey80") +
  xlim(-1, 21) +
  geom_vline(xintercept = 0)
```


## Estimated distribution of the **slope** parameter

```{r echo=FALSE}
ggplot(coefs) +
  geom_density(aes(slope), fill = "grey80") +
  xlim(-0.1, 0.7) +
  geom_vline(xintercept = 0)
```


## Distribution of residuals

```{r echo=FALSE}
res <- data.frame(residual = residuals(m1))
ggplot(res) +
  geom_density(aes(residual), fill = "grey80") +
  geom_vline(xintercept = 0) +
  annotate("text", x = 3, y = 0.015, label = "SD = 4", size = 5)
```


## Degrees of freedom

DF = n - p

n = sample size

p = number of estimated parameters




## R-squared

Proportion of 'explained' variance

$R^{2} = 1 - \frac{Residual Variation}{Total Variation}$



## Adjusted R-squared

Accounts for model complexity (number of parameters)

$1 - (1 - R^2) \frac{n - 1}{n - p - 1}$


## Quiz 

https://pollev.com/franciscorod726



## Retrieving model coefficients

```{r echo = TRUE}
coef(m1)
```

## Confidence intervals for parameters

```{r echo = TRUE}
confint(m1)
```


## Tidy up model coefficients with broom

\footnotesize

```{r}
library("broom")
tidy(m1)
```

\vspace{5mm}

```{r}
glance(m1)
```

https://broom.tidymodels.org/


## Retrieving model parameters with `parameters` package

\footnotesize

```{r}
library("parameters")
parameters(m1)
```

https://easystats.github.io/parameters/



## Understanding the fitted effects with `effects` package

\footnotesize
```{r message=FALSE}
library("effects")
summary(allEffects(m1))
```



# Communicating results

## Avoid dichotomania of statistical significance

```{r echo=FALSE}
include_graphics("images/nature_significance.PNG")
```

- 'Never conclude there is ‘no difference’ or ‘no association’ just because p > 0.05 or CI includes zero'

- Estimate and communicate effect sizes and their uncertainty

- https://doi.org/10.1038/d41586-019-00857-9


## Communicating results

We found a significant relationship between DBH and Height (p<0.05).

We found a ~~significant~~ positive relationship between DBH and Height ~~(p<0.05)~~ (b = 0.61, SE = 0.01).


## Models that describe themselves

```{r results="asis"}
library("report")
report(m1)
```

https://easystats.github.io/report/


## Generating table with model results: `xtable`

```{r echo=TRUE, results='asis'}
library("xtable")
xtable(m1, digits = 2)
```


## Generating table with model results: `texreg`

```{r echo=TRUE, results='asis'}
library("texreg")
texreg(m1, single.row = TRUE)
```


## Generating table with model results: `modelsummary`

```{r results='asis'}
library("modelsummary")
modelsummary(m1, output = "markdown")
```

https://vincentarelbundock.github.io/modelsummary/


## Generating table with model results: `gtsummary`

```{r }
library("gtsummary")
tbl_regression(m1, intercept = TRUE)
```

https://www.danieldsjoberg.com/gtsummary



# Visualising fitted model


## Plot model: `effects` package

```{r echo = TRUE}
library("effects")
plot(allEffects(m1))
```


## Plot model: `visreg`

```{r visreg}
library("visreg")
visreg(m1)
```


## `visreg` can use ggplot2 too

```{r out.width="80%"}
visreg(m1, gg = TRUE) + theme_bw()
```

https://pbreheny.github.io/visreg


## Plot model: `sjPlot`

```{r out.width="70%"}
library("sjPlot")
plot_model(m1, type = "eff")
```

https://strengejacke.github.io/sjPlot


## Plot model: `see`

\footnotesize

```{r out.width="80%"}
library("see")
plot(parameters(m1), show_intercept = TRUE) +
  labs(title = "Height ~ Diameter")   # ggplot2
```


## Plot parameters' estimated distribution: `see`

```{r out.width="70%"}
plot(simulate_parameters(m1)) +
  labs(title = "Density of the slope parameter")
```




# Model checking


## Linear model assumptions

::: nonincremental :::

- **Linearity** (transformations, GAM...)

\vspace{2mm}
  
- **Residuals**:
    - Independent
    - Equal variance
    - Normal

\vspace{2mm}

- Negligible **measurement error** in predictors

:::

```{r out.width="50%", echo = FALSE}
include_graphics("images/lm_resid_assump.png")
```


## Are residuals normal? 

```{r resid_hist, out.width="70%"}
hist(residuals(m1))
```

SD = 4.09


## Model checking: `plot(model)`

```{r echo=FALSE}
def.par <- par(no.readonly = TRUE)
layout(matrix(1:4, nrow=2))
plot(m1)
par(def.par)
```


## Model checking with `performance` package

```{r}
library("performance")
check_model(m1)
```



# Using model for prediction


## How good is the model in predicting tree height?

`fitted` gives expected value for each observation

\vspace{4mm}

```{r}
trees$height.pred <- fitted(m1)
trees$resid <- residuals(m1)
head(trees)
```


## Calibration plot: Observed vs Predicted values

```{r obs_pred, echo=FALSE}
plot(height ~ height.pred, data = trees,
     xlab = "Tree height (predicted)", ylab = "Tree height (observed)", 
     las = 1, xlim = c(10,60), ylim = c(10,60))
abline(a = 0, b = 1)
```




## Making predictions for new data

Q: Expected tree height if DBH = 39 cm?

\vspace{3mm}

```{r}
new.dbh <- data.frame(dbh = c(39))
predict(m1, new.dbh, se.fit = TRUE)
```


## Confidence vs Prediction Intervals

Q: Expected tree height if DBH = 39 cm?

```{r}
predict(m1, new.dbh, interval = "confidence")
```

```{r}
predict(m1, new.dbh, interval = "prediction")
```


## Confidence vs Prediction Intervals

:::::::::::::: {.columns align=center}

::: {.column width="50%"}

```{r echo=FALSE}
library("modelbased")
conf <- estimate_relation(m1, data = NULL)
conf$obs <- trees$height

ggplot(conf) +
  aes(x = dbh, y = obs) +
  geom_point(size = 0.5, alpha = 0.3) +
  geom_ribbon(aes(ymin = CI_low, ymax = CI_high), fill = "blue", alpha = 0.2) +
  geom_line(aes(x = dbh, y = Predicted), colour = "blue", size = 0.5) +
  labs(x = "diameter", y = "height", title = "Confidence interval") +
  theme_minimal(base_size = 18)
```

:::

::: {.column width="50%"}

```{r echo=FALSE}

pred <- estimate_response(m1)
pred$obs <- trees$height

ggplot(pred) +
  aes(x = dbh, y = obs) +
  geom_point(size = 0.5, alpha = 0.5) +
  geom_ribbon(aes(ymin = CI_low, ymax = CI_high), fill = "blue", alpha = 0.2) +
  geom_line(aes(x = dbh, y = Predicted), colour = "blue", size = 2) +
  labs(x = "diameter", y = "height", title = "Prediction interval") +
  theme_minimal(base_size = 18)
```

:::
::::::::::::::



## Workflow

- **Visualise data**

\vspace{5mm}

- **Understand fitted model** (`summary`, `allEffects`...)

\vspace{5mm}

- **Visualise model** (`plot(allEffects)`, `visreg`, `see`, `plot_model`...)

\vspace{5mm}

- **Check model** (`plot`, `check_model`, calibration plot...)

\vspace{5mm}

- **Predict** (`fitted`, `predict`)



`



# Categorical predictors (factors)


## Q: Does tree height vary with sex?

```{r boxplot}
plot(height ~ as.factor(sex), data = trees)
```



## Model height ~ sex

\footnotesize

```{r echo=1}
m2 <- lm(height ~ sex, data = trees)
summary(m2)
```



## Linear model with categorical predictors

```{r eval=FALSE}
m2 <- lm(height ~ sex, data = trees)
```

corresponds to 

$$
  \begin{aligned} 
  Height_{i} = a + b_{male} + \varepsilon _{i} \\
  \varepsilon _{i}\sim N\left( 0,\sigma^2 \right) \\
  \end{aligned} 
$$



## Model height ~ sex

```{r echo=1}
m2 <- lm(height ~ sex, data = trees)
summary(m2)
```


## Quiz

https://pollev.com/franciscorod726


## Let's read the model report...

```{r results='asis'}
report(m2)
```



## Estimated distribution of the **intercept** parameter

Intercept = Height of females

```{r echo=FALSE}
coefs <- as.data.frame(coef(sim(m2)))
names(coefs) <- c("intercept", "slope")

ggplot(coefs) +
  geom_density(aes(intercept), fill = "grey80") +
  xlim(-1, 40) +
  geom_vline(xintercept = 0)
```


## Estimated distribution of the *beta* parameter

*beta* = height difference of males vs females


```{r echo=FALSE}
ggplot(coefs) +
  geom_density(aes(slope), fill = "grey80") +
  xlim(-3, 2) +
  geom_vline(xintercept = 0)
```





## Analysing differences among factor levels

```{r}
library("modelbased")
estimate_means(m2)
```


## Analysing differences among factor levels

\footnotesize

```{r}
estimate_contrasts(m2)
```



## Plot

```{r}
plot(allEffects(m2))
```



## Plot (visreg)

```{r }
visreg(m2)
```


## Plot model (sjPlot)

```{r out.width="60%"}
plot_model(m2, type = "eff")
```


## Model checking: residuals

```{r}
hist(resid(m2))
```


## Model checking: residuals

```{r echo=FALSE}
def.par <- par(no.readonly = TRUE)
layout(matrix(1:4, nrow=2))
plot(m2)
par(def.par)
```


## Model checking

```{r}
library("performance")
check_model(m2)
```





# Q: Does height differ among field sites?

## Plot data first

```{r}
plot(height ~ site, data = trees)
```


## Linear model with categorical predictors

```{r }
m3 <- lm(height ~ site, data = trees)
```


$$
  \begin{aligned} 
  y_{i} = a + b_{site2} + c_{site3} + d_{site4} + e_{site5} +...+ \varepsilon _{i} \\   
  \varepsilon _{i}\sim N\left( 0,\sigma^2 \right) \\
  \end{aligned} 
$$





## Model Height ~ site

**All right here?**

\tiny

```{r echo=1}
m3 <- lm(height ~ site, data = trees)
summary(m3)
```


## Let's check model structure with `equatiomatic`

```{r}
extract_eq(m3)
```

## site is a factor!

```{r}
trees$site <- as.factor(trees$site)
```


## Let's check model structure with `equatiomatic`

```{r }
m3 <- lm(height ~ site, data = trees)
extract_eq(m3)
```


## Model Height ~ site

\scriptsize
```{r echo=FALSE}
m3 <- lm(height ~ site, data = trees)
summary(m3)
```
\normalsize



## Estimated parameter distributions

```{r}
plot(simulate_parameters(m3), stack = FALSE)
```




## Analysing differences among factor levels

```{r}
library("modelbased")
estimate_means(m3)
```


## Analysing differences among factor levels

\tiny

For finer control see `emmeans` package

```{r}
estimate_contrasts(m3)
```



## Presenting model results

```{r }
kable(xtable::xtable(m3), digits = 2)
```


## Estimated tree heights for each site

\scriptsize
```{r}
summary(allEffects(m3))
```
\normalsize






## Plot

```{r}
plot(allEffects(m3))
```



## Plot (visreg)

```{r }
visreg(m3)
```


## Plot model (sjPlot)

```{r out.width="70%"}
plot_model(m3, type = "eff")
```



## Model checking: residuals

```{r echo=FALSE}
def.par <- par(no.readonly = TRUE)
layout(matrix(1:4, nrow = 2))
plot(m3)
par(def.par)
```


## Model checking: residuals

```{r}
check_model(m3)
```





# Combining continuous and categorical predictors


## Predicting tree height based on dbh and site

```{r eval=FALSE}
lm(height ~ site + dbh, data = trees)
```

corresponds to 

$$
  \begin{aligned} 
  y_{i} = a + b_{site2} + c_{site3} + d_{site4} + e_{site5} +...+ k \cdot DBH_{i} + \varepsilon _{i} \\ 
  \varepsilon _{i}\sim N\left( 0,\sigma^2 \right) \\
  \end{aligned} 
$$


## Predicting tree height based on dbh and site

\scriptsize
```{r echo = FALSE}
m4 <- lm(height ~ site + dbh, data = trees)
summary(m4)
```
\normalsize


## Presenting model results

\footnotesize

```{r echo=TRUE}
parameters(m4)
```



## Estimated tree heights for each site

\scriptsize
```{r}
summary(allEffects(m4))
```
\normalsize


## Plot

```{r}
plot(allEffects(m4))
```



## Plot (visreg)

```{r echo=2}
par(mfcol = c(1, 2))
visreg(m4)
dev.off()
```


## Plot model (sjPlot)

```{r echo=TRUE, out.width='30%'}
plot_model(m4, type = "eff")

```

## Plot model (sjPlot)

```{r echo=TRUE}
plot_model(m4, type = "est")
```


## Plot model (see)

```{r}
plot(parameters(m4))
```


## We have fitted model w/ many intercepts and single slope

```{r echo=FALSE, eval=FALSE}
plot(height ~ dbh, data = trees, las = 1)
abline(a = coef(m4)[1], b = coef(m4)[11])
for (i in 2:10) {
  abline(a = coef(m4)[1] + coef(m4)[i], b = coef(m4)[11])
}
```

```{r echo=FALSE}
## ggplot with different colour for each site
ggplot(trees) +
  aes(x = dbh, y = height, colour = site) +
  geom_point() +
  geom_abline(intercept = coef(m4)[1], slope = coef(m4)[11]) +
  geom_abline(intercept = coef(m4)[1] + coef(m4)[2], slope = coef(m4)[11]) +
  geom_abline(intercept = coef(m4)[1] + coef(m4)[3], slope = coef(m4)[11]) +
  geom_abline(intercept = coef(m4)[1] + coef(m4)[8], slope = coef(m4)[11])
```


## Slope is the same for all sites

\scriptsize

```{r}
estimate_slopes(m4)
```



## Model checking: residuals

```{r echo=FALSE}
def.par <- par(no.readonly = TRUE)
layout(matrix(1:4, nrow=2))
plot(m4)
par(def.par)
```


## Model checking: residuals

```{r}
check_model(m4)
```



## How good is this model? Calibration plot

```{r}
trees$height.pred <- fitted(m4)
plot(trees$height.pred, trees$height, xlab = "Tree height (predicted)", ylab = "Tree height (observed)", las = 1, xlim = c(10,60), ylim = c(10,60))
abline(a = 0, b = 1)
```


## *Posterior* predictive checking

Simulating response data from fitted model (`yrep`)

and comparing with observed response (`y`)

```{r out.width="70%"}
performance::pp_check(m4)
```


```{r eval=FALSE, echo=FALSE}
library(bayesplot)
sims <- simulate(m4, nsim = 100)
ppc_dens_overlay(trees$height, yrep = t(as.matrix(sims)))
```




## Using model for prediction

Expected height of 10-cm diameter tree in each site?

```{r}
trees.10cm <- data.frame(site = as.factor(1:10),
                        dbh = 10)
trees.10cm
```

## Using model for prediction

Confidence interval

```{r}
predict(m4, newdata = trees.10cm, interval = "confidence")
```


## Using model for prediction

Prediction interval (accounting for residual variance)

```{r}
predict(m4, newdata = trees.10cm, interval = "prediction")
```


## Using model for prediction

Prediction interval (99%)

```{r}
predict(m4, newdata = trees.10cm, interval = "prediction", 
        level = 0.99)
```





# Q: Does allometric relationship between Height and Diameter vary among sites?



## Model with interactions

\tiny

```{r echo=FALSE}
m5 <- lm(height ~ site*dbh, data = trees)
summary(m5)
```

\normalsize


## Does slope vary among sites?

```{r}
visreg(m5, xvar = "dbh", by = "site")
```





## Extra exercises

- [paperplanes](https://cran.r-project.org/package=paperplanes): How does flight distance differ with age, gender or paper type?

- [mammal sleep](https://ggplot2.tidyverse.org/reference/msleep.html): Are sleep patterns related to diet?

- iris: Predict petal length ~ petal width and species

- [Penguins data](https://cran.r-project.org/package=palmerpenguins): Body mass ~ Flipper length, Bill length ~ Bill depth, differences across sites...

- [racing pigeons](http://blog.yhat.com/posts/7-funny-datasets.html): is speed related to sex?