Model Building
Data Analysis for Psychology in R 3
Josiah King
Department of Psychology
The University of Edinburgh
1 / 17

The poke in the eye

It's a trade-off...

Accurately representing the world

everything that can vary is modelled as varying

Being able to fit the model

in our sample, some things will not vary enough to fit x|g

not enough groups in g
- fit +g instead of (1|g)
predictors on different scales
- mm|g vs km|g
  - can be fixed with scaling
not enough variance in y~x between groups
model estimation hitting boundaries
- variance $\neq 0$
- correlations $\neq 1$ and $\neq -1$

2 / 17

Maximal Structures

"maximal model" = the most complex random effect structure that you can fit to the data

everything that can be modelled as a random effect, is done so

3 / 17

Maximal Structures

"maximal model" = the most complex random effect structure that you can fit to the data

everything that can be modelled as a random effect, is done so

requires sufficient variance at all levels (for both intercepts and slopes where relevant). Which is often not the case.

3 / 17

Maximal Structures

"maximal model" = the most complex random effect structure that you can fit to the data

everything that can be modelled as a random effect, is done so

requires sufficient variance at all levels (for both intercepts and slopes where relevant). Which is often not the case.

Simple(ish) design:
20 participants, each with approx 10 observations

d3 <- read_csv("https://uoepsy.github.io/data/dapr3_mindfuldecline.csv")
maxmodel <- lmer(ACE ~ visit * condition + 
                   (1 + visit | ppt), 
                 data = d3)

isSingular(maxmodel)

## [1] FALSE

Complicated design:
14 Study sites, each with c18 ppts, each with approx 10 observations

d3full <- read_csv("https://uoepsy.github.io/data/dapr3_mindfuldeclineFULL.csv")
maxmodelfull <- lmer(ACE ~ visit * condition + 
                   (1 + visit * condition | sitename) +
                   (1 + visit | sitename:ppt), 
                 data = d3full)

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

isSingular(maxmodelfull)

## [1] TRUE

3 / 17

Model Estimation

Don't report results from a model that doesn't converge. You will probably not be able to trust the estimates.

4 / 17

Model Estimation

Don't report results from a model that doesn't converge. You will probably not be able to trust the estimates.

Check model specification! Do your random effects make sense?

4 / 17

Model Estimation

Don't report results from a model that doesn't converge. You will probably not be able to trust the estimates.

Check model specification! Do your random effects make sense?

Try a different optimiser, adjust the max iterations, or the stopping tolerances

4 / 17

Model Estimation

Don't report results from a model that doesn't converge. You will probably not be able to trust the estimates.

Check model specification! Do your random effects make sense?

Try a different optimiser, adjust the max iterations, or the stopping tolerances

Try all optimisers at once! summary(allFit(model))
- look for the 'messages'
- for those optimisers that successfully converge, compare estimates.

4 / 17

Model Estimation

Don't report results from a model that doesn't converge. You will probably not be able to trust the estimates.

Check model specification! Do your random effects make sense?

Try a different optimiser, adjust the max iterations, or the stopping tolerances

Try all optimisers at once! summary(allFit(model))
- look for the 'messages'
- for those optimisers that successfully converge, compare estimates.

This probably won't help. In most cases, it is just that our model is too complex to fit to our available data.

4 / 17

Deciding on the optimal random effect structure

Keep it maximal

Start with maximal model
Remove random effects with least variance until model converges (see Barr et al., 2013)

Selection based

Use a criterion for model selection (e.g. LRT, AIC, BIC) to choose a random effect structure that is supported by the data (see Matsuchek et al., 2017)

Barr, D. J., Levy, R., Scheepers, C., & Tily, H. J. (2013). Random effects structure for confirmatory hypothesis testing: Keep it maximal. Journal of memory and language, 68(3), 255-278.
Matuschek, H., Kliegl, R., Vasishth, S., Baayen, H., & Bates, D. (2017). Balancing Type I error and power in linear mixed models. Journal of memory and language, 94, 305-315.

