In [1]:

```
#####
# set up the environment
#####
# load packages
library(brms)
options(mc.cores = parallel::detectCores())
# adjust the plot region for viewing on the projector in class
options(repr.plot.width=12, repr.plot.height=8)
```

Random intercept models are the simplest form of random-effects models (aka multilevel models or partial-pooling models). They allow you to say that each group has a different average outcome, but that the relationship between the outcome and the covariates is invariant across groups.

Multilevel (hierarchical) models also provide a framework for interpreting the change in an outcome variable as it relates to covariates at the level of both the individual observation and of the group. We will look at both cases here, in relation to reading test scores among grade-one students.

We will use all of the grade-one data in the Tennessee STAR data set, which has over 6000 students in more than 300 classes.

In [2]:

```
d <- read.csv("https://soci620.netlify.com/data/t_star_student_grade1.csv")
head(d)
```

student_id | teacher_id | class_type | class_size | born | female | race_ethnicity | re_white | re_black | re_asian | re_hispanic | re_native_american | re_other | reading_score | math_score | listening_score | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

<int> | <int> | <int> | <int> | <dbl> | <lgl> | <chr> | <lgl> | <lgl> | <lgl> | <lgl> | <lgl> | <lgl> | <int> | <int> | <int> | |

1 | 1 | 1 | 3 | 23 | 1979.058 | FALSE | White | TRUE | FALSE | FALSE | FALSE | FALSE | FALSE | 516 | 578 | 601 |

2 | 2 | 2 | 2 | 22 | 1980.404 | FALSE | White | TRUE | FALSE | FALSE | FALSE | FALSE | FALSE | 451 | 507 | 584 |

3 | 3 | 3 | 2 | 21 | 1979.307 | TRUE | Black | FALSE | TRUE | FALSE | FALSE | FALSE | FALSE | 483 | 526 | 529 |

4 | 4 | 4 | 3 | 21 | 1980.322 | FALSE | Black | FALSE | TRUE | FALSE | FALSE | FALSE | FALSE | 516 | 505 | 556 |

5 | 5 | 4 | 3 | 21 | 1978.416 | TRUE | Black | FALSE | TRUE | FALSE | FALSE | FALSE | FALSE | 433 | 463 | 504 |

6 | 6 | 5 | 3 | 29 | 1980.648 | TRUE | Black | FALSE | TRUE | FALSE | FALSE | FALSE | FALSE | 444 | 542 | 543 |

To begin, we will use the same random intercept model (with *no* covariates) as we used in the previous lab:

In [3]:

```
# model
model_interceptonly <- bf(
reading_score ~ 1 + (1 | teacher_id)
)
# prior
prior_interceptonly <- c(
prior(cauchy(0,50),class='sigma'),
prior(normal(500,100),class="Intercept"),
prior(cauchy(0,50),class='sd')
)
```

In [4]:

```
fit_interceptonly <- brm(model_interceptonly, data=d,prior=prior_interceptonly, iter = 6000)
summary(fit_interceptonly)
```

Warning message: “Rows containing NAs were excluded from the model.” Compiling Stan program... Start sampling

Family: gaussian Links: mu = identity; sigma = identity Formula: reading_score ~ 1 + (1 | teacher_id) Data: d (Number of observations: 6395) Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1; total post-warmup samples = 12000 Group-Level Effects: ~teacher_id (Number of levels: 334) Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS sd(Intercept) 28.78 1.27 26.34 31.33 1.00 3259 5526 Population-Level Effects: Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS Intercept 521.82 1.69 518.54 525.08 1.00 2255 4041 Family Specific Parameters: Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS sigma 47.28 0.43 46.44 48.14 1.00 20480 8176 Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS and Tail_ESS are effective sample size measures, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat = 1).

Even without predictor variables, this model shows us that there is a significant amount of between-class variation (as estimated in `phi`

).

We will start by seeing how students' age affects their test scores, taking into account inter-class variation.

In [5]:

```
# construct a standardized age variable from the `born` variable
d$age_s <- scale(-d$born)
plot(density(d$age_s,na.rm=TRUE))
```

In [6]:

```
# model
model_age <- bf(
reading_score ~ 1 + age_s + (1 | teacher_id)
)
get_prior(model_age,data=d)
```

