Generalized Linear Models

# Generalized Linear Models
## Wasser Cluster Lunz, Lunz am See
### Aurélien Boyé
### 08/10/2019

---

# Assumptions of Linear Models

1. **Linearity**
2. **Residuals**:
  - Independent
  - Equal variance
  - Normally distributed

1. If nonlinear terms are additive fit with OLS (polynomial regressions)
2. Transform? But think about what it will do to error 
3. Nonlinear Least Squares 
4. Generalized Linear Models & Mixed Models (GLM/GLMM)

---

# Assumptions of Linear Models

1. **Linearity**
2. **Residuals**:
  - Independent
  - Equal variance
  - Normally distributed

1. .comment[**If nonlinear terms are additive fit with OLS (polynomial regressions)**]
2. Transform? But think about what it will do to error 
3. Nonlinear Least Squares 
4. Generalized Linear Models & Mixed Models (GLM/GLMM)

---

# Polynomial regressions

&nbsp;

&nbsp;

<table>
 <thead>
  <tr>
   <th style="text-align:right;"> Degree </th>
   <th style="text-align:left;"> Name </th>
   <th style="text-align:left;"> Example </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:left;"> Constant </td>
   <td style="text-align:left;"> $3$ </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:left;"> Linear </td>
   <td style="text-align:left;"> $x+9$ </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 2 </td>
   <td style="text-align:left;"> Quadratic </td>
   <td style="text-align:left;"> $x^2-x+4$ </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 3 </td>
   <td style="text-align:left;"> Cubic </td>
   <td style="text-align:left;"> $x^3-x^2+5$ </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 4 </td>
   <td style="text-align:left;"> Quartic </td>
   <td style="text-align:left;"> $6x^4-x^3+x-2$ </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 5 </td>
   <td style="text-align:left;"> Quintic </td>
   <td style="text-align:left;"> $x^5-3x^3+x^2+8$ </td>
  </tr>
</tbody>
</table>

---

# Polynomial regressions

&nbsp;

```r
poly_reg <- lm(y ~ x + I(xˆ2) + I(xˆ3), data=data)
# or
poly_reg <- lm(y ~ poly(x, degree = 3), data=data)
```

&nbsp;

- However, note that `x`, `I(xˆ2)`, `I(xˆ3)` will be correlated if they are not centered variables ([Schielzeth, 2010](https://besjournals.onlinelibrary.wiley.com/doi/10.1111/j.2041-210X.2010.00012.x)), which can cause problems.

- The use of `poly()` avoid this by producing orthogonal polynomials.

- In general, if a higher order term is included (e.g. a quadratic term), the lower order terms (e.g. the intercept and linear terms) must also be in the model. This is known as the “principle of marginality”.

---

# Polynomial regressions

&nbsp;

]

---

# Assumptions of Linear Models

1. **Linearity**
2. **Residuals**:
  - Independent
  - Equal variance
  - Normally distributed

1. If nonlinear terms are additive fit with OLS (polynomial regressions)
2. .comment[**Transform? But think about what it will do to error**]
3. Nonlinear Least Squares 
4. Generalized Linear Models & Mixed Models (GLM/GLMM)

---

# Transformations

"When the error structure of data is simple, a transformation (usually a log or power-transformation) can be quite useful to improve the ability of a model to fit to the data by stabilizing variances or by making relationships linear before applying simple linear regression"

"However, a transformation is not guaranteed to solve these problems: there may be a trade-off between homoscedasticity and linearity, or the family of transformations used may not be able to correct one or both of these problems"

]

.footnote[
See the article [here](https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/j.2041-210X.2010.00021.x) and the blog post by Jarrett Byrnes [Do Not Log-Transform Count Data, Bitches!](https://www.r-bloggers.com/do-not-log-transform-count-data-bitches/)
]

---

# Transformations

**Common transformations**
- log(y)
- log(y) and log(x) (power function)
- arcsin(sqrt(y)) for bounded data
- logit for bounded dat `car::logit()`
- Box-Cox Transform `MASS::boxcox()`

**Illustration of the issue with error**
.center[
`$log(MetabolicRate) = log(a) + b \times log(mass) + error$`

implies

`$MetabolicRate = a \times mass^b \times e^{error}$`

but what we often want is

`$MetabolicRate = a \times mass^b + error$`
]

---

# Assumptions of Linear Models

1. **Linearity**
2. **Residuals**:
  - Independent
  - Equal variance
  - Normally distributed

1. If nonlinear terms are additive fit with OLS (polynomial regressions)
2. Transform? But think about what it will do to error
3. .comment[**Nonlinear Least Squares**] 
4. Generalized Linear Models & Mixed Models (GLM/GLMM)

---

# Nonlinear Least Squares fitting

```r
mod <- nls(metab_rate~a*mass^b, data = my_data,
           start = list(a=0.02, b=0.7))
```

- Least square method for fitting
- Very flexible
- .alert[Must specify start values]

---

# Assumptions of Linear Models

1. **Linearity**
2. **Residuals**:
  - Independent
  - Equal variance
  - Normally distributed

1. If nonlinear terms are additive fit with OLS (polynomial regressions)
2. Transform? But think about what it will do to error
3. Nonlinear Least Squares
4. Generalized Linear Models & Mixed Models (GLM/GLMM)

---

# Assumptions of Linear Models

1. **Linearity**
2. **Residuals**:
  - Independent
  - Equal variance
  - Normally distributed

]

1. If nonlinear terms are additive fit with OLS (polynomial regressions)
2. Transform? But think about what it will do to error
3. Nonlinear Least Squares
4. .alert[Generalized Linear Models & Mixed Models (GLM/GLMM)]

---

# Assumptions on the residuals