5 / 17

Deciding on the optimal random effect structure

Keep it maximal

Start with maximal model
Remove random effects with least variance until model converges (see Barr et al., 2013)

Selection based

Use a criterion for model selection (e.g. LRT, AIC, BIC) to choose a random effect structure that is supported by the data (see Matsuchek et al., 2017)

No right answer

Barr, D. J., Levy, R., Scheepers, C., & Tily, H. J. (2013). Random effects structure for confirmatory hypothesis testing: Keep it maximal. Journal of memory and language, 68(3), 255-278.
Matuschek, H., Kliegl, R., Vasishth, S., Baayen, H., & Bates, D. (2017). Balancing Type I error and power in linear mixed models. Journal of memory and language, 94, 305-315.

6 / 17

how to simplify

One thing at a time!

# to extract random effects
VarCorr(maxmodel)

##  Groups   Name        Std.Dev. Corr 
##  ppt      (Intercept) 0.383         
##           visit       0.277    -0.49
##  Residual             0.495

7 / 17

how to simplify

One thing at a time!

# to extract random effects
VarCorr(maxmodel)

##  Groups   Name        Std.Dev. Corr 
##  ppt      (Intercept) 0.383         
##           visit       0.277    -0.49
##  Residual             0.495

Look for:

Variances/standard deviations of 0 (or very small, e.g. 3.56e-05)

Correlations of -1 or 1

7 / 17

how to simplify (2)

One thing at a time!

# to extract random effects
VarCorr(maxmodel)

##  Groups   Name        Std.Dev. Corr 
##  ppt      (Intercept) 0.383         
##           visit       0.277    -0.49
##  Residual             0.495

Consider removing most complex random effects first.

(1 + x1 * x2 | g)
↓
(1 + x1 + x2 | g)
↓
(1 + x1 | g) OR (1 + x2 | g)

8 / 17

how to simplify (3)

One thing at a time!

# to extract random effects
VarCorr(maxmodel)

##  Groups   Name        Std.Dev. Corr 
##  ppt      (Intercept) 0.383         
##           visit       0.277    -0.49
##  Residual             0.495

Categorical predictors with >2 levels are "more complex" (because they require more parameters).

9 / 17

how to simplify (4)

One thing at a time!

# to extract random effects
VarCorr(maxmodel)

##  Groups   Name        Std.Dev. Corr 
##  ppt      (Intercept) 0.383         
##           visit       0.277    -0.49
##  Residual             0.495

If multiple levels of nesting, you'll have fewer groups as you go up the levels (fewer groups to be able to fit things to).

sometimes (but not always) it makes sense to remove higher level random effects first.

(1 + x1 + x2 | g1/g2)
is the same as
(1 + x1 + x2 | g1) + (1 + x1 + x2 | g1:g2)
↓
(1 + x1 | g1) + (1 + x1 + x2 | g1:g2)

10 / 17

how to simplify (5)

One thing at a time!

# to extract random effects
VarCorr(maxmodel)

##  Groups   Name        Std.Dev. Corr 
##  ppt      (Intercept) 0.383         
##           visit       0.277    -0.49
##  Residual             0.495

You will be faced with subjective choices

which simplification can you most easily accept?

11 / 17

random effect correlations

i've kind of been ignoring these so far..

m1 <- lmer(y ~ 1 + x1 + 
             (1 + x1 | clusterid), data = df)
VarCorr(m1)

##  Groups    Name        Std.Dev. Corr
##  clusterid (Intercept) 0.759        
##            x1          0.117    1.00
##  Residual              0.428

12 / 17

random effect correlations (2)

maxmodel <- lmer(ACE ~ visit * condition + 
                   (1 + visit | ppt), 
                 data = d3)
VarCorr(maxmodel)

##  Groups   Name        Std.Dev. Corr 
##  ppt      (Intercept) 0.383         
##           visit       0.277    -0.49
##  Residual             0.495

plot the predictions:

plot(ranef(maxmodel)$ppt)