Warning message: “Rows containing NAs were excluded from the model.”

prior | class | coef | group | resp | dpar | nlpar | bound | source |
---|---|---|---|---|---|---|---|---|

<chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> |

b | default | |||||||

b | age_s | default | ||||||

student_t(3, 514, 57.8) | Intercept | default | ||||||

student_t(3, 0, 57.8) | sd | default | ||||||

sd | teacher_id | default | ||||||

sd | Intercept | teacher_id | default | |||||

student_t(3, 0, 57.8) | sigma | default |

In [7]:

```
prior_age <- c(
prior(cauchy(0,50),class='sigma'),
prior(normal(500,100),class="Intercept"),
prior(cauchy(0,50),class='sd'),
prior(cauchy(0,30),class='b') # for age
)
```

In [8]:

```
fit_age <- brm(model_age, data=d,prior=prior_age, iter = 6000)
summary(fit_age)
```

Family: gaussian Links: mu = identity; sigma = identity Formula: reading_score ~ 1 + age_s + (1 | teacher_id) Data: d (Number of observations: 6392) Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1; total post-warmup samples = 12000 Group-Level Effects: ~teacher_id (Number of levels: 334) Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS sd(Intercept) 28.69 1.29 26.27 31.32 1.00 2954 5454 Population-Level Effects: Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS Intercept 521.63 1.70 518.28 524.96 1.00 2063 4168 age_s -5.36 0.63 -6.59 -4.12 1.00 20700 9131 Family Specific Parameters: Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS sigma 47.02 0.42 46.20 47.86 1.00 20188 9087 Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS and Tail_ESS are effective sample size measures, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat = 1).

Age has a negative relationship with test scores: or every increase of one standard deviation increase in age, the model predicts that student will score about five points lower.

Very often, it is important to be able to look at the relationship between an outcome variable and a covariate that only varies at the group level. This is exactly the case in the STAR data if we are interested in the effect of the experimental treatment. Do smaller classes lead to better test scores? A model that estimates this effect, allowing between-class variation inn average scores, be written mathematically (omitting priors for now) as

$$ \begin{aligned} S_{ik} &\sim \mathrm{Norm}(\mu_{ik},\sigma)\\ \mu_{ik} &= \beta_{0k} + \beta_{1}Age_i\\ \\ \beta_{0k} &= \gamma_{0} + \gamma_1Small_k+ \eta_{0k}\\ \eta_{0k} &\sim \mathrm{Norm}(0, \phi_0) \end{aligned}$$To determine how to specify this model in `brms`

, note that we can substitute $\beta_{0k}$ into the individual-level linear model and rearrange, giving us

The classroom-level covariate can be specified as a student-level covariate in this case, so the changes to the formula are quite minimal:

In [9]:

```
# create an indicator variable for 'small' classes
d$small <- d$class_type == 1
# specify the model
model_age_size <- bf(
reading_score ~ 1 + age_s + small + (1 | teacher_id)
)
# look at what kind of prior we need
get_prior(model_age_size,data=d)
```

Warning message: “Rows containing NAs were excluded from the model.”

prior | class | coef | group | resp | dpar | nlpar | bound | source |
---|---|---|---|---|---|---|---|---|

<chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> |

b | default | |||||||

b | age_s | default | ||||||

b | smallTRUE | default | ||||||

student_t(3, 514, 57.8) | Intercept | default | ||||||

student_t(3, 0, 57.8) | sd | default | ||||||

sd | teacher_id | default | ||||||

sd | Intercept | teacher_id | default | |||||

student_t(3, 0, 57.8) | sigma | default |

In [10]:

```
prior_age_size <- c(
prior(cauchy(0,50),class='sigma'),
prior(normal(500,100),class="Intercept"),
prior(cauchy(0,50),class='sd'),
prior(cauchy(0,30),coef='age_s'), # for age
prior(cauchy(0,50),coef='smallTRUE') # for small classes
)
```

In [11]:

```
fit_age_size <- brm(model_age_size, data=d,prior=prior_age_size, iter = 6000)
```

Error in summary(fit_age_smize): object 'fit_age_smize' not found Traceback: 1. summary(fit_age_smize)

In [12]:

```
summary(fit_age_size)
```