```r
mites <- read.csv('data/mites.csv')
head(mites,3)
#   Galumna pa totalabund        prop SubsDens WatrCont Substrate Shrub
# 1       8  1        140 0.057142857    39.18   350.15   Sphagn1   Few
# 2       3  1        268 0.011194030    54.99   434.81    Litter   Few
# 3       1  1        186 0.005376344    46.07   371.72 Interface   Few
#      Topo
# 1 Hummock
# 2 Hummock
# 3 Hummock
```

This dataset is a subset of the 'Oribatid mite dataset'

.small[
> `library(vegan); data(mite)` `$\Rightarrow$` 70 moss and mite samples
>
> 5 environmental measurements and abundance of the mite *Galumna sp.*
]

**Goal**: Model the abundance (`abund`), occurrence (`pa`), and proportion (`prop`) of Galumna as a function of the 5 environmental variables.

---

# Assumptions on the residuals

A negative relationship between `Galumna` and water `content`?

```r
  par(mfrow = c(1, 3), cex = 1.4)
  plot(Galumna ~ WatrCont, data = mites, xlab = 'Water content', ylab='Abundance')
  boxplot(WatrCont ~ pa, data = mites, xlab='Presence/Absence', ylab = 'Water content')
  plot(prop ~ WatrCont, data = mites, xlab = 'Water content', ylab='Proportion')
```

---

# Assumptions on the residuals

Fit linear models to test whether `abund`, `pa`, and/or `prop` varies as a function of water content.

```r
lm.abund <- lm(Galumna ~ WatrCont, data = mites)
lm.pa <- lm(pa ~ WatrCont, data = mites)
lm.prop <- lm(prop ~ WatrCont, data = mites)
```

--
.pull-left[

```r
# p-values
summary(lm.abund)$coefficients[, 4]
#  (Intercept)     WatrCont 
# 3.981563e-08 1.206117e-05
summary(lm.pa)$coefficients[, 4]
#  (Intercept)     WatrCont 
# 6.030252e-12 4.676755e-08
summary(lm.prop)$coefficients[, 4]
#  (Intercept)     WatrCont 
# 4.977432e-08 1.665437e-05
```
]
.pull-right[
&nbsp;

Significant relationship in all models!

---

# Assumptions on the residuals

Look at that error for the abundance !

```r
plot(Galumna ~ WatrCont, data = mites)
abline(lm.abund)
```

<img src="02_GLM_GLMM_files/figure-html/unnamed-chunk-12-1.png" style="display: block; margin: auto;" />
]

```r
par(mfrow = c(2, 2), cex = 1.4)
plot(lm.abund)
```

<img src="02_GLM_GLMM_files/figure-html/unnamed-chunk-13-1.png" style="display: block; margin: auto;" />
]

---

# Assumptions on the residuals

It's not better for the proportion `prop`:

```r
plot(prop ~ WatrCont, data = mites)
abline(lm.prop)
```

<img src="02_GLM_GLMM_files/figure-html/unnamed-chunk-14-1.png" style="display: block; margin: auto;" />
]

```r
par(mfrow = c(2, 2), cex = 1.4)
plot(lm.prop)
```

<img src="02_GLM_GLMM_files/figure-html/unnamed-chunk-15-1.png" style="display: block; margin: auto;" />
]

---

# Assumptions on the residuals

It's getting even worse for presence/absence `pa`:

```r
plot(pa ~ WatrCont, data = mites)
abline(lm.pa)
```

<img src="02_GLM_GLMM_files/figure-html/unnamed-chunk-16-1.png" style="display: block; margin: auto;" />
]

```r
par(mfrow = c(2, 2), cex = 1.4)
plot(lm.pa)
```

<img src="02_GLM_GLMM_files/figure-html/unnamed-chunk-17-1.png" style="display: block; margin: auto;" />
]

---

# Why be normal?

## *Your data is ok;*
## *it's the model that's wrong*

---

# Different ways to write the same model

#### Option 1

`$y_{i} = \beta_{0} + \beta_{1} x_{i} + \epsilon_{i}$`

`$\epsilon_{i} \sim N(0,\sigma^2)$`

]

#### Option 2

`$y_{i} \sim N(\mu_{i},\sigma^2)$`

`$\mu_{i} = \beta_{0} + \beta_{1} x_{i}$`

`$\epsilon_{i} \sim N(0,\sigma^2)$`

]

&nbsp;

This means that `$y_i$` is drawn from a normal distribution with parameters `$\mu$` (which depends on `$x_i$`) and `$\sigma$` (which has the same value for all `$Y$`s)

---

# Lets predict Galumna abundance

We need regression coefficients ( `$\beta$`) and `$\sigma$`:

```r
coef(lm.abund)
#  (Intercept)     WatrCont 
#  3.439348672 -0.006044788
summary(lm.abund)$sigma
# [1] 1.513531
```

What are the parameters of the normal distribution used to model `$y$` when water content = 300?

`$y_i \sim N(\mu = \beta_0 + \beta_1 X_i, \sigma^2)$`

`$\mu = 3.44 + (-0.006 \times 300) = 1.63$`

`$\sigma = 1.51$`

---

# Lets predict Galumna abundance

- At `$x = 300$`, residuals should follow a normal distribution with `$\mu = 1.63$` and `$\sigma^2 = 1.51$`.

- At `$x = 400$`, we get `$\mu = 1.02$` and `$\sigma^2 = 1.51$`, etc.

<br>
Graphically, this is our model:

--
.pull-left[
.center[
<img src="images/modelPredic.png" width="650" style="display: block; margin: auto;" />
]]
--
.pull-right[
**Problems**:
- `$\sigma^2$` is not homogeneous, yet `lm()` forces a constant `$\sigma^2$`
- Predicted values should be integers and positive
]

---

# Different ways to write the same model

#### Option 1

`$y_{i} = \beta_{0} + \beta_{1} x_{i} + \epsilon_{i}$`

`$\epsilon_{i} \sim N(0,\sigma^2)$`

]

#### Option 2