13 / 17

random effect correlations (3)

maxmodelZC <- lmer(ACE ~ visit * condition + 
                   (1 + visit || ppt), 
                 data = d3)
VarCorr(maxmodelZC)

##  Groups   Name        Std.Dev.
##  ppt      (Intercept) 0.343   
##  ppt.1    visit       0.263   
##  Residual             0.498

removing correlations between random effects simplifies your model.

plot the predictions:

plot(ranef(maxmodelZC)$ppt)

cor(ranef(maxmodelZC)$ppt)

##             (Intercept)   visit
## (Intercept)      1.0000 -0.2384
## visit           -0.2384  1.0000

14 / 17

random effect correlations (4)

Different data! (here i have simulated the data so that people do have a very strong correlation in their intercepts and slopes)

maxmodel_COR <- lmer(ACE ~ visit * condition + 
                   (1 + visit | ppt), 
                 data = d3COR)
VarCorr(maxmodel_COR)

##  Groups   Name        Std.Dev. Corr 
##  ppt      (Intercept) 0.206         
##           visit       0.364    -1.00
##  Residual             0.105

maxmodel_CORsuppress <- lmer(ACE ~ visit * condition + 
                   (1 + visit || ppt), 
                 data = d3COR)
VarCorr(maxmodel_CORsuppress)

##  Groups   Name        Std.Dev.
##  ppt      (Intercept) 0.199   
##  ppt.1    visit       0.361   
##  Residual             0.107

ggplot(d3COR, aes(x=visit, y=ACE))+
  geom_line(aes(group=ppt, col=condition), alpha=.7)+
  theme(legend.position = "bottom")

15 / 17

random effect correlations (4)

maxmodel_COR <- lmer(ACE ~ visit * condition + 
                   (1 + visit | ppt), 
                 data = d3COR)
VarCorr(maxmodel_COR)

##  Groups   Name        Std.Dev. Corr 
##  ppt      (Intercept) 0.206         
##           visit       0.364    -1.00
##  Residual             0.105

maxmodel_CORsuppress <- lmer(ACE ~ visit * condition + 
                   (1 + visit || ppt), 
                 data = d3COR)
VarCorr(maxmodel_CORsuppress)

##  Groups   Name        Std.Dev.
##  ppt      (Intercept) 0.199   
##  ppt.1    visit       0.361   
##  Residual             0.107

plot(ranef(maxmodel_COR)$ppt)
plot(ranef(maxmodel_CORsuppress)$ppt)

cor(ranef(maxmodel_CORsuppress)$ppt)

##             (Intercept)   visit
## (Intercept)      1.0000 -0.9678
## visit           -0.9678  1.0000

15 / 17

Summary

random effect structures can get complicated quite quickly
there may multiple ways in things can be clustered.
- clusters can themselves be clustered (nesting).
- we can have different ways of clustering the same things (crossed)
the maximal random effect structure is the most complex structure we can fit to the data.
- it often leads to problems with model convergence
building lmers is a balancing act between accounting for different sources of variation and attempting to fitting a model that can be supported by our data.

16 / 17

End17 / 17

Model Building
Data Analysis for Psychology in R 3
Josiah King
Department of Psychology
The University of Edinburgh
1 / 17

The poke in the eye

It's a trade-off...

Accurately representing the world

everything that can vary is modelled as varying

Being able to fit the model

in our sample, some things will not vary enough to fit x|g

not enough groups in g
- fit +g instead of (1|g)
predictors on different scales
- mm|g vs km|g
  - can be fixed with scaling
not enough variance in y~x between groups
model estimation hitting boundaries
- variance $\neq 0$
- correlations $\neq 1$ and $\neq -1$

2 / 17

Maximal Structures

"maximal model" = the most complex random effect structure that you can fit to the data

everything that can be modelled as a random effect, is done so

3 / 17

Maximal Structures

"maximal model" = the most complex random effect structure that you can fit to the data

everything that can be modelled as a random effect, is done so

