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

\Large

- What is the relationship between DBH and height?

\vspace{5mm}

- Do taller trees have bigger trunks?

\vspace{5mm}

- 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(x = dbh, y = 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)
```


```{r}
equatiomatic::extract_eq(m1, use_coefs = TRUE)
```


## To preview LaTeX:

```{r eval=FALSE}
library(texPreview)
tex_preview(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

\scriptsize

```{r echo=FALSE, message=FALSE}
parameters::parameters(m1)[1,]
```

```{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

\scriptsize

```{r echo=FALSE, message=FALSE}
parameters::parameters(m1)[2,]
```

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

```{r echo=FALSE, eval=FALSE}
library(easystats)
plot(simulate_parameters(m1)) +
  labs(title = "Density of the slope parameter")
```

## Distribution of residuals

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


## Degrees of freedom

\LARGE

DF = n - p

n = sample size

p = number of estimated parameters




## R-squared

\LARGE

Proportion of 'explained' variance

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



## Adjusted R-squared

\LARGE

Accounts for model complexity 

(number of parameters)

\vspace{5mm}

$R^2_{adj} = 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/


::: hide :::

## 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)**.

- (add p-value if you wish)


## Models that describe themselves

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

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



::: hide :::
## 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: `gtsummary`

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

https://www.danieldsjoberg.com/gtsummary




## Generating table with model results: `modelsummary`

```{r results='asis'}
library("modelsummary")
modelsummary(m1, output = "markdown")  # Word, PDF, PowerPoint, png...
```

https://modelsummary.com/


## Generating table with model results: `modelsummary`

```{r results='asis'}
modelsummary(m1, fmt = 2, 
             estimate = "{estimate} ({std.error})", 
             statistic = NULL,
             gof_map = c("nobs", "r.squared", "rmse"),
             output = "markdown")  # Word, PDF, PowerPoint, png...
```




# 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


## easystats

```{r}
library("easystats")
plot(estimate_expectation(m1))
```


## Plot model: `sjPlot`

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

https://strengejacke.github.io/sjPlot


## ggeffects

\scriptsize

```{r}
library("ggeffects")
```

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

::: {.column width="40%"}
```{r}
mydf <- ggpredict(m1, terms = "dbh")
dplyr::glimpse(mydf, width = 40)
```

:::

::: {.column width="60%" }
```{r}
ggplot(mydf, aes(x, predicted)) +
  geom_line() +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), 
              alpha = 0.1)