`$y_{i} \sim N(\mu_{i},\sigma^2)$`

`$\mu_{i} = \beta_{0} + \beta_{1} x_{i}$`

`$\epsilon_{i} \sim N(0,\sigma^2)$`

]

&nbsp;

This means that `$y_i$` is drawn from a normal distribution with parameters `$\mu$` (which depends on `$x_i$`) and `$\sigma$` (which has the same value for all `$Y$`s)

---

# Generalized linear models

Generalized linear models (GLMs) build on this approach. In general, they can be written as :

`$y_i \sim D(\mu_i,\phi)~~~~~$`   **(1)**

`$g(\mu_i) = \eta_i~~~~~~~~~~~$` **(2)**

`$\eta_i = \sum_{j=1}^s \beta_j x_{ij}~$`  **(3)**

]

where `$D()$` is a distribution, with mean `$\mu_i$` and the variance depending on `$\phi$`

- The error distribution `$D()$` is from the exponential family (*e.g.* Normal, Poisson, Binomial...) :

]

`$f(y|\theta,\phi) = exp \left( {\frac{y\theta - b(\theta)}{a(\phi)}+c(y, \phi)} \right)$`

]

&nbsp;

.small[
where `$a(·)$`, `$b(·)$` and `$c(·)$` are functions, `$\theta$` is the **natural parameter** (the parameter of interest) , and `$\phi$` is a scale parameter (known or seen as a nuisance)
]

---

# Some common distributions

**Normal distribution**

`$$f(x)={\frac {1}{\sigma {\sqrt {2\pi }}}}\operatorname {e} ^{-{\frac {1}{2}}\left({\frac {x-\mu }{\sigma }}\right)^{2}}$$`
]

Two parameters:

- `$\mu = E(x)$`
- `$\sigma^2 = Var(x)$`

---

# Some common distributions

**Normal distribution**

&nbsp;

---

# Some common distributions

**Poisson distribution**

&nbsp;

`$$P(X=k) = \frac{\lambda^k}{k!}e^{-\lambda}$$`

with a single parameter

`$$\lambda=E(X)=Var(X)$$`

---

# Some common distributions

**Poisson distribution**

---

# Some common distributions

**Bernouilli distribution**

`$$P(X=k)=p^{k}(1-p)^{1-k}$$` with `$$k\in \{0,1\}$$`

- Only two possible outcomes in its range: success (`1`) or failure (`0`)
- One parameter, `$p$`, the probability of success

---

# Some common distributions

**Bernouilli distribution**

&nbsp;

We can use the Bernouilli distribution to calculate the probability Galumna present (`1`) vs. absent (`0`)

---

# Some common distributions

**Binomial distribution**

`$$P(X=k)={n \choose k}\,p^{k}(1-p)^{n-k}$$`

When there are multiple trials (each with a success/failure), the Bernoulli distribution expands into the binomial
- Additional parameter, n, for number of trials
- Predicts the probability of observing a given proportion of successes, p, out of a known total number of trials, `$n$`

---

# Some common distributions

**Binomial distribution**

&nbsp;

---

# Generalized linear models

Generalized linear models (GLMs) build on this approach. In general, they can be written as :

`$y_i \sim D(\mu_i,\phi)~~~~~$`   **(1)**

`$g(\mu_i) = \eta_i~~~~~~~~~~~$` **(2)**

`$\eta_i = \sum_{j=1}^s \beta_j x_{ij}~$`  **(3)**

]

where `$D()$` is a distribution, with mean `$\mu_i$` and the variance depending on `$\phi$`

- The variance is a [function](https://en.wikipedia.org/wiki/Variance_function) of the fitted values `$var(Y_i) = \phi V(\mu_i)$` where the dispersion parameter `$\phi$` is a constant (that may be estimated from the data e.g. for the Normal distribution)
  - For the Normal distribution `$var(Y_i) = \phi \times 1$` as `$V(\mu_i)=1$` : it is constant with `$\phi = \sigma^2$` estimated from the data
  - For the Poisson distribution `$var(Y_i) =  1 \times \mu_i$` as `$V(\mu_i)=\mu_i$` : the variance increases with predicted values

]

---

# Generalized linear models

Generalized linear models (GLMs) build on this approach. In general, they can be written as :

`$y_i \sim D(\mu_i,\phi)~~~~~$`   **(1)**

`$g(\mu_i) = \eta_i~~~~~~~~~~~$` **(2)**

`$\eta_i = \sum_{j=1}^s \beta_j x_{ij}~$`  **(3)**

]

&nbsp;

- The model can be separated into a stochastic (1) and a deterministic (3) part, linked together by equation (2)

- Equation (3) is the linear part of the model: the covariates affect `$\eta_i$`, rather than `$\mu_i$`, in a linear way `$\Rightarrow$` `$\mu_i$` is seen as an invertible and smooth function of the linear predictor `$g(\mu_i) = \eta_i$` or `$\mu_i = g^{-1}(\eta_i)$`

- `$g()$` is called the **the link function**

---

# The link functions

&nbsp;

`$y_i \sim D(\mu_i,\phi)~~~~~$`   **(1)**

`$g(\mu_i) = \eta_i~~~~~~~~~~~$` **(2)**

`$\eta_i = \sum_{j=1}^s \beta_j x_{ij}~$`  **(3)**

]

&nbsp;

**Some common link functions `$g()$`**

- Identity: `$\mu = \eta~~~$` .comment[*e.g.*] `$~~~\mu = \alpha + \beta x$`
- Log: `$log(\mu) = \eta~~~$` .comment[*e.g.*] `$~~~\mu = e^{\alpha + \beta x}$`
- Logit: `$logit(\mu) = \eta~~~$` .comment[*e.g.*] `$~~~\mu = \frac{e^{\alpha + \beta x}}{1 +e^{\alpha + \beta x}}$`
- Inverse: `$\frac{1}{\mu} = \eta~~~$` .comment[*e.g.*] `$~~~\mu = (\alpha + \beta x)^{-1}$`