requires sufficient variance at all levels (for both intercepts and slopes where relevant). Which is often not the case.

3 / 17

Maximal Structures

"maximal model" = the most complex random effect structure that you can fit to the data

everything that can be modelled as a random effect, is done so

requires sufficient variance at all levels (for both intercepts and slopes where relevant). Which is often not the case.

Simple(ish) design:
20 participants, each with approx 10 observations

d3 <- read_csv("https://uoepsy.github.io/data/dapr3_mindfuldecline.csv")
maxmodel <- lmer(ACE ~ visit * condition + 
                   (1 + visit | ppt), 
                 data = d3)

isSingular(maxmodel)

## [1] FALSE

Complicated design:
14 Study sites, each with c18 ppts, each with approx 10 observations

d3full <- read_csv("https://uoepsy.github.io/data/dapr3_mindfuldeclineFULL.csv")
maxmodelfull <- lmer(ACE ~ visit * condition + 
                   (1 + visit * condition | sitename) +
                   (1 + visit | sitename:ppt), 
                 data = d3full)

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

isSingular(maxmodelfull)

## [1] TRUE

3 / 17

Model Estimation

Don't report results from a model that doesn't converge. You will probably not be able to trust the estimates.

4 / 17

Model Estimation

Don't report results from a model that doesn't converge. You will probably not be able to trust the estimates.

Check model specification! Do your random effects make sense?

4 / 17

Model Estimation

Don't report results from a model that doesn't converge. You will probably not be able to trust the estimates.

Check model specification! Do your random effects make sense?

Try a different optimiser, adjust the max iterations, or the stopping tolerances

4 / 17

Model Estimation

Don't report results from a model that doesn't converge. You will probably not be able to trust the estimates.

Check model specification! Do your random effects make sense?

Try a different optimiser, adjust the max iterations, or the stopping tolerances

Try all optimisers at once! summary(allFit(model))
- look for the 'messages'
- for those optimisers that successfully converge, compare estimates.

4 / 17

Model Estimation

Don't report results from a model that doesn't converge. You will probably not be able to trust the estimates.

Check model specification! Do your random effects make sense?

Try a different optimiser, adjust the max iterations, or the stopping tolerances

Try all optimisers at once! summary(allFit(model))
- look for the 'messages'
- for those optimisers that successfully converge, compare estimates.

This probably won't help. In most cases, it is just that our model is too complex to fit to our available data.

4 / 17

Deciding on the optimal random effect structure

Keep it maximal

Start with maximal model
Remove random effects with least variance until model converges (see Barr et al., 2013)

Selection based

Use a criterion for model selection (e.g. LRT, AIC, BIC) to choose a random effect structure that is supported by the data (see Matsuchek et al., 2017)

Barr, D. J., Levy, R., Scheepers, C., & Tily, H. J. (2013). Random effects structure for confirmatory hypothesis testing: Keep it maximal. Journal of memory and language, 68(3), 255-278.
Matuschek, H., Kliegl, R., Vasishth, S., Baayen, H., & Bates, D. (2017). Balancing Type I error and power in linear mixed models. Journal of memory and language, 94, 305-315.

5 / 17

Deciding on the optimal random effect structure

Keep it maximal

Start with maximal model
Remove random effects with least variance until model converges (see Barr et al., 2013)

Selection based

Use a criterion for model selection (e.g. LRT, AIC, BIC) to choose a random effect structure that is supported by the data (see Matsuchek et al., 2017)

No right answer

Barr, D. J., Levy, R., Scheepers, C., & Tily, H. J. (2013). Random effects structure for confirmatory hypothesis testing: Keep it maximal. Journal of memory and language, 68(3), 255-278.
Matuschek, H., Kliegl, R., Vasishth, S., Baayen, H., & Bates, D. (2017). Balancing Type I error and power in linear mixed models. Journal of memory and language, 94, 305-315.

6 / 17

how to simplify

One thing at a time!

# to extract random effects
VarCorr(maxmodel)

