Agenda

• More equation practice

• Finishing up our intro to Bayes

  • Specific focus on MCMC

• Implementation with {brms}

• Logistic regression review

• Extending to multilevel logistic regression models with lme4

The data

Read in the following data

library(tidyverse)
library(lme4)
d <- read_csv(here::here("data", "three-lev.csv"))
dCopy Code

## # A tibble: 7,230 x 12
##    schid       sid  size lowinc mobility female black hispanic retained grade  year   math
##    <dbl>     <dbl> <dbl>  <dbl>    <dbl>  <dbl> <dbl>    <dbl>    <dbl> <dbl> <dbl>  <dbl>
##  1  2020 273026452   380   40.3     12.5      0     0        1        0     2   0.5  1.146
##  2  2020 273026452   380   40.3     12.5      0     0        1        0     3   1.5  1.134
##  3  2020 273026452   380   40.3     12.5      0     0        1        0     4   2.5  2.3  
##  4  2020 273030991   380   40.3     12.5      0     0        0        0     0  -1.5 -1.303
##  5  2020 273030991   380   40.3     12.5      0     0        0        0     1  -0.5  0.439
##  6  2020 273030991   380   40.3     12.5      0     0        0        0     2   0.5  2.43 
##  7  2020 273030991   380   40.3     12.5      0     0        0        0     3   1.5  2.254
##  8  2020 273030991   380   40.3     12.5      0     0        0        0     4   2.5  3.873
##  9  2020 273059461   380   40.3     12.5      0     0        1        0     0  -1.5 -1.384
## 10  2020 273059461   380   40.3     12.5      0     0        1        0     1  -0.5  0.338
## # … with 7,220 more rowsCopy Code

Model 1

Fit the following model

$$\begin{array}{cc}{\mathrm{math}}_i& \sim N({\alpha }_{j\left[i\right]},{\sigma }^2)\\ {\alpha }_j& \sim N({\gamma }_0^{\alpha }+{\gamma }_1^{\alpha }(\mathrm{mobility}),{\sigma }_{{\alpha }_j}^2)\text{, for schid j = 1,}\dots \text{,J}\end{array}$$

lmer(math ~ mobility + (1|schid), data = d)Copy Code

Model 2

Fit the following model

$$\begin{array}{cc}{\mathrm{math}}_i& \sim N({\alpha }_{j\left[i\right],k\left[i\right]}+{\beta }_{1j\left[i\right]}(\mathrm{year}),{\sigma }^2)\\ \left(\begin{array}{c}\begin{array}{cc}& {\alpha }_j\\ & {\beta }_{1j}\end{array}\end{array}\right)& \sim N(\left(\begin{array}{c}\begin{array}{cc}& {\gamma }_0^{\alpha }+{\gamma }_1^{\alpha }\left(\mathrm{female}\right)\\ & {\mu }_{{\beta }_{1j}}\end{array}\end{array}\right),\left(\begin{array}{cc}{\sigma }_{{\alpha }_j}^2& {\rho }_{{\alpha }_j{\beta }_{1j}}\\ {\rho }_{{\beta }_{1j}{\alpha }_j}& {\sigma }_{{\beta }_{1j}}^2\end{array}\right))\text{, for sid j = 1,}\dots \text{,J}\\ {\alpha }_k& \sim N({\gamma }_0^{\alpha }+{\gamma }_1^{\alpha }(\mathrm{mobility}),{\sigma }_{{\alpha }_k}^2)\text{, for schid k = 1,}\dots \text{,K}\end{array}$$

lmer(math ~ year + female + mobility + 
             (year|sid) +  (1|schid), 
     data = d)Copy Code

Model 3

Fit the following model

$$\begin{array}{cc}{\mathrm{math}}_i& \sim N({\alpha }_{j\left[i\right],k\left[i\right]}+{\beta }_{1j\left[i\right]}(\mathrm{year}),{\sigma }^2)\\ \left(\begin{array}{c}\begin{array}{cc}& {\alpha }_j\\ & {\beta }_{1j}\end{array}\end{array}\right)& \sim N(\left(\begin{array}{c}\begin{array}{cc}& {\gamma }_0^{\alpha }+{\gamma }_1^{\alpha }\left(\mathrm{black}\right)+{\gamma }_2^{\alpha }\left(\mathrm{hispanic}\right)+{\gamma }_3^{\alpha }\left(\mathrm{female}\right)\\ & {\gamma }_0^{{\beta }_1}+{\gamma }_1^{{\beta }_1}\left(\mathrm{black}\right)+{\gamma }_2^{{\beta }_1}\left(\mathrm{hispanic}\right)\end{array}\end{array}\right),\left(\begin{array}{cc}{\sigma }_{{\alpha }_j}^2& {\rho }_{{\alpha }_j{\beta }_{1j}}\\ {\rho }_{{\beta }_{1j}{\alpha }_j}& {\sigma }_{{\beta }_{1j}}^2\end{array}\right))\text{, for sid j = 1,}\dots \text{,J}\\ {\alpha }_k& \sim N({\gamma }_0^{\alpha }+{\gamma }_1^{\alpha }(\mathrm{mobility}),{\sigma }_{{\alpha }_k}^2)\text{, for schid k = 1,}\dots \text{,K}\end{array}$$