---

# The .comment[canonical] link functions

`$y_i \sim D(\mu_i,\phi)~~~~~$`   **(1)**

`$g(\mu_i) = \eta_i~~~~~~~~~~~$` **(2)**

`$\eta_i = \sum_{j=1}^s \beta_j x_{ij}~$`  **(3)**

]
]

with `$D()$` in the form of

&nbsp;

`$f(y|\theta,\phi) = exp \left( {\frac{y\theta - b(\theta)}{a(\phi)}+c(y, \phi)} \right)$`

]

&nbsp;

&nbsp;

`$g$` is called a **canonical link function** when

`$$\gamma'^{-1} \circ g^{-1} = (g\circ \gamma')^{-1} = I$$`

so that `$g$` connects `$\eta$`, `$\mu$`, and `$\theta$` such that `$\theta \equiv \eta$`

---

# Distributions, Canonical Links, and Dispersion

&nbsp;

<table>
 <thead>
  <tr>
   <th style="text-align:left;"> Distribution </th>
   <th style="text-align:left;"> Canonical Link </th>
   <th style="text-align:left;"> Variance Function </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> Normal </td>
   <td style="text-align:left;"> identity </td>
   <td style="text-align:left;"> ϕ </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Poisson </td>
   <td style="text-align:left;"> log </td>
   <td style="text-align:left;"> μ </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Quasipoisson </td>
   <td style="text-align:left;"> log </td>
   <td style="text-align:left;"> μϕ </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Binomial </td>
   <td style="text-align:left;"> logit </td>
   <td style="text-align:left;"> μ(1−μ) </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Quasibinomial </td>
   <td style="text-align:left;"> logit </td>
   <td style="text-align:left;"> μ(1 − μ)ϕ </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Negative Binomial </td>
   <td style="text-align:left;"> log </td>
   <td style="text-align:left;"> μ+kμ^2 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Gamma </td>
   <td style="text-align:left;"> inverse </td>
   <td style="text-align:left;"> μ2 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Inverse Normal </td>
   <td style="text-align:left;"> 1/μ2 </td>
   <td style="text-align:left;"> μ3 </td>
  </tr>
</tbody>
</table>

---

# Regressions are special cases of GLM

**Gneralized Linear Models**

`$y_i \sim D(\mu_i,\phi)~~~~~$`

`$g(\mu_i) = \eta_i~~~~~~~~~~~$`

`$\eta_i = \sum_{j=1}^s \beta_j x_{ij}~$`

]
]

**Normal distribution with identity link**

`$y_{i} \sim N(\mu_{i},\sigma^2)$`

`$\mu_{i} = \eta_{i}$`

`$\mu_{i} = \eta_{i} = \beta_{0} + \beta_{1} x_{i}$` 
]
]

---

# Common GLM

&nbsp;

---

# Distribution and other links

&nbsp;

<table>
 <thead>
  <tr>
   <th style="text-align:left;"> Distribution </th>
   <th style="text-align:left;"> Links </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> Normal </td>
   <td style="text-align:left;"> identity, log, inverse </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Poisson </td>
   <td style="text-align:left;"> log, identity, sqrt </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Quasipoisson </td>
   <td style="text-align:left;"> log, identity, sqrt </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Binomial </td>
   <td style="text-align:left;"> logit, probit, cauchit, log, log-log </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Quasibinomial </td>
   <td style="text-align:left;"> logit, probit, cauchit, log, log-log </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Negative Binomial </td>
   <td style="text-align:left;"> log, identity, sqrt </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Gamma </td>
   <td style="text-align:left;"> inverse, identity, log </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Inverse Normal </td>
   <td style="text-align:left;"> 1/μ2, inverse, identity, log </td>
  </tr>
</tbody>
</table>

---

---

# Exercise 1

&nbsp;

.large[
.center[
.comment[
Build a model of the presence of *Galumna* sp. as a function of water content and topography
]
]
]

&nbsp;

**Tip** :

- The `glm()` function
  - `glm_reg <- glm(formula, data, family)`

- Need to specify two things (`family`):
  1. probability distribution
  2. a link function

---

# Exercise 1: Solution

```r
mites <- read.csv("data/mites.csv", header = TRUE)
str(mites)
# 'data.frame':	70 obs. of  9 variables:
#  $ Galumna   : int  8 3 1 1 2 1 1 1 2 5 ...
#  $ pa        : int  1 1 1 1 1 1 1 1 1 1 ...
#  $ totalabund: int  140 268 186 286 199 209 162 126 123 166 ...
#  $ prop      : num  0.05714 0.01119 0.00538 0.0035 0.01005 ...
#  $ SubsDens  : num  39.2 55 46.1 48.2 23.6 ...
#  $ WatrCont  : num  350 435 372 360 204 ...
#  $ Substrate : Factor w/ 7 levels "Barepeat","Interface",..: 4 3 2 4 4 4 4 2 3 4 ...
#  $ Shrub     : Factor w/ 3 levels "Few","Many","None": 1 1 1 1 1 1 1 2 2 2 ...
#  $ Topo      : Factor w/ 2 levels "Blanket","Hummock": 2 2 2 2 2 2 2 1 1 2 ...
```

---

# Exercise 1: Solution

```r
logit.reg <- glm(pa ~ WatrCont + Topo, data=mites, family = binomial(link = "logit"))

summary(logit.reg)
# 
# Call:
# glm(formula = pa ~ WatrCont + Topo, family = binomial(link = "logit"), 
#     data = mites)
# 
# Deviance Residuals: 
#     Min       1Q   Median       3Q      Max  
# -2.0387  -0.5589  -0.1594   0.4112   2.0252  
# 
# Coefficients:
#              Estimate Std. Error z value Pr(>|z|)    
# (Intercept)  4.464402   1.670622   2.672 0.007533 ** 
# WatrCont    -0.015813   0.004535  -3.487 0.000489 ***
# TopoHummock  2.090757   0.735348   2.843 0.004466 ** 
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# (Dispersion parameter for binomial family taken to be 1)
# 
#     Null deviance: 91.246  on 69  degrees of freedom
# Residual deviance: 48.762  on 67  degrees of freedom
# AIC: 54.762
# 
# Number of Fisher Scoring iterations: 6
```