##  Groups   Name        Std.Dev. Corr 
##  ppt      (Intercept) 0.383         
##           visit       0.277    -0.49
##  Residual             0.495

7 / 17

how to simplify

One thing at a time!

# to extract random effects
VarCorr(maxmodel)

##  Groups   Name        Std.Dev. Corr 
##  ppt      (Intercept) 0.383         
##           visit       0.277    -0.49
##  Residual             0.495

Look for:

Variances/standard deviations of 0 (or very small, e.g. 3.56e-05)

Correlations of -1 or 1

7 / 17

how to simplify (2)

One thing at a time!

# to extract random effects
VarCorr(maxmodel)

##  Groups   Name        Std.Dev. Corr 
##  ppt      (Intercept) 0.383         
##           visit       0.277    -0.49
##  Residual             0.495

Consider removing most complex random effects first.

(1 + x1 * x2 | g)
↓
(1 + x1 + x2 | g)
↓
(1 + x1 | g) OR (1 + x2 | g)

8 / 17

how to simplify (3)

One thing at a time!

# to extract random effects
VarCorr(maxmodel)

##  Groups   Name        Std.Dev. Corr 
##  ppt      (Intercept) 0.383         
##           visit       0.277    -0.49
##  Residual             0.495

Categorical predictors with >2 levels are "more complex" (because they require more parameters).

9 / 17

how to simplify (4)

One thing at a time!

# to extract random effects
VarCorr(maxmodel)

##  Groups   Name        Std.Dev. Corr 
##  ppt      (Intercept) 0.383         
##           visit       0.277    -0.49
##  Residual             0.495

If multiple levels of nesting, you'll have fewer groups as you go up the levels (fewer groups to be able to fit things to).

sometimes (but not always) it makes sense to remove higher level random effects first.

(1 + x1 + x2 | g1/g2)
is the same as
(1 + x1 + x2 | g1) + (1 + x1 + x2 | g1:g2)
↓
(1 + x1 | g1) + (1 + x1 + x2 | g1:g2)

10 / 17

how to simplify (5)

One thing at a time!

# to extract random effects
VarCorr(maxmodel)

##  Groups   Name        Std.Dev. Corr 
##  ppt      (Intercept) 0.383         
##           visit       0.277    -0.49
##  Residual             0.495

You will be faced with subjective choices

which simplification can you most easily accept?

11 / 17

random effect correlations

i've kind of been ignoring these so far..

m1 <- lmer(y ~ 1 + x1 + 
             (1 + x1 | clusterid), data = df)
VarCorr(m1)

##  Groups    Name        Std.Dev. Corr
##  clusterid (Intercept) 0.759        
##            x1          0.117    1.00
##  Residual              0.428

12 / 17

random effect correlations (2)

maxmodel <- lmer(ACE ~ visit * condition + 
                   (1 + visit | ppt), 
                 data = d3)
VarCorr(maxmodel)

##  Groups   Name        Std.Dev. Corr 
##  ppt      (Intercept) 0.383         
##           visit       0.277    -0.49
##  Residual             0.495

plot the predictions:

plot(ranef(maxmodel)$ppt)

13 / 17

random effect correlations (3)

maxmodelZC <- lmer(ACE ~ visit * condition + 
                   (1 + visit || ppt), 
                 data = d3)
VarCorr(maxmodelZC)

##  Groups   Name        Std.Dev.
##  ppt      (Intercept) 0.343   
##  ppt.1    visit       0.263   
##  Residual             0.498

removing correlations between random effects simplifies your model.

plot the predictions:

plot(ranef(maxmodelZC)$ppt)

cor(ranef(maxmodelZC)$ppt)

##             (Intercept)   visit
## (Intercept)      1.0000 -0.2384
## visit           -0.2384  1.0000

14 / 17

random effect correlations (4)

Different data! (here i have simulated the data so that people do have a very strong correlation in their intercepts and slopes)

maxmodel_COR <- lmer(ACE ~ visit * condition + 
                   (1 + visit | ppt), 
                 data = d3COR)