lmer(math ~ year * black + year * hispanic + female + mobility + 
             (year|sid) +  (1|schid), 
     data = d)Copy Code

Model 4

Fit the following model

$$\begin{array}{cc}{\mathrm{math}}_i& \sim N({\alpha }_{j\left[i\right],k\left[i\right]}+{\beta }_{1j\left[i\right],k\left[i\right]}(\mathrm{year}),{\sigma }^2)\\ \left(\begin{array}{c}\begin{array}{cc}& {\alpha }_j\\ & {\beta }_{1j}\end{array}\end{array}\right)& \sim N(\left(\begin{array}{c}\begin{array}{cc}& {\gamma }_0^{\alpha }+{\gamma }_{1k\left[i\right]}^{\alpha }\left(\mathrm{female}\right)\\ & {\gamma }_0^{{\beta }_1}+{\gamma }_1^{{\beta }_1}\left(\mathrm{female}\right)\end{array}\end{array}\right),\left(\begin{array}{cc}{\sigma }_{{\alpha }_j}^2& {\rho }_{{\alpha }_j{\beta }_{1j}}\\ {\rho }_{{\beta }_{1j}{\alpha }_j}& {\sigma }_{{\beta }_{1j}}^2\end{array}\right))\text{, for sid j = 1,}\dots \text{,J}\\ \left(\begin{array}{c}\begin{array}{cc}& {\alpha }_k\\ & {\beta }_{1k}\\ & {\gamma }_{1k}\end{array}\end{array}\right)& \sim N(\left(\begin{array}{c}\begin{array}{cc}& {\gamma }_0^{\alpha }+{\gamma }_1^{\alpha }\left(\mathrm{mobility}\right)+{\gamma }_2^{\alpha }\left(\mathrm{lowinc}\right)\\ & {\gamma }_0^{{\beta }_1}+{\gamma }_1^{{\beta }_1}\left(\mathrm{mobility}\right)\\ & {\mu }_{{\gamma }_{1k}}\end{array}\end{array}\right),\left(\begin{array}{ccc}{\sigma }_{{\alpha }_k}^2& {\rho }_{{\alpha }_k{\beta }_{1k}}& {\rho }_{{\alpha }_k{\gamma }_{1k}}\\ {\rho }_{{\beta }_{1k}{\alpha }_k}& {\sigma }_{{\beta }_{1k}}^2& {\rho }_{{\beta }_{1k}{\gamma }_{1k}}\\ {\rho }_{{\gamma }_{1k}{\alpha }_k}& {\rho }_{{\gamma }_{1k}{\beta }_{1k}}& {\sigma }_{{\gamma }_{1k}}^2\end{array}\right))\text{, for schid k = 1,}\dots \text{,K}\end{array}$$

lmer(math ~ year * female + year * mobility + lowinc +
             (year|sid) +  (year + female|schid), 
     data = d)Copy Code

Model 5

Fit the following model

Don't worry if you run into convergene warnings

