Multilevel Generalized Linear Models

# Multilevel Generalized Linear Models
## cont’d
### 04.13.22

---

## [Click here for PDF of slides](26-multilevel-glms-pt2.pdf)

---

## Announcements

- Quiz 04 open due **Fri, April 15 at 11:59pm**

- Final project - optional draft due 
  - **Fri, Apr 15 at 11:59pm** 
  - final report due **Wed, Apr 27 at 11:59pm**

- Please fill out course evaluations!

- [Click here](https://sta310-sp22.github.io/slides/supplemental-notes/multilevel-model-questions.html) for answers to questions about multilevel models submitted on Quiz 03. Thanks to Jose for putting this document together!

---

## Learning goals

- Exploratory data analysis for multilevel data with non-normal response variable

- Write One, Level Two and composite models for mutlilevel GLM

- Fit and interpret multilevel GLM

- Understand crossed random effects and incorporate them in the multilevel model

---

## Data: College Basketball referees

The dataset [`basketball0910.csv`](data/basketball0910.csv) contains data on 4972 fouls in 340 NCAA basketball games from the Big Ten, ACC, and Big East conferences during the 2009-2010 season. The goal is to determine whether the data from this season support a conclusion from [Anderson and Pierce (2009)](https://www.tandfonline.com/doi/abs/10.1080/02640410902729733) that referees tend to "even out" foul calls in a game. The variables we'll focus on are

- **`foul.home`**: foul was called on home team (1: yes, 0: no)
- **`foul.diff`**: difference in fouls before current foul was called (home - visitor)
- **`game`**: Unique game ID number
- **`visitor`**: visiting team abbreviation 
- **`home`**: home team abbreviation

See [BMLR: Section 11.3.1](https://bookdown.org/roback/bookdown-BeyondMLR/ch-GLMM.html#explore-glmm) for full codebook.

---

## Data: College basketball referees

<table>
 <thead>
  <tr>
   <th style="text-align:right;"> game </th>
   <th style="text-align:left;"> visitor </th>
   <th style="text-align:left;"> hometeam </th>
   <th style="text-align:right;"> foul.num </th>
   <th style="text-align:right;"> foul.home </th>
   <th style="text-align:right;"> foul.vis </th>
   <th style="text-align:right;"> foul.diff </th>
   <th style="text-align:left;"> foul.type </th>
   <th style="text-align:right;"> time </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:left;"> IA </td>
   <td style="text-align:left;"> MN </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:left;"> Personal </td>
   <td style="text-align:right;"> 14.167 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:left;"> IA </td>
   <td style="text-align:left;"> MN </td>
   <td style="text-align:right;"> 2 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -1 </td>
   <td style="text-align:left;"> Personal </td>
   <td style="text-align:right;"> 11.433 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:left;"> IA </td>
   <td style="text-align:left;"> MN </td>
   <td style="text-align:right;"> 3 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:left;"> Personal </td>
   <td style="text-align:right;"> 10.233 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:left;"> IA </td>
   <td style="text-align:left;"> MN </td>
   <td style="text-align:right;"> 4 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:left;"> Personal </td>
   <td style="text-align:right;"> 9.733 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:left;"> IA </td>
   <td style="text-align:left;"> MN </td>
   <td style="text-align:right;"> 5 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:left;"> Shooting </td>
   <td style="text-align:right;"> 7.767 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:left;"> IA </td>
   <td style="text-align:left;"> MN </td>
   <td style="text-align:right;"> 6 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> -1 </td>
   <td style="text-align:left;"> Shooting </td>
   <td style="text-align:right;"> 5.567 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:left;"> IA </td>
   <td style="text-align:left;"> MN </td>
   <td style="text-align:right;"> 7 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -2 </td>
   <td style="text-align:left;"> Shooting </td>
   <td style="text-align:right;"> 2.433 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:left;"> IA </td>
   <td style="text-align:left;"> MN </td>
   <td style="text-align:right;"> 8 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -1 </td>
   <td style="text-align:left;"> Offensive </td>
   <td style="text-align:right;"> 1.000 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 2 </td>
   <td style="text-align:left;"> MI </td>
   <td style="text-align:left;"> MIST </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:left;"> Shooting </td>
   <td style="text-align:right;"> 18.983 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 2 </td>
   <td style="text-align:left;"> MI </td>
   <td style="text-align:left;"> MIST </td>
   <td style="text-align:right;"> 2 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -1 </td>
   <td style="text-align:left;"> Personal </td>
   <td style="text-align:right;"> 17.200 </td>
  </tr>
</tbody>
</table>

---

## Composite model

**Composite model**

.eq[
`$$\log\Big(\frac{p_{ij}}{1 - p_{ij}}\Big) = \alpha_0 + \beta_0 ~ \text{foul.diff}_{ij} + [u_i + v_i ~ \text{foul.diff}_{ij}]$$`
`$$\left[ \begin{array}{c}
            u_i \\ v_i
          \end{array}  \right] \sim N \left( \left[
          \begin{array}{c}
            0 \\ 0
          \end{array} \right], \left[
          \begin{array}{cc}
            \sigma_u^{2} & \sigma_{uv} \\
            \sigma_{uv} & \sigma_{v}^{2}
          \end{array} \right] \right)$$`
          
]

---

## Fit the model in R

Use the `glmer` function in the **lme4** package to fit multilevel GLMs.

```r
model1 <- glmer(foul.home ~ foul.diff + (foul.diff|game), 
                data = basketball, family = binomial)
```

```
## boundary (singular) fit: see help('isSingular')
```

<table>
 <thead>
  <tr>
   <th style="text-align:left;"> effect </th>
   <th style="text-align:left;"> group </th>
   <th style="text-align:left;"> term </th>
   <th style="text-align:right;"> estimate </th>
   <th style="text-align:right;"> std.error </th>
   <th style="text-align:right;"> statistic </th>
   <th style="text-align:right;"> p.value </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> fixed </td>
   <td style="text-align:left;"> NA </td>
   <td style="text-align:left;"> (Intercept) </td>
   <td style="text-align:right;"> -0.157 </td>
   <td style="text-align:right;"> 0.046 </td>
   <td style="text-align:right;"> -3.382 </td>
   <td style="text-align:right;"> 0.001 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> fixed </td>
   <td style="text-align:left;"> NA </td>
   <td style="text-align:left;"> foul.diff </td>
   <td style="text-align:right;"> -0.285 </td>
   <td style="text-align:right;"> 0.038 </td>
   <td style="text-align:right;"> -7.440 </td>
   <td style="text-align:right;"> 0.000 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> ran_pars </td>
   <td style="text-align:left;"> game </td>
   <td style="text-align:left;"> sd__(Intercept) </td>
   <td style="text-align:right;"> 0.542 </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
  </tr>
  <tr>
   <td style="text-align:left;"> ran_pars </td>
   <td style="text-align:left;"> game </td>
   <td style="text-align:left;"> cor__(Intercept).foul.diff </td>
   <td style="text-align:right;"> -1.000 </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
  </tr>
  <tr>
   <td style="text-align:left;"> ran_pars </td>
   <td style="text-align:left;"> game </td>
   <td style="text-align:left;"> sd__foul.diff </td>
   <td style="text-align:right;"> 0.035 </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
  </tr>
</tbody>
</table>

---

## Boundary constraints

- The estimates of the parameters `$\alpha_0, \beta_0, \sigma_u, \sigma_v, \rho_{uv}$` are those that maximize the likelihood of observing the data

- The fixed effects, e.g., `$\alpha_0$` and `$\beta_0$`, can take any values, but the terms associated with the error terms are constrained to a set of "allowable" values

`$$\sigma_u \geq 0 \hspace{10mm} \sigma_v \geq 0 \hspace{10mm} -1 \leq \rho_{uv} \leq 1$$`

- Because of these boundaries, a "constrained" search is used to find the MLEs.

- The warning message `"## boundary (singular) fit"`, means the estimate of one or more terms was set at the maximum (or minimum) allowable value, not the value it would've been if an unconstrained search were allowable

---

## Illustrating boundary constraints

Contour plots from a hypothetical likelihood `$L(\beta_0, \sigma^2)$`

<br>

<div class="figure" style="text-align: center">
<img src="data:image/png;base64,#26-multilevel-glms-pt2_files/figure-html/boundary-1.png" alt="From BMLR Figure 10.14" width="100%" />
<p class="caption">From BMLR Figure 10.14</p>
</div>
]

- In the unconstrained search, the likelihood `$L(\beta_0, \sigma^2)$` is maximized at `$\hat{\beta}_0 = 4, \hat{\sigma}^2 = -2$`

- In reality `$\sigma^2$` must be non-negative, so the search for the MLE is restricted to the region such that `$\sigma^2 \geq 0$`.

- The constrained likelihood is maximized at `$\hat{\beta}_0 = 4, \hat{\sigma}^2 = 0$`
]

---

**Original model**

.question[
1. Which term(s) should we remove to try to address boundary constraint? 
2. Refit the model with these terms removed. 
3. Which model do you choose? Use AIC to help make your choice. 
]

---

## Interpret model coefficients

Interpret the following coefficients in the selected model (if applicable):

- `$\hat{\alpha}_0$`

- `$\hat{\beta}_0$`

- `$\hat{\sigma}_u$`

- `$\hat{\sigma}_v$`

- `$\hat{\rho}_{uv}$`

---

## Crossed random effects

- The Level Two covariates are the home team and visiting team

- There is some evidence in the EDA that there may be differences in the probability of a foul depending on the home team

- We will account for this difference by treating home team and visiting team as random effects in the model

- *Issue*: Home and visiting team are not nested within game, since a single home and visiting team can be in multiple games

- The random effects for game, home team, and visiting team are **crossed random effects**

---

## Notation

`$Y_{i[gh]j}$`: Indicator of whether the `$j^{th}$` foul in Game `$i$` was called on home team `$h$` instead of visiting team `$g$`

`$$Y_{i[gh]j} \sim Bernoulli(p_{i[gh]j})$$`

where `$p_{i[gh]j}$` is the true probability a foul  in Game `$i$` was called on home team `$h$` instead of visiting team `$g$`

---

## Models

**Level One**

`$$\log\Big(\frac{p_{i[gh]j}}{1 - p_{i[gh]j}}\Big) = a_i + b_i ~ \text{foul.diff}_{ij}$$`

**Level Two**

`$$\begin{aligned}&a_i = \alpha_0 + u_i + v_h + w_g\\
&\beta_i = \beta_0\end{aligned}$$`

<br>

`$$u_i \sim N(0, \sigma^2_u) \hspace{10mm} v_h \sim N(0, \sigma^2_v) \hspace{10mm} w_g \sim N(0, \sigma^2_w)$$`

---

## Composite model

`$$\log\Big(\frac{p_{i[gh]j}}{1 - p_{i[gh]j}}\Big) = \alpha_0 + \beta_0 ~ \text{foul.diff}_{ij} + [u_i + v_h + w_g]$$`

`$$u_i \sim N(0, \sigma^2_u) \hspace{10mm} v_h \sim N(0, \sigma^2_v) \hspace{10mm} w_g \sim N(0, \sigma^2_w)$$`
]