```

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





## modelsummary

```{r}
modelplot(m1)
```


## Plot model parameters with easystats (`see` package)

\footnotesize

```{r out.width="80%"}
library("easystats")
plot(parameters(m1), show_intercept = TRUE, show_labels = TRUE)
```


## Plot parameters' estimated distribution

```{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` (easystats)

```{r out.width="80%"}
library("easystats")
check_model(m1)
```

\tiny https://easystats.github.io/performance/articles/check_model.html



## A dashboard to explore the full model

```{r eval=FALSE}
library("easystats")
model_dashboard(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_expectation(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_prediction(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)
```

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


# Making predictions with easystats

## Estimate expected values

```{r echo=1}
pred <- estimate_expectation(m1)
head(pred)
```


## Expected values given DBH

```{r}
plot(estimate_expectation(m1))
```



## Calibration plot: observed vs predicted

\footnotesize

```{r}
pred$height.obs <- trees$height
plot(height.obs ~ Predicted, data = pred, xlim = c(15, 60), ylim = c(15, 60))
abline(a = 0, b = 1)
```


## Estimate prediction interval

Accounting for residual variation!

```{r}
pred <- estimate_prediction(m1)
head(pred)
```

## Confidence vs Prediction interval

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

::: {.column width="50%"}
```{r}
plot(estimate_expectation(m1))
```
:::

::: {.column width="50%" }
```{r}
plot(estimate_prediction(m1))
```
:::
::::::::::::::


## Make predictions for new data

```{r}
estimate_expectation(m1, data = data.frame(dbh = 39))
```

```{r}
estimate_prediction(m1, data = data.frame(dbh = 39))
```




## Workflow

- **Visualise data**

\vspace{5mm}

- **Understand fitted model** (`summary`)

\vspace{5mm}

- **Visualise model** (`visreg`...)

\vspace{5mm}

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

\vspace{5mm}

- **Predict** (`predict`, `estimate_expectation`, `estimate_prediction`)







# Categorical predictors (factors)


## Q: Does tree height vary with sex?

```{r boxplot}
boxplot(height ~ 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

\footnotesize

```{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**

\scriptsize

```{r echo=FALSE, message=FALSE}
parameters(m2)[1,]
```


```{r echo=FALSE, out.width="65%"}
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

\scriptsize

```{r echo=FALSE, message=FALSE}
parameters(m2)[2,]
```

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





## Analysing differences among factor levels

```{r}
library("easystats")  # modelbased package
estimate_means(m2)
```


## Analysing differences among factor levels

\footnotesize

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


# Visualising the fitted model


## Plot (effects)

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



## Plot (visreg)

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


## Plot (easystats)

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


## Plot model (sjPlot)

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


# Model checking

## 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("easystats")
check_model(m2)
```


## Model dashboard

```{r eval=FALSE}
model_dashboard(m2)
```


# Q: Does height differ among field sites?

## Quiz 

https://pollev.com/franciscorod726

## 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?**

\footnotesize

```{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)
```




## Estimated tree heights for each site

```{r}
estimate_means(m3)
```


## Plot estimated tree heights for each site

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

## Analysing differences among factor levels

\tiny

For finer control see `emmeans` package

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


## Analysing differences among factor levels

How different are site 2 and site 9?

```{r}
library("marginaleffects")
hypotheses(m3, "site2 = site9")
```



## Presenting model results

\footnotesize

```{r}
parameters(m3)
```


## Presenting model results

\scriptsize

```{r }
modelsummary(m3, estimate  = "{estimate} ({std.error})", statistic = NULL, 
             fmt = 1, gof_map = NA, coef_rename = paste0("site", 1:10), output = "markdown")
```



## Presenting model results

```{r }
library("gtsummary")
tbl_regression(m3)
```


## Plot

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



## Plot (visreg)

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

## Plot (easystats)

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


## Plot model (sjPlot)

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


## Plot model (modelsummary)

```{r out.width="70%"}
modelplot(m3)
```


## Plot model (easystats)

```{r out.width="70%"}
plot(parameters(m3), show_intercept = TRUE)
```


## Fit model without intercept

\scriptsize

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


## Model without intercept

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

## 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

```{r}
estimate_means(m4)
```


## Fit model without intercept

\scriptsize

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

## Plot

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



## Plot (visreg)

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

## Plot (visreg)

```{r}
visreg(m4, xvar = "dbh", by = "site", overlay = TRUE, band = FALSE)
```


## Plot model (sjPlot)

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

```

## Plot model (sjPlot)

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


## Plot model (easystats)

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


## Plot model (easystats)

Keeping sites only, dropping "dbh"

```{r}
plot(parameters(m4, drop = "dbh"))
```


## Plot model (modelsummary)

```{r}
modelplot(m4)
```


## Plot model (modelsummary)

Keeping sites only, dropping "dbh"

```{r}
modelplot(m4, coef_omit = "dbh")
```


## What happened to site 8?

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

::: {.column width="50%"}
```{r}
visreg(m3)
```
:::

::: {.column width="50%" }
```{r}
visreg(m4, xvar = "site")
```
:::
::::::::::::::





## What happened to site 8?

site 8 has the largest diameters

```{r echo = FALSE}
ggplot(trees) +
  geom_boxplot(aes(site, dbh))
```


## What happened to site 8?

\footnotesize

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

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

**DBH**

```{r}
aggregate(trees$dbh ~ trees$site, FUN = mean)
```
:::

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

```{r}
aggregate(trees$height ~ trees$site, FUN = mean)
```
:::
::::::::::::::






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

```{r echo=FALSE}
visreg(m4, xvar = "dbh", by = "site", overlay = TRUE, band = FALSE)
```


```{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, eval=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

\footnotesize

```{r}
parameters(m4, keep = "dbh")
```



## 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)
```


## How good is this model? Calibration plot (easystats)

```{r out.width="80%"}
pred <- estimate_expectation(m4)
pred$obs <- trees$height
plot(obs ~ Predicted, data = pred, xlim = c(15, 60), ylim = c(15, 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::check_predictions(m4)
```


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


# Predicting heights of new trees

## 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)
```


# Predicting heights of new trees (easystats)

## 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

Expected height of 10-cm DBH trees at each site

\footnotesize

```{r echo=1}
pred <- estimate_expectation(m4, data = trees.10cm)
pred
```

## Using model for prediction

Prediction intervals (accounting for residual variance)

\footnotesize

```{r echo=1}
pred <- estimate_prediction(m4, data = trees.10cm)
pred
```



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

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

```{r echo=FALSE}
df <- data.frame(dbh = seq(10, 50, by = 1), 
                 height = seq(20, 60, by = 1))
ggplot(df) +
  aes(dbh, height) +
  geom_blank() +
  geom_abline(intercept = 25, slope = 0.6) +
  geom_abline(intercept = 40, slope = 0.1, colour = "steelblue") +
  geom_abline(intercept = 50, slope = -0.3, colour = "orangered")
  
```



## 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")
```


## Does slope vary among sites?

\footnotesize

```{r}
visreg(m5, xvar = "dbh", by = "site", overlay = TRUE, band = FALSE)
```


## Does slope vary among sites?

\scriptsize

```{r}
library("marginaleffects")
hypotheses(m5, "`site9:dbh` = `site10:dbh`")
```


## Examining fitted model with {modelStudio}

```{r eval=FALSE}
library("modelStudio")
m5.explain <- DALEX::explain(
  m5, 
  data = trees, 
  y = trees$height)
modelStudio(m5.explain)
```


## 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](https://web.archive.org/web/20200201153653/http://blog.yhat.com/posts/7-funny-datasets.html): is speed related to sex?