$$\begin{array}{cc}{\mathrm{math}}_i& \sim N({\alpha }_{j\left[i\right],k\left[i\right]}+{\beta }_{1j\left[i\right],k\left[i\right]}(\mathrm{year}),{\sigma }^2)\\ \left(\begin{array}{c}\begin{array}{cc}& {\alpha }_j\\ & {\beta }_{1j}\end{array}\end{array}\right)& \sim N(\left(\begin{array}{c}\begin{array}{cc}& {\gamma }_0^{\alpha }+{\gamma }_{1k\left[i\right]}^{\alpha }\left(\mathrm{female}\right)\\ & {\gamma }_{1k[i0}^{{\beta }_1}+{\gamma }_{1k\left[i\right]}^{{\beta }_1}\left(\mathrm{female}\right)\end{array}\end{array}\right),\left(\begin{array}{cc}{\sigma }_{{\alpha }_j}^2& {\rho }_{{\alpha }_j{\beta }_{1j}}\\ {\rho }_{{\beta }_{1j}{\alpha }_j}& {\sigma }_{{\beta }_{1j}}^2\end{array}\right))\text{, for sid j = 1,}\dots \text{,J}\\ \left(\begin{array}{c}\begin{array}{cc}& {\alpha }_k\\ & {\beta }_{1k}\\ & {\gamma }_{1k}\\ & {\gamma }_{1k}^{{\beta }_1}\end{array}\end{array}\right)& \sim N(\left(\begin{array}{c}\begin{array}{cc}& {\gamma }_0^{\alpha }+{\gamma }_1^{\alpha }\left(\mathrm{mobility}\right)+{\gamma }_2^{\alpha }\left(\mathrm{lowinc}\right)\\ & {\mu }_{{\beta }_{1k}}\\ & {\mu }_{{\gamma }_{1k}}\\ & {\mu }_{{\gamma }_{1k}^{{\beta }_1}}\end{array}\end{array}\right),\left(\begin{array}{cccc}{\sigma }_{{\alpha }_k}^2& 0& 0& 0\\ 0& {\sigma }_{{\beta }_{1k}}^2& 0& 0\\ 0& 0& {\sigma }_{{\gamma }_{1k}}^2& 0\\ 0& 0& 0& {\sigma }_{{\gamma }_{1k}^{{\beta }_1}}^2\end{array}\right))\text{, for schid k = 1,}\dots \text{,K}\end{array}$$

lmer(math ~ year * female + mobility + lowinc +
       (year|sid) +  (year * female||schid), 
     data = d)Copy Code

Model 6

Fit the following model

$$\begin{array}{cc}{\mathrm{math}}_i& \sim N({\alpha }_{j\left[i\right]}+{\beta }_{1j\left[i\right]}(\mathrm{year}),{\sigma }^2)\\ \left(\begin{array}{c}\begin{array}{cc}& {\alpha }_j\\ & {\beta }_{1j}\end{array}\end{array}\right)& \sim N(\left(\begin{array}{c}\begin{array}{cc}& {\gamma }_0^{\alpha }+{\gamma }_1^{\alpha }\left(\mathrm{female}\right)\\ & {\mu }_{{\beta }_{1j}}\end{array}\end{array}\right),\left(\begin{array}{cc}{\sigma }_{{\alpha }_j}^2& {\rho }_{{\alpha }_j{\beta }_{1j}}\\ {\rho }_{{\beta }_{1j}{\alpha }_j}& {\sigma }_{{\beta }_{1j}}^2\end{array}\right))\text{, for sid j = 1,}\dots \text{,J}\\ {\gamma }_{1k}& \sim N({\mu }_{{\gamma }_{1k}},{\sigma }_{{\gamma }_{1k}}^2)\text{, for schid k = 1,}\dots \text{,K}\end{array}$$

lmer(math ~ year + female + mobility * size +
             (year|sid) +  (0 + mobility|schid), 
     data = d)Copy Code

Model 7

Fit the following model

$$\begin{array}{cc}{\mathrm{math}}_i& \sim N({\alpha }_{j\left[i\right]}+{\beta }_{1j\left[i\right]}(\mathrm{year}),{\sigma }^2)\\ \left(\begin{array}{c}\begin{array}{cc}& {\alpha }_j\\ & {\beta }_{1j}\end{array}\end{array}\right)& \sim N(\left(\begin{array}{c}\begin{array}{cc}& {\gamma }_0^{\alpha }+{\gamma }_1^{\alpha }\left(\mathrm{female}\right)+{\gamma }_2^{\alpha }\left(\mathrm{mobility}\right)+{\gamma }_3^{\alpha }\left(\mathrm{size}\right)\\ & {\mu }_{{\beta }_{1j}}\end{array}\end{array}\right),\left(\begin{array}{cc}{\sigma }_{{\alpha }_j}^2& {\rho }_{{\alpha }_j{\beta }_{1j}}\\ {\rho }_{{\beta }_{1j}{\alpha }_j}& {\sigma }_{{\beta }_{1j}}^2\end{array}\right))\text{, for sid j = 1,}\dots \text{,J}\end{array}$$

lmer(math ~ year + female + mobility + size +
       (year|sid), 
     data = d)Copy Code

Model 8

Last one

$$\begin{array}{cc}{\mathrm{math}}_i& \sim N({\alpha }_{j\left[i\right],k\left[i\right]}+{\beta }_{1j\left[i\right],k\left[i\right]}(\mathrm{year}),{\sigma }^2)\\ \left(\begin{array}{c}\begin{array}{cc}& {\alpha }_j\\ & {\beta }_{1j}\end{array}\end{array}\right)& \sim N(\left(\begin{array}{c}\begin{array}{cc}& {\gamma }_0^{\alpha }+{\gamma }_1^{\alpha }\left(\mathrm{female}\right)\\ & {\mu }_{{\beta }_{1j}}\end{array}\end{array}\right),\left(\begin{array}{cc}{\sigma }_{{\alpha }_j}^2& {\rho }_{{\alpha }_j{\beta }_{1j}}\\ {\rho }_{{\beta }_{1j}{\alpha }_j}& {\sigma }_{{\beta }_{1j}}^2\end{array}\right))\text{, for sid j = 1,}\dots \text{,J}\\ \left(\begin{array}{c}\begin{array}{cc}& {\alpha }_k\\ & {\beta }_{1k}\end{array}\end{array}\right)& \sim N(\left(\begin{array}{c}\begin{array}{cc}& {\gamma }_0^{\alpha }+{\gamma }_1^{\alpha }\left(\mathrm{lowinc}\right)\\ & {\gamma }_0^{{\beta }_1}+{\gamma }_1^{{\beta }_1}\left(\mathrm{lowinc}\right)\end{array}\end{array}\right),\left(\begin{array}{cc}{\sigma }_{{\alpha }_k}^2& {\rho }_{{\alpha }_k{\beta }_{1k}}\\ {\rho }_{{\beta }_{1k}{\alpha }_k}& {\sigma }_{{\beta }_{1k}}^2\end{array}\right))\text{, for schid k = 1,}\dots \text{,K}\end{array}$$

lmer(math ~ year * lowinc + female + 
       (year|sid) + (year |schid), 
     data = d)Copy Code

Updated beliefs

• Remember our posterior is defined by

$$\text{posterior}=\frac{\text{likelihood}\times \text{prior}}{\text{average likelihood}}$$ When we just had a single parameter to estimate, $\mu$, this was tractable with grid search.

With even simple linear regression, however, we have three parameters: $\alpha$, $\beta$, and $\sigma$

Our Bayesian model then becomes considerably more complicated:

$$P(\alpha ,\beta ,\sigma \mid x)=\frac{P(x\mid \alpha ,\beta ,\sigma )P(\alpha ,\beta ,\sigma )}{\iiint P(x\mid \alpha ,\beta ,\sigma )P(\alpha ,\beta ,\sigma )d\alpha d\beta d\sigma }$$

Estimation

Rather than trying to compute the integrals, we draw samples from the joint posterior distribution.

This sounds a bit like magic - how do we do this?

Multiple different algorithms, but all use some form of Markov-Chain Monte-Carlo sampling

Conceptually

• Imagine the posterior as a hill

• We start with the parameters set to random numbers

  • Estimate the posterior with this values (our spot on the hill)

• Try to "walk" around in a way that we get a complete "picture" of the hill from the samples

• Use these samples as your posterior distribution

An example

Example from Ravenzwaaij et al.

Let's go back to an example where we're trying to estimate the mean with a know standard deviation. $$N(\mu ,15)$$ Start with an initial guess, say 110

MCMC

One algorithm:

• Generate a proposal by taking a sample from your best guess according to some distribution, e.g., $\text{proposal}\sim N(110,10)$. Suppose we get a value of 108

• Observe the actual data. Let's say we have only one point: 100

• Compare the "height" of the proposal distribution to height of the current distribution, relative to observed data: $N(100|110,15)$, $N(100|108,15)$

In code

dnorm(110, 100, 15)Copy Code

## [1] 0.02129653Copy Code

dnorm(108, 100, 15)Copy Code

## [1] 0.02307026Copy Code

• If probability of target distribution is higher, accept

• If probability is lower, randomly select between the two with probability equal to the "heights" (probability of the two distributions)

• If proposal is accepted - it's the next sample in the chain

• Otherwise, the next sample is a copy of current sample

This completes one iteration

Next iteration starts by generating a new proposal distribution

Stop when you have enough samples to have an adequate "picture"

Example in code

Let's take 5000 samples

set.seed(42) # for reproducibility
samples <- c(
  110, # initial guess
  rep(NA, 4999) # space for subsequent samples to fill in
)Copy Code

Loop

for(i in 2:5000) {
  # generate proposal distribution
  proposal <- rnorm(1, mean = samples[i - 1], sd = 10)
  # calculate current/proposal distribution likelihoood
  prob_current <- dnorm(samples[i - 1], 100, 15)
  prob_proposal <- dnorm(proposal, 100, 15)
  # compute the probability ratio
  prob_ratio <- prob_proposal / prob_current
  # Determine which to select
  if(prob_ratio > runif(1)) {
    samples[i] <- proposal # accept
  } else {
    samples[i] <- samples[i - 1] # reject
  }
}Copy Code

Plot

Iteration history

tibble(iteration = 1:5000,
       value = samples) %>% 
  ggplot(aes(iteration, value)) +
  geom_line() +
  geom_hline(yintercept = 100, 
             color = "magenta",
             size = 3)Copy Code