Week 5: Mixed-effect models

class: center, inverse
background-image: url("SumerianStand2.jpg")
background-size: contain
---

</style>

## Press record

---
## Intended aims

- Motivate usage of mixed-effect models, argument differences between types of pooling and explain how multilevel models differ from linear regression models

- Build a mixed-effect model in R environment with different random adjustments, estimate and interpret the coefficients

- Derive specific level estimate and evaluate the fit of the model

---
## Regression
Gaussian distribution:
`$y_{i} \sim \mathcal{N^{iid}}(\mu_{i},\sigma)$` 
`$\mu_{i} = \alpha+\beta*x_i$` 
Model mean and standard deviation 
 
 
a) Errors: `$\mathcal{N}^{iid}(0,\sigma)$` 
b) Linearity and additivity 
c) Validity 
d) Lack of perfect multicolinearity

---

## Correlated observations

Nested structures: observations nested hierarchically under general groups (eg. answers of pupils in the UK) 
 - Pupils are nested in ther classes < schools < counties 
 
Repeated measurements: observations clustered within the participants 
- Participants answering on a number of questions in an experiment 
- NBA players contributing multiple observations of their performance throughout their career

Clustering of the observations that are not nested: 
- Questionnaire that collects profession of people and their country 
---

## Multilevel models

Generalization of regression `$(y_i=\alpha+\beta*x_i+\epsilon_i)$`, where we allow adjustments of: 
- Intercept: `$y_i = \alpha_{j[i]}+\beta*x_i+\epsilon_i$` 
- Slope: `$y_i=\alpha*\beta_{j[i]}*x_i+\epsilon_i$` 
- Both: `$y_i=\alpha_{j[i]}*\beta_{j[i]}*x_i+\epsilon_i$` 
 
