Agenda

• Review some notation

• Unstructured VCV Matrices and alternatives

Learning Objectives

• Gain a deeper understanding of how the residual structure is different in multilevel models

• Understand that there are methods for changing the residual structure, and understand when and why this might be preferable

• Be able to implement alternative methods using {nlme}

Data

Let's start by using the sleepstudy data from {lme4}.

library(lme4)
head(sleepstudy)Copy Code

##   Reaction Days Subject
## 1 249.5600    0     308
## 2 258.7047    1     308
## 3 250.8006    2     308
## 4 321.4398    3     308
## 5 356.8519    4     308
## 6 414.6901    5     308Copy Code

Translate

Translate the following model into lme4::lmer() syntax

$$\begin{array}{cc}{\mathrm{Reaction}}_i& \sim N({\alpha }_{j\left[i\right]},{\sigma }^2)\\ {\alpha }_j& \sim N({\mu }_{{\alpha }_j},{\sigma }_{{\alpha }_j}^2)\text{, for Subject j = 1,}\dots \text{,J}\end{array}$$

lmer(Reaction ~ 1 + (1|Subject), data = sleepstudy)Copy Code

More complicated

Translate this equation into lme4::lmer() syntax

$$\begin{array}{cc}{\mathrm{Reaction}}_i& \sim N({\alpha }_{j\left[i\right]}+{\beta }_1(\mathrm{Days}),{\sigma }^2)\\ {\alpha }_j& \sim N({\mu }_{{\alpha }_j},{\sigma }_{{\alpha }_j}^2)\text{, for Subject j = 1,}\dots \text{,J}\end{array}$$

lmer(Reaction ~ Days + (1|Subject), data = sleepstudy)Copy Code

Even more complicated

Translate this equation into lme4::lmer() syntax

$$\begin{array}{cc}{\mathrm{Reaction}}_i& \sim N({\alpha }_{j\left[i\right]}+{\beta }_{1j\left[i\right]}(\mathrm{Days}),{\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}& {\mu }_{{\alpha }_j}\\ & {\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 Subject j = 1,}\dots \text{,J}\end{array}$$

lmer(Reaction ~ Days + (Days|Subject), data = sleepstudy)Copy Code

New data

sim3 <- read_csv(here::here("data", "sim3level.csv"))
sim3Copy Code

## # A tibble: 570 x 6
##        Math ActiveTime ClassSize Classroom School StudentID
##       <dbl>      <dbl>     <dbl>     <dbl> <chr>      <dbl>
##  1 55.42604 0.06913359        18         1 Sch1           1
##  2 54.34306 0.08462063        18         1 Sch1           2
##  3 61.42570 0.1299456         18         1 Sch1           3
##  4 56.12271 0.7461320         18         1 Sch1           4
##  5 53.34900 0.03887918        18         1 Sch1           5
##  6 57.99773 0.6856354         18         1 Sch1           6
##  7 54.46326 0.1439774         18         1 Sch1           7
##  8 64.20517 0.8910800         18         1 Sch1           8
##  9 54.49117 0.08963612        18         1 Sch1           9
## 10 57.87599 0.03773272        18         1 Sch1          10
## # … with 560 more rowsCopy Code

Data obtained here

A bit of a problem

sim3 %>% 
  count(Classroom, School)Copy Code

## # A tibble: 30 x 3
##    Classroom School     n
##        <dbl> <chr>  <int>
##  1         1 Sch1      18
##  2         1 Sch2      18
##  3         1 Sch3      19
##  4         2 Sch1      14
##  5         2 Sch2      18
##  6         2 Sch3      24
##  7         3 Sch1      20
##  8         3 Sch2      15
##  9         3 Sch3      20
## 10         4 Sch1      14
## # … with 20 more rowsCopy Code

Make classroom unique

We could handle this with our model syntax, but nested IDs like this always makes me nervous. Let's make them unique.

sim3 <- sim3 %>% 
  mutate(class_id = paste0("class", Classroom, ":", School))
sim3Copy Code

## # A tibble: 570 x 7
##        Math ActiveTime ClassSize Classroom School StudentID class_id   
##       <dbl>      <dbl>     <dbl>     <dbl> <chr>      <dbl> <chr>      
##  1 55.42604 0.06913359        18         1 Sch1           1 class1:Sch1
##  2 54.34306 0.08462063        18         1 Sch1           2 class1:Sch1
##  3 61.42570 0.1299456         18         1 Sch1           3 class1:Sch1
##  4 56.12271 0.7461320         18         1 Sch1           4 class1:Sch1
##  5 53.34900 0.03887918        18         1 Sch1           5 class1:Sch1
##  6 57.99773 0.6856354         18         1 Sch1           6 class1:Sch1
##  7 54.46326 0.1439774         18         1 Sch1           7 class1:Sch1
##  8 64.20517 0.8910800         18         1 Sch1           8 class1:Sch1
##  9 54.49117 0.08963612        18         1 Sch1           9 class1:Sch1
## 10 57.87599 0.03773272        18         1 Sch1          10 class1:Sch1
## # … with 560 more rowsCopy Code

New model

Translate this equation into code

$$\begin{array}{cc}{\mathrm{Math}}_i& \sim N({\alpha }_{j\left[i\right],k\left[i\right]}+{\beta }_{1j\left[i\right]}(\mathrm{ActiveTime}),{\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}& {\mu }_{{\alpha }_j}\\ & {\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 class\_id j = 1,}\dots \text{,J}\\ {\alpha }_k& \sim N({\mu }_{{\alpha }_k},{\sigma }_{{\alpha }_k}^2)\text{, for School k = 1,}\dots \text{,K}\end{array}$$

lmer(Math ~ ActiveTime + (ActiveTime|class_id) + (1|School),
     data = sim3)Copy Code

This one

What about this 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{ActiveTime}),{\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{ClassSize}\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 class\_id 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}& {\mu }_{{\alpha }_k}\\ & {\mu }_{{\beta }_{1k}}\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 School k = 1,}\dots \text{,K}\end{array}$$