# Exercise 1: Solution

```r
par(mfrow=c(2,2))
plot(logit.reg)
```

---

```r
ggplot(aes(x = WatrCont, y = pa), data=mites) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = "glm", method.args = list(family = binomial(link = "logit"))) +
  labs(title = "Binomial logit") +
  theme_linedraw()
```

---

```r
ggplot(aes(x = WatrCont, y = pa), data=mites) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = "glm", method.args = list(family = binomial(link = "probit"))) +
  labs(title = "Binomial probit") +
  theme_linedraw()
```

---

```r
ggplot(aes(x = WatrCont, y = pa), data=mites) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = "glm", method.args = list(family = binomial(link = "cauchit"))) +
  labs(title = "Binomial cauchit") +
  theme_linedraw()
```

---

```r
ggplot(aes(x = WatrCont, y = pa), data=mites) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = "glm", method.args = list(family = binomial(link = "cloglog"))) +
  labs(title = "Binomial log-log") +
  theme_linedraw()
```

---

# Which one to use ?

&nbsp;

- **Coefficient of discrimination (D)** assess the predictive power of a logistic regression in terms of classification success or failure

- **Test de Hosmer-Lemeshow**: compare the number of observed *vs* expected values (similar to a `$Chi^2$` test)
  - Non significant p-value indicate adequate fit

---

# Coefficient D

```r
logit.reg <- glm(pa ~ WatrCont + Topo, data=mites, family = binomial(link = "logit"))

cauchit.reg <- glm(pa ~ WatrCont + Topo, data=mites, family = binomial(link = "cauchit"))

binomTools::Rsq(object = logit.reg)
# 
# R-square measures and the coefficient of discrimination, 'D':
# 
#      R2mod     R2res     R2cor         D
#  0.5205221 0.5024101 0.5025676 0.5114661
# 
# Number of binomial observations:  70
# Number of binary observation:  70
# Average group size:  1
binomTools::Rsq(object = cauchit.reg)
# 
# R-square measures and the coefficient of discrimination, 'D':
# 
#      R2mod     R2res     R2cor         D
#  0.5045138 0.5004369 0.5010976 0.5024754
# 
# Number of binomial observations:  70
# Number of binary observation:  70
# Average group size:  1
```
]

---

# Hosmer-Lemeshow test (logit)

```r
binomTools::HLtest(binomTools::Rsq(logit.reg))
# 
# Hosmer and Lemeshow's goodness-of-fit test for logistic regression models
# with method = deciles
# Observed counts:
#  1 2 3 4 5 6 7 8 9 10
#  7 7 7 6 5 6 3 2 2  0
#  0 0 0 1 2 1 4 5 5  7
# 
# Expected counts:
#  1   2   3   4   5   6   7   8 9  10
#  7 6.9 6.7 6.3 5.8 5.2 3.7 2.3 1 0.2
#  0 0.1 0.3 0.7 1.2 1.8 3.3 4.7 6 6.8
# 
# Pearson residuals:
#     1    2    3    4    5    6    7    8    9   10
#   0.0  0.0  0.1 -0.1 -0.3  0.4 -0.3 -0.2  1.0 -0.5
#  -0.1 -0.3 -0.6  0.3  0.7 -0.6  0.4  0.1 -0.4  0.1
# 
# Chi-square statistic:  3.421693  with  8  df
# P-value:  0.9051814 
# 
# Warning: Test may not be appropriate due to low expected frequencies
```
]

---

# Hosmer-Lemeshow test (cauchit)

```r
binomTools::HLtest(binomTools::Rsq(cauchit.reg))
# 
# Hosmer and Lemeshow's goodness-of-fit test for logistic regression models
# with method = deciles
# Observed counts:
#  1 2 3 4 5 6 7 8 9 10
#  7 7 7 6 5 6 3 2 2  0
#  0 0 0 1 2 1 4 5 5  7
# 
# Expected counts:
#    1   2   3   4   5   6   7   8   9  10
#  6.7 6.6 6.4 6.2 5.9 5.4 3.8 1.7 0.9 0.5
#  0.3 0.4 0.6 0.8 1.1 1.6 3.2 5.3 6.1 6.5
# 
# Pearson residuals:
#     1    2    3    4    5    6    7    8    9   10
#   0.1  0.2  0.2 -0.1 -0.4  0.3 -0.4  0.2  1.2 -0.7
#  -0.5 -0.6 -0.8  0.3  0.9 -0.5  0.4 -0.1 -0.5  0.2
# 
# Chi-square statistic:  5.231125  with  8  df
# P-value:  0.7326116 
# 
# Warning: Test may not be appropriate due to low expected frequencies
```
]

---

```r
ggplot(aes(x = WatrCont, y = pa), data=mites) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = "glm", method.args = list(family = poisson(link = "log"))) +
  labs(title = "Poisson log : Inappropriate distribution and link for this type data ") +
  theme_linedraw()
```

---

# Exercise 1: Interpreting the output

```r
summary(logit.reg)$coefficients
#                Estimate  Std. Error   z value     Pr(>|z|)
# (Intercept)  4.46440199 1.670622482  2.672299 0.0075333598
# WatrCont    -0.01581255 0.004535069 -3.486728 0.0004889684
# TopoHummock  2.09075654 0.735348234  2.843220 0.0044660283
```

- The output indicates that both water content and topography are significant