**Why add additional random effects?**

- Get more precise estimates of fixed effects

- Can make comparisons of game-to-game and team-to-team variability

- Can get estimated random effects for each team and use them to compare odds of a foul on the home team for different teams

---

## Model

```r
model2 <- glmer(foul.home ~ foul.diff + (1|game) + (1|hometeam) + (1 | visitor),
               data = basketball, family = binomial)
```

```r
tidy(model2) %>% kable(digits = 3)
```

<table>
 <thead>
  <tr>
   <th style="text-align:left;"> effect </th>
   <th style="text-align:left;"> group </th>
   <th style="text-align:left;"> term </th>
   <th style="text-align:right;"> estimate </th>
   <th style="text-align:right;"> std.error </th>
   <th style="text-align:right;"> statistic </th>
   <th style="text-align:right;"> p.value </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> fixed </td>
   <td style="text-align:left;"> NA </td>
   <td style="text-align:left;"> (Intercept) </td>
   <td style="text-align:right;"> -0.188 </td>
   <td style="text-align:right;"> 0.063 </td>
   <td style="text-align:right;"> -2.967 </td>
   <td style="text-align:right;"> 0.003 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> fixed </td>
   <td style="text-align:left;"> NA </td>
   <td style="text-align:left;"> foul.diff </td>
   <td style="text-align:right;"> -0.264 </td>
   <td style="text-align:right;"> 0.039 </td>
   <td style="text-align:right;"> -6.795 </td>
   <td style="text-align:right;"> 0.000 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> ran_pars </td>
   <td style="text-align:left;"> game </td>
   <td style="text-align:left;"> sd__(Intercept) </td>
   <td style="text-align:right;"> 0.414 </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
  </tr>
  <tr>
   <td style="text-align:left;"> ran_pars </td>
   <td style="text-align:left;"> hometeam </td>
   <td style="text-align:left;"> sd__(Intercept) </td>
   <td style="text-align:right;"> 0.261 </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
  </tr>
  <tr>
   <td style="text-align:left;"> ran_pars </td>
   <td style="text-align:left;"> visitor </td>
   <td style="text-align:left;"> sd__(Intercept) </td>
   <td style="text-align:right;"> 0.152 </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
  </tr>