lmer(Math ~ ActiveTime + ClassSize + 
       (ActiveTime|class_id) + (ActiveTime|School),
     data = sim3)Copy Code

Last one

What about this 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{ActiveTime}),{\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{ClassSize}\right)\\ & {\gamma }_0^{{\beta }_1}+{\gamma }_1^{{\beta }_1}\left(\mathrm{ClassSize}\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 class\_id 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}& {\mu }_{{\alpha }_k}\\ & {\mu }_{{\beta }_{1k}}\\ & {\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 School k = 1,}\dots \text{,K}\end{array}$$

lmer(Math ~ ActiveTime * ClassSize + 
       (ActiveTime|class_id) + (ActiveTime + ClassSize|School),
     data = sim3)Copy Code

Data

Willett, 1988

• $n$ = 35 people
• Each completed a cognitive inventory on "opposites naming"
• At first time point, participants also completed a general cognitive measure

Read in data

willett <- read_csv(here::here("data", "willett-1988.csv"))
willettCopy Code

## # A tibble: 140 x 4
##       id  time   opp   cog
##    <dbl> <dbl> <dbl> <dbl>
##  1     1     0   205   137
##  2     1     1   217   137
##  3     1     2   268   137
##  4     1     3   302   137
##  5     2     0   219   123
##  6     2     1   243   123
##  7     2     2   279   123
##  8     2     3   302   123
##  9     3     0   142   129
## 10     3     1   212   129
## # … with 130 more rowsCopy Code

Standard OLS

• We have four observations per participant.

• If we fit a standard OLS model, it would look like this

bad <- lm(opp ~ time, data = willett)
summary(bad)Copy Code

## 
## Call:
## lm(formula = opp ~ time, data = willett)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -88.374 -25.584   1.186  28.926  64.746 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  164.374      5.035   32.65   <2e-16 ***
## time          26.960      2.691   10.02   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 35.6 on 138 degrees of freedom
## Multiple R-squared:  0.421,    Adjusted R-squared:  0.4168 
## F-statistic: 100.3 on 1 and 138 DF,  p-value: < 2.2e-16Copy Code

Assumptions

As we discussed previously, this model looks like this

$$\mathrm{opp}=\alpha +{\beta }_1\left(\mathrm{time}\right)+ϵ$$

where

$$ϵ\sim (0,\sigma )$$

Individual level residuals

We can expand our notation, so it looks like a multivariate normal distribution

$$\left(\begin{array}{c}{ϵ}_1\\ {ϵ}_2\\ {ϵ}_3\\ ⋮\\ {ϵ}_n\end{array}\right)\sim MVN(\left[\begin{array}{c}0\\ 0\\ 0\\ ⋮\\ 0\end{array}\right],\left[\begin{array}{ccccc}{\sigma }_{ϵ}& 0& 0& \dots & 0\\ 0& {\sigma }_{ϵ}& 0& 0& 0\\ 0& 0& {\sigma }_{ϵ}& 0& 0\\ ⋮& 0& 0& \ddots & ⋮\\ 0& 0& 0& \dots & {\sigma }_{ϵ}\end{array}\right])$$

This is where the $i.i.d.$ part comes in. The residuals are assumed $i$ndependent and $i$dentically $d$istributed.

Multilevel model

Very regularly, there are reasons to believe the $i.i.d.$ assumption is violated. Consider our current case, with 4 time points for each individual.

• Is an observation for one time point for one individual independent from the other observations for that individual?

• Rather than estimating a single residual variance, we estimate an additional components associated with individuals, leading to a block diagonal structure