- **But how do we interpret the slope coefficients?**

---

# Exercise 1: Interpreting the output

&nbsp;

**Remember we used a logit transformation on the expected values!**

&nbsp;

To properly interpret the regression parameters, we have to use a **'reverse' function**

- The natural exponential function to obtain the **odds**: `$e^x$`
  - `$Odds = \frac{p}{1-p}$` with `$p$` the probability of success
  - Beta coefficient in logistic regressions correspond to the log of the odds ratio for a given predictor

- The inverse logit function to obtain the **probabilities**:

`$$logit^{-1} = \frac{1}{1 + \frac{1}{e^x}}$$`

---

# Exercise 1: Interpreting the output

- **On the odds scale**

```r
exp(logit.reg$coefficient[1])
# (Intercept) 
#    86.86907
exp(logit.reg$coefficient[2])
#  WatrCont 
# 0.9843118
```

> - The odds of the species being present when water content is 0 is 86.86 
- One unit increase in water content will multiply the odds of the species being present by 0.988 `$\Rightarrow$` decreases the odds

---

# Exercise 1: Interpreting the output

- **On the probability scale**

.small[You cannot just add the proba ! See [here](https://sebastiansauer.github.io/convert_logit2prob/) for further explanations]

```r
1 / (1 + 1/exp(logit.reg$coefficient[1]))
# (Intercept) 
#   0.9886194
1 / (1 + 1/exp(logit.reg$coefficient[1] + 1 * logit.reg$coefficient[2]))
# (Intercept) 
#   0.9884401
1 / (1 + 1/exp(logit.reg$coefficient[1] + 100 * logit.reg$coefficient[2]))
# (Intercept) 
#    0.947007
1 / (1 + 1/exp(logit.reg$coefficient[1] + 500 * logit.reg$coefficient[2]))
# (Intercept) 
#  0.03101209
```

>- The probability of presence when water content is 0 is 0.98
>- One unit increase in water content will slightly decrease the probability of the species occuring
>- When water content is 100, the proba of the species occuring is 0.947
]
---

# Exercise 1: Interpreting the output

**Note**:

- Coefficients are rarely interpreted quantitatively in GLM as in LM because it is not as straightforward but see [here](http://environmentalcomputing.net/interpreting-coefficients-in-glms/) for more examples with other distributions

---

# Predictive Power and goodness of fit

Get the pseudo-R², the analogue of the `$R^2$` for models fitted by maximum likelihood:

`$$\text{pseudo-R}^2 = \frac{\text{null deviance - residual deviance}}{\text{null deviance}}$$`

<br>

`$\text{pseudo-R}^2 = \text{variance explained by the model}$`

---

# Predictive Power and goodness of fit

Comparing deviance of your model (residual deviance) to the deviance of a null model (null deviance)

The **null model** is a model without explanatory variables

```R
null.model <- glm(Response.variable ~ 1, family = binomial)
```

---

# Predictive Power and goodness of fit

In R, we can extract the residual and null deviances directly from the glm object:

```r
objects(logit.reg)
#  [1] "aic"               "boundary"          "call"             
#  [4] "coefficients"      "contrasts"         "control"          
#  [7] "converged"         "data"              "deviance"         
# [10] "df.null"           "df.residual"       "effects"          
# [13] "family"            "fitted.values"     "formula"          
# [16] "iter"              "linear.predictors" "method"           
# [19] "model"             "null.deviance"     "offset"           
# [22] "prior.weights"     "qr"                "R"                
# [25] "rank"              "residuals"         "terms"            
# [28] "weights"           "xlevels"           "y"
```

```r
pseudoR2 <- (logit.reg$null.deviance - logit.reg$deviance) / logit.reg$null.deviance
pseudoR2
# [1] 0.4655937
```

---

# Predictive Power and goodness of fit

```r
# For binomial model only; based on Tjur (2009)
# devtools::install_version("binomTools", version = "1.0-1")
binomTools::Rsq(object = logit.reg) 
# 
# R-square measures and the coefficient of discrimination, 'D':
# 
#      R2mod     R2res     R2cor         D
#  0.5205221 0.5024101 0.5025676 0.5114661
# 
# Number of binomial observations:  70
# Number of binary observation:  70
# Average group size:  1

# Based on Nakagawa and Schielzeth (2013) & Jaeger et al. (2016)
r2glmm::r2beta(logit.reg, partial = FALSE, method = 'nsj')
#   Effect   Rsq upper.CL lower.CL
# 1  Model 0.478    0.625     0.33
r2glmm::r2beta(logit.reg, partial = TRUE, method = 'nsj')
#        Effect   Rsq upper.CL lower.CL
# 1       Model 0.478    0.625    0.330
# 2    WatrCont 0.364    0.528    0.203
# 3 TopoHummock 0.044    0.179    0.000
```

---

# Predictive Power and goodness of fit

```r
# Revised statistics based on Nakagawa et al. (2017) paper
MuMIn::r.squaredGLMM(logit.reg)
#                   R2m       R2c
# theoretical 0.7093845 0.7093845
# delta       0.6483495 0.6483495
```

---

# Exercise 2: count data

&nbsp;

.large[
.center[
.comment[
Build a model of the **abundance** of *Galumna* sp. as a function of water content and topography
]
]
]

---

# Exercise 2: Solution

```r
glm.poisson <- glm(Galumna ~ WatrCont + Topo, data=mites,
family = poisson(link = "log"))
```

---

```r
summary(glm.poisson)
# 
# Call:
# glm(formula = Galumna ~ WatrCont + Topo, family = poisson(link = "log"), 
#     data = mites)
# 
# Deviance Residuals: 
#     Min       1Q   Median       3Q      Max  
# -1.8977  -0.9760  -0.6431  -0.1651   4.0743  
# 
# Coefficients:
#              Estimate Std. Error z value Pr(>|z|)    
# (Intercept)  1.937184   0.449489   4.310 1.63e-05 ***
# WatrCont    -0.006736   0.001146  -5.877 4.17e-09 ***
# TopoHummock  0.615095   0.280823   2.190   0.0285 *  
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# (Dispersion parameter for poisson family taken to be 1)
# 
#     Null deviance: 168.25  on 69  degrees of freedom
# Residual deviance: 102.98  on 67  degrees of freedom
# AIC: 176.14
# 
# Number of Fisher Scoring iterations: 6
```

---

# Exercise 2: Solution

```r
ggplot(aes(x = WatrCont, y = Galumna), data=mites) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = "glm", method.args = list(family = poisson(link = "log"))) +
  labs(title = "Poisson log") +
  theme_linedraw()
```

---

# Overdispersion in count data

- The residual deviance is defined as twice the difference between the log likelihood of a model that provides a perfect fit and the log likelihood of our model

`$$\text{Res dev} = 2 \, log(L(y;\,y)) - 2 \, log(L(y;\, \mu))$$`

- In a Poisson GLM, the residual deviance should equal the residual degrees of freedoms

- When the residual deviance is higher than the residual degrees of freedom we say that the model is **overdispersed**

- Occurs when the variance in the data is even higher than the mean, hence the Poisson distribution is not the best choice (many zeros, many very high values, missing covariates, etc)

---

# Use other distributions

&nbsp;

<table>
 <thead>
  <tr>
   <th style="text-align:left;"> Distribution </th>
   <th style="text-align:left;"> Canonical Link </th>
   <th style="text-align:left;"> Variance Function </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> Normal </td>
   <td style="text-align:left;"> identity </td>
   <td style="text-align:left;"> θ </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Poisson </td>
   <td style="text-align:left;"> log </td>
   <td style="text-align:left;"> μ </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Quasipoisson </td>
   <td style="text-align:left;"> log </td>
   <td style="text-align:left;"> μθ </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Binomial </td>
   <td style="text-align:left;"> logit </td>
   <td style="text-align:left;"> μ(1−μ) </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Quasibinomial </td>
   <td style="text-align:left;"> logit </td>
   <td style="text-align:left;"> μ(1 − μ)θ </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Negative Binomial </td>
   <td style="text-align:left;"> log </td>
   <td style="text-align:left;"> μ+kμ^2 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Gamma </td>
   <td style="text-align:left;"> inverse </td>
   <td style="text-align:left;"> μ2 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Inverse Normal </td>
   <td style="text-align:left;"> 1/μ2 </td>
   <td style="text-align:left;"> μ3 </td>
  </tr>
</tbody>
</table>

---

# Use other distributions

&nbsp;

<table>
 <thead>
  <tr>
   <th style="text-align:left;"> Distribution </th>
   <th style="text-align:left;"> Canonical Link </th>
   <th style="text-align:left;"> Variance Function </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> Normal </td>
   <td style="text-align:left;"> identity </td>
   <td style="text-align:left;"> θ </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Poisson </td>
   <td style="text-align:left;"> log </td>
   <td style="text-align:left;"> μ </td>
  </tr>
  <tr>
   <td style="text-align:left;font-weight: bold;color: white !important;background-color: #D7261E !important;"> Quasipoisson </td>
   <td style="text-align:left;font-weight: bold;color: white !important;background-color: #D7261E !important;"> log </td>
   <td style="text-align:left;font-weight: bold;color: white !important;background-color: #D7261E !important;"> μθ </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Binomial </td>
   <td style="text-align:left;"> logit </td>
   <td style="text-align:left;"> μ(1−μ) </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Quasibinomial </td>
   <td style="text-align:left;"> logit </td>
   <td style="text-align:left;"> μ(1 − μ)θ </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Negative Binomial </td>
   <td style="text-align:left;"> log </td>
   <td style="text-align:left;"> μ+kμ^2 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Gamma </td>
   <td style="text-align:left;"> inverse </td>
   <td style="text-align:left;"> μ2 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Inverse Normal </td>
   <td style="text-align:left;"> 1/μ2 </td>
   <td style="text-align:left;"> μ3 </td>
  </tr>
</tbody>
</table>

---

# Quasi-Poisson GLM

The **systematic part** and the **link function** remain the same as in a Poisson GLM but the variance changes with

`$$Var(Y_I) = \phi \times \mu_i$$`

instead of

`$$Var(Y_I) = \mu_i$$`

`$\phi$` is the dispersion parameter. It will be estimated prior to estimate parameters. Correcting for overdispersion will not affect parameter estimates but will affect their **significance**. Indeed, the standard errors of the parameters are multiplied by `$\sqrt{\phi}$`

&nbsp;

---

# Quasi-Poisson GLM

Create a new GLM using the 'quasipoisson' family or update the previous one

```r
glm.quasipoisson <- glm(Galumna ~ WatrCont + Topo, data=mites,
family = quasipoisson)

glm.quasipoisson <- update(glm.poisson,family = quasipoisson)
```

When comparing the outputs with a poisson model see:

- changes in the *se* of the coefficients
- the AIC `$\Rightarrow$` The quasi-Poisson is not a full maximum likelihood (ML) model but a quasi-ML model, so that the AIC is not defined

---

```r
summary(glm.quasipoisson)
# 
# Call:
# glm(formula = Galumna ~ WatrCont + Topo, family = quasipoisson, 
#     data = mites)
# 
# Deviance Residuals: 
#     Min       1Q   Median       3Q      Max  
# -1.8977  -0.9760  -0.6431  -0.1651   4.0743  
# 
# Coefficients:
#              Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  1.937184   0.666256   2.908 0.004935 ** 
# WatrCont    -0.006736   0.001699  -3.965 0.000181 ***
# TopoHummock  0.615095   0.416250   1.478 0.144172    
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# (Dispersion parameter for quasipoisson family taken to be 2.197067)
# 
#     Null deviance: 168.25  on 69  degrees of freedom
# Residual deviance: 102.98  on 67  degrees of freedom
# AIC: NA
# 
# Number of Fisher Scoring iterations: 6
```

---

```r
ggplot(aes(x = WatrCont, y = Galumna), 
    data=mites) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = "glm", 
    method.args = list(
      family = "quasipoisson")) +
  labs(title = "Quasipoisson Model") +
  theme_linedraw()
```

]

```r
ggplot(aes(x = WatrCont, y = Galumna), 
    data=mites) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = "glm", 
    method.args = list(
      family = "poisson")) +
  labs(title = "Poisson Model") +
  theme_linedraw()
```

]