</tbody>
</table>

---

About what percent of the variability in the intercepts is due to...

- game-to-game differences?
- differences among home teams? 
- differences among visiting teams?

]

---

## Keep the crossed random effects?

- Given a large proportion of the variability in the intercepts is explained by game-to-game differences, we can assess if the random effects for home team and visiting team are providing useful information.

- To do so, we will compare the following models

```r
modela <-  glmer(foul.home ~ foul.diff + (1|game), 
                 data = basketball, family = binomial)

modelb <- glmer(foul.home ~ foul.diff + (1|game) + (1 | hometeam) + (1|visitor), 
                data = basketball, family = binomial)
```

---

## Keep the crossed random effects?

**Likelihood ratio test using** `${\color{#87037B}\chi^2}$` **distribution** (potentially unreliable when testing variance components)

```r
anova(modela, modelb, test = "Chisq") %>% kable(digits = 3)
```

<table>
 <thead>
  <tr>
   <th style="text-align:left;">   </th>
   <th style="text-align:right;"> npar </th>
   <th style="text-align:right;"> AIC </th>
   <th style="text-align:right;"> BIC </th>
   <th style="text-align:right;"> logLik </th>
   <th style="text-align:right;"> deviance </th>
   <th style="text-align:right;"> Chisq </th>
   <th style="text-align:right;"> Df </th>
   <th style="text-align:right;"> Pr(&gt;Chisq) </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> modela </td>
   <td style="text-align:right;"> 3 </td>
   <td style="text-align:right;"> 6792.540 </td>
   <td style="text-align:right;"> 6812.075 </td>
   <td style="text-align:right;"> -3393.270 </td>
   <td style="text-align:right;"> 6786.540 </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
  </tr>
  <tr>
   <td style="text-align:left;"> modelb </td>
   <td style="text-align:right;"> 5 </td>
   <td style="text-align:right;"> 6780.466 </td>
   <td style="text-align:right;"> 6813.024 </td>
   <td style="text-align:right;"> -3385.233 </td>
   <td style="text-align:right;"> 6770.466 </td>
   <td style="text-align:right;"> 16.074 </td>
   <td style="text-align:right;"> 2 </td>
   <td style="text-align:right;"> 0 </td>
  </tr>