VarCorr(maxmodel_COR)

##  Groups   Name        Std.Dev. Corr 
##  ppt      (Intercept) 0.206         
##           visit       0.364    -1.00
##  Residual             0.105

maxmodel_CORsuppress <- lmer(ACE ~ visit * condition + 
                   (1 + visit || ppt), 
                 data = d3COR)
VarCorr(maxmodel_CORsuppress)

##  Groups   Name        Std.Dev.
##  ppt      (Intercept) 0.199   
##  ppt.1    visit       0.361   
##  Residual             0.107

ggplot(d3COR, aes(x=visit, y=ACE))+
  geom_line(aes(group=ppt, col=condition), alpha=.7)+
  theme(legend.position = "bottom")

15 / 17

random effect correlations (4)

maxmodel_COR <- lmer(ACE ~ visit * condition + 
                   (1 + visit | ppt), 
                 data = d3COR)
VarCorr(maxmodel_COR)

##  Groups   Name        Std.Dev. Corr 
##  ppt      (Intercept) 0.206         
##           visit       0.364    -1.00
##  Residual             0.105

maxmodel_CORsuppress <- lmer(ACE ~ visit * condition + 
                   (1 + visit || ppt), 
                 data = d3COR)
VarCorr(maxmodel_CORsuppress)

##  Groups   Name        Std.Dev.
##  ppt      (Intercept) 0.199   
##  ppt.1    visit       0.361   
##  Residual             0.107

plot(ranef(maxmodel_COR)$ppt)
plot(ranef(maxmodel_CORsuppress)$ppt)

cor(ranef(maxmodel_CORsuppress)$ppt)

##             (Intercept)   visit
## (Intercept)      1.0000 -0.9678
## visit           -0.9678  1.0000

15 / 17

Summary

random effect structures can get complicated quite quickly
there may multiple ways in things can be clustered.
- clusters can themselves be clustered (nesting).
- we can have different ways of clustering the same things (crossed)
the maximal random effect structure is the most complex structure we can fit to the data.
- it often leads to problems with model convergence
building lmers is a balancing act between accounting for different sources of variation and attempting to fitting a model that can be supported by our data.

16 / 17

End17 / 17

↑, ←, Pg Up, k	Go to previous slide
↓, →, Pg Dn, Space, j	Go to next slide
Home	Go to first slide
End	Go to last slide
Number + Return	Go to specific slide
b / m / f	Toggle blackout / mirrored / fullscreen mode
c	Clone slideshow
p	Toggle presenter mode
t	Restart the presentation timer
?, h	Toggle this help
o	Tile View: Overview of Slides

Model Building

Data Analysis for Psychology in R 3

Josiah King

Department of PsychologyThe University of Edinburgh

The poke in the eye

Maximal Structures

Maximal Structures

Maximal Structures

Model Estimation

Model Estimation

Model Estimation

Model Estimation

Model Estimation

Deciding on the optimal random effect structure

Deciding on the optimal random effect structure

how to simplify

how to simplify

how to simplify (2)

how to simplify (3)

how to simplify (4)

how to simplify (5)

random effect correlations

random effect correlations (2)

random effect correlations (3)

random effect correlations (4)

random effect correlations (4)

Summary

End

The poke in the eye

Help

Model Building

Model Building

Data Analysis for Psychology in R 3

Josiah King

Department of PsychologyThe University of Edinburgh

The poke in the eye

Maximal Structures

Maximal Structures

Maximal Structures

Model Estimation

Model Estimation

Model Estimation

Model Estimation

Model Estimation

Deciding on the optimal random effect structure

Deciding on the optimal random effect structure

how to simplify

how to simplify

how to simplify (2)

how to simplify (3)

how to simplify (4)

how to simplify (5)

random effect correlations

random effect correlations (2)

random effect correlations (3)

random effect correlations (4)

random effect correlations (4)

Summary

End

Department of Psychology
The University of Edinburgh

Department of Psychology
The University of Edinburgh