---

# Degree of overdispersion

The degree of overdispersion can be computed as

`$$\phi \, = \, \frac{\text{Residual deviance}}{\text{Residual degrees of freedom}}$$`

> In our case : `$102.98/67 = 1.537$`

Depending on the degree of overdispersion, solutions may be :

&nbsp;

1. Correct for overdispersion using **quasi-Poisson GLM**
2. Choose another distribution: **the negative binomial**
]

---

# Negative binomial GLM

Negative binomial GLMs are favor when overdispersion is high

- It has **two parameters** `$\mu$` and `$k$`. `$k$` controls for the dispersion parameter (smaller `$k$` indicates higher dispersion)
- It corresponds to a combination of two distributions (**Poisson** and **gamma**)
- It assumes that the `$Y_i$` are Poisson distributed with the mean `$\mu$` assumed to follow a gamma distribution

`$$E(Y_i) = \mu_i$$`

`$$Var(Y_i) = \mu_i + \frac{\mu^2_i}{k}$$`

---

# Fitting a negative binomial in R

Negative Binomial is not in the `glm()` function so you need the `MASS` package

```r
install.packages('MASS')
```

```r
glm.negbin <- MASS::glm.nb(Galumna ~ WatrCont + Topo, data=mites)
```

```r
summary(glm.negbin)
```