</tbody>
</table>

<br>

**Parametric bootstrap** (long computational time!)

---

## Keep the crossed random effects?

**AIC or BIC**

```r
glance(modela) %>% kable(digits = 3)
```

<table>
 <thead>
  <tr>
   <th style="text-align:right;"> sigma </th>
   <th style="text-align:right;"> logLik </th>
   <th style="text-align:right;"> AIC </th>
   <th style="text-align:right;"> BIC </th>
   <th style="text-align:right;"> deviance </th>
   <th style="text-align:right;"> df.residual </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> -3393.27 </td>
   <td style="text-align:right;"> 6792.54 </td>
   <td style="text-align:right;"> 6812.075 </td>
   <td style="text-align:right;"> 6397.136 </td>
   <td style="text-align:right;"> 4969 </td>
  </tr>
</tbody>
</table>

```r
glance(modelb) %>% kable(digits = 3)
```

<table>
 <thead>
  <tr>
   <th style="text-align:right;"> sigma </th>
   <th style="text-align:right;"> logLik </th>
   <th style="text-align:right;"> AIC </th>
   <th style="text-align:right;"> BIC </th>
   <th style="text-align:right;"> deviance </th>
   <th style="text-align:right;"> df.residual </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> -3385.233 </td>
   <td style="text-align:right;"> 6780.466 </td>
   <td style="text-align:right;"> 6813.024 </td>
   <td style="text-align:right;"> 6420.537 </td>
   <td style="text-align:right;"> 4967 </td>
  </tr>