Visualisation of the idea: [Link](http://mfviz.com/hierarchical-models/)
---

## How to understand multilevel models

1. __Complete pooling__ - take into account all observations, without categorical information 
2. __No pooling__ - adjust regression for each category 
3. __Partial pooling__ - use categories to adjust individual estimates

```r
knitr::kable(head(Babies), format = 'html')
```

<table>
 <thead>
 <tr>
 <th style="text-align:right;"> Babies_id </th>
 <th style="text-align:right;"> T_0s </th>
 <th style="text-align:right;"> T_1s </th>
 <th style="text-align:right;"> Age </th>
 <th style="text-align:right;"> Surounding </th>
 <th style="text-align:right;"> Weight </th>
 <th style="text-align:right;"> CalorieIntake </th>
 </tr>
 </thead>
<tbody>
 <tr>
 <td style="text-align:right;"> 1 </td>
 <td style="text-align:right;"> 27.47432 </td>
 <td style="text-align:right;"> -40.30564 </td>
 <td style="text-align:right;"> 10 </td>
 <td style="text-align:right;"> -5.375273 </td>
 <td style="text-align:right;"> 14.65666 </td>
 <td style="text-align:right;"> 665.2148 </td>
 </tr>
 <tr>
 <td style="text-align:right;"> 1 </td>
 <td style="text-align:right;"> 27.47432 </td>
 <td style="text-align:right;"> -40.30564 </td>
 <td style="text-align:right;"> 13 </td>
 <td style="text-align:right;"> 3.168467 </td>
 <td style="text-align:right;"> 12.91480 </td>
 <td style="text-align:right;"> 781.1619 </td>
 </tr>
 <tr>
 <td style="text-align:right;"> 1 </td>
 <td style="text-align:right;"> 27.47432 </td>
 <td style="text-align:right;"> -40.30564 </td>
 <td style="text-align:right;"> 2 </td>
 <td style="text-align:right;"> 43.613152 </td>
 <td style="text-align:right;"> 11.96138 </td>
 <td style="text-align:right;"> 850.2196 </td>
 </tr>
 <tr>
 <td style="text-align:right;"> 1 </td>
 <td style="text-align:right;"> 27.47432 </td>
 <td style="text-align:right;"> -40.30564 </td>
 <td style="text-align:right;"> 5 </td>
 <td style="text-align:right;"> 14.991462 </td>
 <td style="text-align:right;"> 30.50425 </td>
 <td style="text-align:right;"> 777.9477 </td>
 </tr>
 <tr>
 <td style="text-align:right;"> 1 </td>
 <td style="text-align:right;"> 27.47432 </td>
 <td style="text-align:right;"> -40.30564 </td>
 <td style="text-align:right;"> 3 </td>
 <td style="text-align:right;"> 12.254607 </td>
 <td style="text-align:right;"> 29.12190 </td>
 <td style="text-align:right;"> 703.8290 </td>
 </tr>
 <tr>
 <td style="text-align:right;"> 1 </td>
 <td style="text-align:right;"> 27.47432 </td>
 <td style="text-align:right;"> -40.30564 </td>
 <td style="text-align:right;"> 3 </td>
 <td style="text-align:right;"> 4.587097 </td>
 <td style="text-align:right;"> 22.35996 </td>
 <td style="text-align:right;"> 741.4472 </td>
 </tr>
</tbody>
</table>
---
## Complete pooling

Linear model where we ignore information about the countries:

```r
mod1CP<-lm(CalorieIntake~Age, data=Babies)
summary(mod1CP)
```

```
## 
## Call:
## lm(formula = CalorieIntake ~ Age, data = Babies)
## 
## Residuals:
## Min 1Q Median 3Q Max 
## -1544.49 -226.93 79.93 295.13 1035.00 
## 
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) 
## (Intercept) 470.846 20.356 23.13 <2e-16 ***
## Age 49.603 1.143 43.38 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 428.1 on 1998 degrees of freedom
## Multiple R-squared: 0.485,	Adjusted R-squared: 0.4848 
## F-statistic: 1882 on 1 and 1998 DF, p-value: < 2.2e-16
```

pre[class] {
 max-height: 100px;
}
</style>

## Complete pooling

Looking at residuals:

```r
par(mfrow=c(1,1), bty='n',mar = c(5, 4, .1, .1), cex=1.1, pch=16)
plot(resid(mod1CP)[1:200], ylab='Residuals', xlab='Index')
```

---

## No pooling

Linear model where we estimates parameters for each country separately:

```r
mod1NP<-lm(CalorieIntake~Age+factor(Babies_id)-1, data=Babies)
summary(mod1NP)
```

```
## 
## Call:
## lm(formula = CalorieIntake ~ Age + factor(Babies_id) - 1, data = Babies)
## 
## Residuals:
## Min 1Q Median 3Q Max 
## -809.29 -160.35 -6.41 162.10 855.07 
## 
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) 
## Age 50.2790 0.6783 74.124 < 2e-16 ***
## factor(Babies_id)1 102.5872 20.9058 4.907 9.99e-07 ***
## factor(Babies_id)2 684.2652 20.6108 33.199 < 2e-16 ***
## factor(Babies_id)3 862.3240 21.2219 40.634 < 2e-16 ***
## factor(Babies_id)4 -305.3031 20.8831 -14.620 < 2e-16 ***
## factor(Babies_id)5 487.4414 20.9251 23.295 < 2e-16 ***
## factor(Babies_id)6 132.7962 21.0883 6.297 3.72e-10 ***
## factor(Babies_id)7 659.7740 20.8901 31.583 < 2e-16 ***
## factor(Babies_id)8 656.3745 20.7328 31.659 < 2e-16 ***
## factor(Babies_id)9 679.1863 20.4171 33.266 < 2e-16 ***
## factor(Babies_id)10 642.8545 20.7242 31.019 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 253.2 on 1989 degrees of freedom
## Multiple R-squared: 0.9668,	Adjusted R-squared: 0.9666 
## F-statistic: 5259 on 11 and 1989 DF, p-value: < 2.2e-16
```

---

## No pooling

```r
par(mfrow=c(1,1), bty='n',mar = c(5, 4, .1, .1), cex=1.1, pch=16)
plot(resid(mod1NP)[1:200], ylab='Residuals',xlab='Index')
```

---

## Trade-off between the two approaches

_Complete pooling_ ignores variation between the categorical outcomes

_No pooling_ can overfit the data within each group, especially in the case when there are few observations per group or extreme cases

__Multilevel model (partial pooling)__

Multilevel estimates for a given country `$j$` can be approximated as a weighted average of the mean of observations in the country (the unpooled estimate, `$y_j$`) and the mean over all countries (completely pooled estimate `$y_{all}$`)

---

## Shrinkage estimator

Averages from countries with smaller sample size cary less information, and the weighting pulls multilevel estimates closer to the overall average. 
If the `$n_j=0$` then overal average

Averages from countries with larger sample size carry more information, and the multilevel estimates are close to the country averages. 
If `$n_j->\infty$` then country average

Every other time, multilevel estimates are between the two extremes

---

## Multilevel model: intercept

`$$y_i \sim N(\alpha_{j[i]}+\beta*x_i, \sigma_y)$$` 
Estimation: `$\alpha_j \sim N(\mu_\alpha, \sigma_\alpha)$` 
--

```r
#install.packages('lme4')
require(lme4)

mult1<-lmer(CalorieIntake~Age+(1|Babies_id), data=Babies)
summary(mult1)
```

```
## Linear mixed model fit by REML ['lmerMod']
## Formula: CalorieIntake ~ Age + (1 | Babies_id)
##    Data: Babies
## 
## REML criterion at convergence: 27858.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1994 -0.6327 -0.0245  0.6410  3.3694 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Babies_id (Intercept) 132275   363.7   
##  Residual               64118   253.2   
## Number of obs: 2000, groups:  Babies_id, 10
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 460.2558   115.6422    3.98
## Age          50.2773     0.6783   74.12
## 
## Correlation of Fixed Effects:
##     (Intr)
## Age -0.092
```

---

## Multilevel model: intercept and slope adjustments

`$$y_i \sim N(\alpha_{j[i]}+\beta_{j[i]}*x_i, \sigma^2_y)$$` 
Distributional assumptions:

<img src="insl.png" width="40%" style="display: block; margin: auto;" />
--

```r
mult2<-lmer(CalorieIntake~Age+(1+Age|Babies_id), data=Babies)
summary(mult2)
```

```
## Linear mixed model fit by REML ['lmerMod']
## Formula: CalorieIntake ~ Age + (1 + Age | Babies_id)
##    Data: Babies
## 
## REML criterion at convergence: 25474.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7037 -0.6672  0.0194  0.6576  4.0965 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev. Corr 
##  Babies_id (Intercept) 35937.8  189.57        
##            Age           701.9   26.49   -0.50
##  Residual              18938.3  137.62        
## Number of obs: 2000, groups:  Babies_id, 10
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  469.579     60.307   7.786
## Age           50.048      8.386   5.968
## 
## Correlation of Fixed Effects:
##     (Intr)
## Age -0.497
```

---
## Comparison of residuals

```r
par(mfrow=c(1,3), bty='n',mar = c(5, 4,1, .1), cex=1.1, pch=16)
plot(resid(mod1CP)[1:200], ylab='Residuals',xlab='Index', main='Complete pooling')
plot(resid(mod1NP)[1:200], ylab='Residuals',xlab='Index', main='No pooling')
plot(resid(mult2)[1:200], ylab='Residuals',xlab='Index', main='Partial pooling')
```

<img src="Week5_files/figure-html/unnamed-chunk-12-1.png" style="display: block; margin: auto;" />
---

## Fixed and random effects

Many different definitions: [LINK](https://stats.stackexchange.com/questions/4700/what-is-the-difference-between-fixed-effect-random-effect-and-mixed-effect-mode)

```r
fixef(mult2)
```

```
## (Intercept)         Age 
##   469.57946    50.04775
```

```r
ranef(mult2)
```

```
## $Babies_id
##    (Intercept)        Age
## 1    273.96266 -40.083989
## 2    -75.56727  19.518621
## 3     15.73010  22.659262
## 4   -130.81904 -40.356098
## 5    362.47248 -21.416831
## 6   -129.73763 -12.337377
## 7   -105.14316  18.855336
## 8     94.77461   6.176257
## 9   -209.91822  29.259636
## 10   -95.75452  17.725183
## 
## with conditional variances for "Babies_id"
```
---
## Intra-class correlation (ICC)

How strongly units withing the group correlate with each other

`$\frac{\sigma^2_{\alpha}}{\sigma^2_{\alpha}+\sigma^2_y} = \frac{132275}{132275+64118}=0.673522$`

Works for only intercept random effect models

```r
mod1<-lmer(CalorieIntake~Age+(1|Babies_id), data=Babies)
performance::icc(mod1)
```

```
## # Intraclass Correlation Coefficient
## 
##      Adjusted ICC: 0.674
##   Conditional ICC: 0.354
```

---

## Why multilevel modelling

1. Helps us with better inference about fixed effects 
2. Interest in between group effects 
3. Better inferences about levels of group not included in the sample 
4. Deals better with the non-normality of outcomes 
5. Deals better with outliers (shrinks them towards the mean) 
6. Deals better with unequal number of observations per group

When is multilevel model equal to linear regression:

- Very little group-level variation - complete pooling

- Too much group-level variation - no pooling

???
How many groups/observations? - understanding partial pooling as a compromise between complete and no pooling indicates that the discussion about exact number of groups and observations is misguided. With few groups, we get classical regression - complete pooling and our multilevel model should be at least as good as linear regression. Having only two observations per group or even one observation in most of the groups might result in imprecise estimate of between group variance, but it should provide partial information that allows better estimation of parameters.

---
class: inverse, middle, center
# Random effect structure 
---
## Specification of the structure

In all R packages, random effect structure can be specified using:

object <- lmer(Dependent ~ Fixed (predictor) + __(1+|...)__, data=data)

- (1|Factor) = intercept adjustments for each level of factor 
- (1+Predictor|Factor) = Intercept for each level of factor and slope adjustment for Predictor for each level of factor 
- (0 + Predictor|Factor) or (Predictor - 1|Factor) = Slope adjustment for Predictor for each level of factor

???
LMERs cheat sheet: 
- [LINK 1](https://stats.stackexchange.com/a/13173) 
- [LINK 2](https://stats.stackexchange.com/a/61466)

---

## Nested effects

<img src="nested.png" width="80%" style="display: block; margin: auto;" />
 
 
Random effect specification: (1|Country/Region) 
Responses are nested within regions, that are further nested within Countries

---

## Crossed random effects

Random effect specification: (1|Country) + (1|Region) 
Responses have multimembership status, they are nested in a combination of countries

In the experimental situations (psycholinguistics): this will be the case. Responses in experimental paradigm are nested in two random structures simultaneously: participants and stimuli

---

## How to decide what to include

There are no recipes, except that you need to think about your model and what it does

There is some advice from the experts and it is that we should keep it maximal! [LINK](https://talklab.psy.gla.ac.uk/KeepItMaximalR2.pdf) 
However, this is not always possible

???
What happens when you have singular or convergence issues: [link](https://stats.stackexchange.com/questions/110004/how-scared-should-we-be-about-convergence-warnings-in-lme4)

Additional information on random effects: [link](https://rpsychologist.com/r-guide-longitudinal-lme-lmer)

---

## Significance testing: fixed effects

```r
#install.packages('lmerTest')
require(lmerTest)
mult2<-lmer(CalorieIntake~Age+(1+Age|Babies_id), data=Babies)
print(summary(mult2), cor=F)
```

```
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CalorieIntake ~ Age + (1 + Age | Babies_id)
##    Data: Babies
## 
## REML criterion at convergence: 25474.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7037 -0.6672  0.0194  0.6576  4.0965 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev. Corr 
##  Babies_id (Intercept) 35937.8  189.57        
##            Age           701.9   26.49   -0.50
##  Residual              18938.3  137.62        
## Number of obs: 2000, groups:  Babies_id, 10
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  469.579     60.307   8.995   7.786 2.75e-05 ***
## Age           50.048      8.386   9.001   5.968 0.000211 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

???
REML: [LINK](https://towardsdatascience.com/maximum-likelihood-ml-vs-reml-78cf79bef2cf)

Satterweight approximation of degrees of freedom: [LINK](https://link.springer.com/article/10.3758/s13428-016-0809-y)
---

## Significance testing: random effects

We are not interested in understanding exactly which levels of random factors are different in comparison to the mean

We are interested in making correct estimates of our coefficients

However, we can compare whether random structure adds more information by comparing the nested models

---
class: inverse, middle, center
# Practical aspect
---

## Data

Relation between digital-screen use and the mental well-being of adolescents: [LINK](https://journals.sagepub.com/doi/10.1177/0956797616678438) 
.pull-left[
Outcome:
- Mental well-being: mwbi

Predictors:
- TV weekdays: Watchwk_adj 
- TV weekends: Watchwe_adj 
- Computer weekdays: Compwk_adj 
- Computer weekends: Compwe_adj 
- Smartphones weekday: Smartwk_adj 
- Smartphones weekends: Smartwe_adj 
- Total: sum of those
]
.pull-right[
Ethnicity: Ethnicg 
Region: Region 
Local area: LANAME

Factors:
male 
minority 
deprived 
]

---

## Reading the data in R

```r
#install.packages('foreign')
require(foreign)
mwell<-read.spss('data.sav', to.data.frame = T)
dim(mwell)
```

```
## [1] 120115     74
```

```r
mwell$Total=mwell$Watchwe_adj+mwell$Watchwk_adj+mwell$Comphwe_adj+mwell$Comphwk_adj+mwell$Smartwe_adj+mwell$Smartwk_adj
```

---

## Missing data

Outcome:

```r
table(is.na(mwell$mwb))
```

```
## 
##  FALSE   TRUE 
## 102580  17535
```

Predictor:

```r
table(is.na(mwell$Watchwk))
```

```
## 
##  FALSE   TRUE 
## 117638   2477
```

---

## Visualisations of the raw data

```r
par(mfrow=c(1,2), bty='n',mar = c(5, 4, .1, .1), cex=1.1,pch=16)
plot(density(mwell$mwb, na.rm=TRUE), main='')
plot(density(mwell$Total, na.rm = T), main='')
```

---

## Subsetting the data - excluding NAs

```r
mwell2=mwell[!is.na(mwell$mwb) & !is.na(mwell$Total),]
dim(mwell2)
```

```
## [1] 98853    75
```
---

## Scatter plots

```r
cor(mwell2$mwb, mwell2$Total)
```

```
## [1] -0.1724534
```

```r
par(mfrow=c(1,1), bty='n',mar = c(5, 4, .1, .1), cex=1.1,pch=16)
plot(mwell2$Total[1:500], mwell2$mwb[1:500])
```

---

## Observations for each category

```r
table(mwell2$Ethnicg)
```

```
## 
##            White            Mixed            Asian            Black 
##            76428             3827             9537             4413 
## Other or unknown 
##             4648
```

```r
table(mwell2$REGION)
```

```
## 
## East Midlands                                                            
##                                                                     6440 
## East of England                                                          
##                                                                     7612 
## London                                                                   
##                                                                    18524 
## North East                                                               
##                                                                     6826 
## North West                                                               
##                                                                    15276 
## South East                                                               
##                                                                    13160 
## South West                                                               
##                                                                    10105 
## West Midlands                                                            
##                                                                    10328 
## Yorkshire and the Humber                                                 
##                                                                    10582
```

---

## Observation for each category

```r
table(mwell2$LANAME)[1:10]
```

```
## 
## Barking and Dagenham                                                                                                        
##                                                                                                                         623 
## Barnet                                                                                                                      
##                                                                                                                         721 
## Barnsley                                                                                                                    
##                                                                                                                         645 
## Bath and North East Somerset                                                                                                
##                                                                                                                         549 
## Bedford                                                                                                                     
##                                                                                                                         575 
## Bexley                                                                                                                      
##                                                                                                                         678 
## Birmingham                                                                                                                  
##                                                                                                                         785 
## Blackburn with Darwen                                                                                                       
##                                                                                                                         491 
## Blackpool                                                                                                                   
##                                                                                                                         478 
## Bolton                                                                                                                      
##                                                                                                                         716
```

---

## Building the model: Random structure 1

```
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00241727 (tol = 0.002, component 1)
```

```r
anova(MWmod1, MWmod2, MWmod3)
```

```
## refitting model(s) with ML (instead of REML)
```

```
## Data: mwell2
## Models:
## MWmod1: mwb ~ (1 | LANAME)
## MWmod2: mwb ~ (1 | LANAME) + (1 | Ethnicg)
## MWmod3: mwb ~ (1 | REGION/LANAME) + (1 | Ethnicg)
##        npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)   
## MWmod1    3 726331 726359 -363162   726325                        
## MWmod2    4 726326 726364 -363159   726318 7.1022  1   0.007699 **
## MWmod3    5 726323 726371 -363157   726313 4.2095  1   0.040198 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

---

## Building the model: Fixed structure 1

```r
MWmod2a<-lmer(mwb~Total+(1|LANAME)+(1|Ethnicg), data=mwell2)
print(summary(MWmod2a), cor=F)
```

```
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mwb ~ Total + (1 | LANAME) + (1 | Ethnicg)
## Data: mwell2
## 
## REML criterion at convergence: 723288.3
## 
## Scaled residuals: 
## Min 1Q Median 3Q Max 
## -4.0004 -0.6215 0.0686 0.6822 3.0084 
## 
## Random effects:
## Groups Name Variance Std.Dev.
## LANAME (Intercept) 0.2160 0.4647 
## Ethnicg (Intercept) 0.2936 0.5418 
## Residual 87.9915 9.3804 
## Number of obs: 98853, groups: LANAME, 150; Ethnicg, 5
## 
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|) 
## (Intercept) 5.109e+01 2.594e-01 4.583e+00 196.96 3.53e-10 ***
## Total -2.054e-01 3.695e-03 9.801e+04 -55.59 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

---

## Building the model: Fixed structure 2

```r
MWmod2b<-lmer(mwb~Total+male+(1|LANAME)+(1|Ethnicg), data=mwell2)
MWmod2c<-lmer(mwb~Total*male+(1|LANAME)+(1|Ethnicg), data=mwell2)
anova(MWmod2a, MWmod2b, MWmod2c)
```

```
## refitting model(s) with ML (instead of REML)
```

```
## Data: mwell2
## Models:
## MWmod2a: mwb ~ Total + (1 | LANAME) + (1 | Ethnicg)
## MWmod2b: mwb ~ Total + male + (1 | LANAME) + (1 | Ethnicg)
## MWmod2c: mwb ~ Total * male + (1 | LANAME) + (1 | Ethnicg)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq) 
## MWmod2a 5 723288 723335 -361639 723278 
## MWmod2b 6 717197 717254 -358592 717185 6093.02 1 < 2.2e-16 ***
## MWmod2c 7 716978 717044 -358482 716964 221.25 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

---

## Results: Fixed structure 2

```r
print(summary(MWmod2c), cor=F)
```

```
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mwb ~ Total * male + (1 | LANAME) + (1 | Ethnicg)
## Data: mwell2
## 
## REML criterion at convergence: 716985.8
## 
## Scaled residuals: 
## Min 1Q Median 3Q Max 
## -4.1959 -0.6181 0.0586 0.6695 3.3167 
## 
## Random effects:
## Groups Name Variance Std.Dev.
## LANAME (Intercept) 0.1892 0.4349 
## Ethnicg (Intercept) 0.3367 0.5802 
## Residual 82.5528 9.0859 
## Number of obs: 98853, groups: LANAME, 150; Ethnicg, 5
## 
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|) 
## (Intercept) 4.884e+01 2.835e-01 5.101e+00 172.26 8.42e-11 ***
## Total -1.979e-01 4.928e-03 9.864e+04 -40.17 < 2e-16 ***
## maleYes 2.921e+00 1.328e-01 9.881e+04 21.99 < 2e-16 ***
## Total:maleYes 1.084e-01 7.285e-03 9.880e+04 14.88 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

---

## Building the model: Random structure 2

```r
MWmod3c<-lmer(mwb~Total*male+(1+Total|LANAME)+(1|Ethnicg), data=mwell2)
```

```
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient
```

```
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues
```

```
## Warning: Model failed to converge with 1 negative eigenvalue: -8.9e+00
```

```r
MWmod3c<-lmer(mwb~Total*male+(1|LANAME:male)+(1|Ethnicg), data=mwell2)
anova(MWmod2c, MWmod3c)
```

```
## refitting model(s) with ML (instead of REML)
```

```
## Data: mwell2
## Models:
## MWmod2c: mwb ~ Total * male + (1 | LANAME) + (1 | Ethnicg)
## MWmod3c: mwb ~ Total * male + (1 | LANAME:male) + (1 | Ethnicg)
##         npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## MWmod2c    7 716978 717044 -358482   716964                    
## MWmod3c    7 716995 717062 -358491   716981     0  0
```

---

## Building the model: Random structure 3

```r
MWmod3c<-lmer(mwb~Total*male+(1|LANAME)+(1|Ethnicg:male), data=mwell2)
anova(MWmod2c, MWmod3c)
```

```
## refitting model(s) with ML (instead of REML)
```

```
## Data: mwell2
## Models:
## MWmod2c: mwb ~ Total * male + (1 | LANAME) + (1 | Ethnicg)
## MWmod3c: mwb ~ Total * male + (1 | LANAME) + (1 | Ethnicg:male)
##         npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## MWmod2c    7 716978 717044 -358482   716964                    
## MWmod3c    7 716955 717022 -358471   716941 22.41  0
```

---

## Results: Full structure

```r
summary(MWmod3c)
```

```
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mwb ~ Total * male + (1 | LANAME) + (1 | Ethnicg:male)
## Data: mwell2
## 
## REML criterion at convergence: 716960.3
## 
## Scaled residuals: 
## Min 1Q Median 3Q Max 
## -4.2035 -0.6184 0.0574 0.6703 3.2764 
## 
## Random effects:
## Groups Name Variance Std.Dev.
## LANAME (Intercept) 0.1899 0.4357 
## Ethnicg:male (Intercept) 0.3687 0.6072 
## Residual 82.5235 9.0842 
## Number of obs: 98853, groups: LANAME, 150; Ethnicg:male, 10
## 
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|) 
## (Intercept) 4.897e+01 2.974e-01 9.627e+00 164.627 < 2e-16 ***
## Total -1.980e-01 4.934e-03 9.851e+04 -40.131 < 2e-16 ***
## maleYes 2.619e+00 4.167e-01 9.263e+00 6.285 0.000127 ***
## Total:maleYes 1.098e-01 7.308e-03 9.854e+04 15.023 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
## (Intr) Total maleYs
## Total -0.307 
## maleYes -0.702 0.218 
## Total:malYs 0.207 -0.673 -0.297
```

---

## Visualisation of the random structure: 1

```r
require(sjPlot)
plot_model(MWmod3c, type='re', sort.est='sort.all', grid=FALSE)[1]
```

```
## [[1]]
```

---

## Visualisation of the random structure: 2

```r
plot_model(MWmod3c, type='re', sort.est='sort.all', grid=FALSE)[2]
```

```
## [[1]]
```

---

## Visualisation of the fixed effects

```r
plot_model(MWmod3c, type='int')
```

---

## Additional information - here be dragons

Significance of the random structure

```r
ranova(MWmod3c)
```

```
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## mwb ~ Total + male + (1 | LANAME) + (1 | Ethnicg:male) + Total:male
## npar logLik AIC LRT Df Pr(>Chisq) 
## <none> 7 -358480 716974 
## (1 | LANAME) 6 -358523 717059 86.708 1 < 2.2e-16 ***
## (1 | Ethnicg:male) 6 -358525 717061 88.990 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

Explained variance - R2:

```r
#install.packages('MuMIn')
require(MuMIn)
r.squaredGLMM(MWmod3c)
```

```
##             R2m        R2c
## [1,] 0.08319687 0.08936068
```

---

## Prediction of the model

```r
mwell2$predicted=predict(MWmod3c)
par(mfrow=c(1,1), bty='n',mar = c(5, 4, .1, .1), cex=1.1,pch=16)
plot(mwell2$predicted, mwell2$mwb)
```

---
## Important aspects: theory

- Why and when would we want to use multilevel models 
- What can we include in our random structure - types of adjustments 
- Multilevel models versus complete and no pooling 
- Understanding how to specify random structure

---

## Important aspects: practice

- Building linear mixed-effect model
- Specifying random intercepts and slopes
- Comparing nested models 
- Interpreting coefficients from linear mixed-effect models

---
## Literature

Chapter 11, 12 and 13 of "Data Analysis Using Regression and Multilevel/Hierarchical Models" by Andrew Gelman and Jennifer Hill

---

# Thank you for your attention