---

```
# 
# Call:
# MASS::glm.nb(formula = Galumna ~ WatrCont + Topo, data = mites, 
#     init.theta = 0.9528403246, link = log)
# 
# Deviance Residuals: 
#     Min       1Q   Median       3Q      Max  
# -1.4964  -0.8113  -0.5068  -0.1635   2.5071  
# 
# Coefficients:
#              Estimate Std. Error z value Pr(>|z|)    
# (Intercept)  2.346171   0.746735   3.142  0.00168 ** 
# WatrCont    -0.008303   0.001863  -4.456 8.34e-06 ***
# TopoHummock  0.832092   0.423744   1.964  0.04957 *  
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# (Dispersion parameter for Negative Binomial(0.9528) family taken to be 1)
# 
#     Null deviance: 88.685  on 69  degrees of freedom
# Residual deviance: 51.295  on 67  degrees of freedom
# AIC: 159.76
# 
# Number of Fisher Scoring iterations: 1
# 
# 
#               Theta:  0.953 
#           Std. Err.:  0.413 
# 
#  2 x log-likelihood:  -151.763
```
]]
.pull-right2[

<br><br><br><br><br><br><br><br><br><br><br><br><br>
`theta` `$= k$`
]

---

# In summary

---

<table>
 <thead>
  <tr>
   <th style="text-align:left;"> Names </th>
   <th style="text-align:left;"> Distribution </th>
   <th style="text-align:left;"> Link function </th>
   <th style="text-align:left;"> Applications </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> Regression, ANOVA </td>
   <td style="text-align:left;"> Normal </td>
   <td style="text-align:left;"> Identity </td>
   <td style="text-align:left;"> Normal homoscedastic unbounded data </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Logistic Regression </td>
   <td style="text-align:left;"> Binomial </td>
   <td style="text-align:left;"> logit </td>
   <td style="text-align:left;"> Presence/absence data </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Logistic Regression </td>
   <td style="text-align:left;"> Binomial </td>
   <td style="text-align:left;"> log-log </td>
   <td style="text-align:left;"> Presence/absence dat with many 0s or many 1s </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Probit Analysis </td>
   <td style="text-align:left;"> Binomial </td>
   <td style="text-align:left;"> probit </td>
   <td style="text-align:left;"> Proportions </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Probit Analysis </td>
   <td style="text-align:left;"> Quasininomial </td>
   <td style="text-align:left;"> probit </td>
   <td style="text-align:left;"> Proportions with overdispersion </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Log-linear model </td>
   <td style="text-align:left;"> Poisson </td>
   <td style="text-align:left;"> log </td>
   <td style="text-align:left;"> Count data </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Log-linear model </td>
   <td style="text-align:left;"> Quasipoisson </td>
   <td style="text-align:left;"> log </td>
   <td style="text-align:left;"> Count data with overdispersion </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Negative binomial </td>
   <td style="text-align:left;"> Negative binomial </td>
   <td style="text-align:left;"> log </td>
   <td style="text-align:left;"> Count data with strong overdispersion </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Ordinal Regression </td>
   <td style="text-align:left;"> Multinomial </td>
   <td style="text-align:left;"> logit </td>
   <td style="text-align:left;"> Multiclass data (factors with multiple states) </td>
  </tr>
</tbody>
</table>

]

---

# GLMM

---

# Generalized Linear Mixed Models

Extension of GLMs to account for additional structure in the data

Follows similar steps introduced for Linear Mixed Models (LMM)

1. LMMs incorporate random effects
2. GLMs handle non-normal data (letting errors take on different distribution families)

```r
library(lme4)

mod <- glmer(y ~ x1*x2 + x3 + (1|site/quadrat), data = df, family ="poison")
```

.footnote[
Too slow ? Try different [optimizers](http://svmiller.com/blog/2018/06/mixed-effects-models-optimizer-checks/)
]

---

## Ressources

&nbsp;

]