</tbody>
</table>

---

## Estimated random effects for each team

- We will use  `model B` with (crossed random effects) to get the estimated random effect for each team.

- These are **empirical Bayes estimates** ("shrinkage estimates"). 
  - Combine individual-specific information with information from all teams
  - "Shrinks" the individual estimates toward the group averages

<br>

.footnote[.small[See [On the Use of Empirical Bayes Estimates
as Measures of Individual Traits](https://sakai.duke.edu/access/content/group/ec265469-bdb3-47a7-beb0-317956b6f86f/empirical-bayes.pdf)] for more detail on empirical Bayes estimates.]

---

## Estimated random effects for each team

We can get these effects using the **`ranef`** function in the lmer R package.

```r
reffect_game <- ranef(modelb)$game %>% select(`(Intercept)`) %>% pull()
reffect_home <- ranef(modelb)$hometeam %>% select(`(Intercept)`) %>% pull()
reffect_visitor <- ranef(modelb)$visitor %>% select(`(Intercept)`) %>% pull()
team_names <- rownames(ranef(modelb)$visitor)
reffect_team <- tibble(team = team_names, 
                       reffect_home = reffect_home,
                     reffect_visitor = reffect_visitor)
```

```r
reffect_team %>%
  slice(1:5)
```

```
## # A tibble: 5 × 3
##   team   reffect_home reffect_visitor
##   <chr>         <dbl>           <dbl>
## 1 BC          -0.0900          0.0355
## 2 CIN          0.261          -0.0944
## 3 CLEM        -0.105          -0.110 
## 4 CT          -0.290           0.0360
## 5 DEPAUL       0.416          -0.154
```

]

---

## Distribution of random home team effects

```r
ggplot(data = reffect_team, aes(x = reffect_home)) + 
  geom_histogram(binwidth = .05, color = "black", fill = "steelblue") + 
  labs(x = "Random home team effects", 
       title = "Random home team effects from Model B") 
```

---

## Estimated home random effects by team

.panelset.sideways[
.panel[.panel-name[Plot]

<img src="data:image/png;base64,#26-multilevel-glms-pt2_files/figure-html/unnamed-chunk-15-1.png" width="100%" style="display: block; margin: auto;" />
]

```r
var <- attr(ranef(modelb)$hometeam, "postVar")
reffect_predict <- tibble(Intercepts = reffect_home,
              SD = 2*sqrt(var[,,1:length(var)]),
              team_names = reffect_team$team)

ggplot(data = reffect_predict, aes(fct_reorder(team_names, Intercepts), 
                                   Intercepts)) + 
  geom_point() + 
  geom_hline(yintercept = 0) + 
  geom_errorbar(aes(ymin = Intercepts - SD,
                    ymax = Intercepts + SD),
                width=0,color="black") + 
  labs(title = "Estimated random home team effects", 
       x = "Home teams", 
       y = "Estimated random effects") + 
 theme(axis.text.y = element_text(size = 7)) + 
  coord_flip()
```
]
]
]

---

## Modeling next steps

**Do the data provide evidence that referees tend to "even out" foul calls?**

- Adjust for additional covariates (score differential, type of foul, time left in first half)

- Hypothesize that the effect of foul differential may depend on other covariates, so consider potential interaction terms with foul differential

- Consider other random effects associated with game, home team, and visiting team

See [Section 11.7: A final model for examining referee bias](https://bookdown.org/roback/bookdown-BeyondMLR/ch-GLMM.html#sec:finalmodel-glmm) for final model chosen by the authors.

---

## Acknowledgements

- [BMLR: Section 10.5 - Encountering boundary constraints](https://bookdown.org/roback/bookdown-BeyondMLR/ch-3level.html#sec:boundary)

- [Chapter 11: Multilevel generalized linear models  structure](https://bookdown.org/roback/bookdown-BeyondMLR/ch-GLMM.html) 
  - Sections 11.1 - 11.4
  
- Anderson, Kyle J., and David A. Pierce. 2009. “Officiating Bias: The Effect of Foul Differential on Foul Calls in NCAA Basketball.” Journal of Sports Sciences 27 (7): 687–94. https://doi.org/10.1080/02640410902729733.