Lecture 5: Mixed-effect regression

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

</style>

## Press record

---

## Corrections from the previous lecture

---

## R code

[LINK](https://nvaci.github.io/Lecture_5_code/Lecture_5Rcode.html)

---

## Regression
.pull-left[
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^2)$` 
b) Linearity and additivity 
c) Validity 
d) Lack of perfect multicolinearity 
]
.pull-right[
We rarely have that: [Link](http://mfviz.com/hierarchical-models/)
]
---

## 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$`

---

## Babies

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

<table>
 <thead>
 <tr>
 <th style="text-align:left;"> Country_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;"> Weight </th>
 <th style="text-align:left;"> Gender </th>
 <th style="text-align:right;"> RunningDist </th>
 </tr>
 </thead>
<tbody>
 <tr>
 <td style="text-align:left;"> 1 </td>
 <td style="text-align:right;"> 67.71321 </td>
 <td style="text-align:right;"> 4.521216 </td>
 <td style="text-align:right;"> 10 </td>
 <td style="text-align:right;"> 3861.479 </td>
 <td style="text-align:left;"> Girls </td>
 <td style="text-align:right;"> 710.8318 </td>
 </tr>
 <tr>
 <td style="text-align:left;"> 1 </td>
 <td style="text-align:right;"> 67.71321 </td>
 <td style="text-align:right;"> 4.521216 </td>
 <td style="text-align:right;"> 13 </td>
 <td style="text-align:right;"> 3679.318 </td>
 <td style="text-align:left;"> Boys </td>
 <td style="text-align:right;"> 783.3042 </td>
 </tr>
 <tr>
 <td style="text-align:left;"> 1 </td>
 <td style="text-align:right;"> 67.71321 </td>
 <td style="text-align:right;"> 4.521216 </td>
 <td style="text-align:right;"> 2 </td>
 <td style="text-align:right;"> 4428.437 </td>
 <td style="text-align:left;"> Girls </td>
 <td style="text-align:right;"> 517.1590 </td>
 </tr>
 <tr>
 <td style="text-align:left;"> 1 </td>
 <td style="text-align:right;"> 67.71321 </td>
 <td style="text-align:right;"> 4.521216 </td>
 <td style="text-align:right;"> 5 </td>
 <td style="text-align:right;"> 4156.922 </td>
 <td style="text-align:left;"> Girls </td>
 <td style="text-align:right;"> 588.9336 </td>
 </tr>
 <tr>
 <td style="text-align:left;"> 1 </td>
 <td style="text-align:right;"> 67.71321 </td>
 <td style="text-align:right;"> 4.521216 </td>
 <td style="text-align:right;"> 3 </td>
 <td style="text-align:right;"> 4611.231 </td>
 <td style="text-align:left;"> Boys </td>
 <td style="text-align:right;"> 545.8017 </td>
 </tr>
 <tr>
 <td style="text-align:left;"> 1 </td>
 <td style="text-align:right;"> 67.71321 </td>
 <td style="text-align:right;"> 4.521216 </td>
 <td style="text-align:right;"> 3 </td>
 <td style="text-align:right;"> 3436.836 </td>
 <td style="text-align:left;"> Boys </td>
 <td style="text-align:right;"> 541.9700 </td>
 </tr>
</tbody>
</table>

```r
table(Babies$Country_id)
```

```
## 
##  1  2  3  4  5  6  7  8  9 10 
## 20 20 20 20 20 20 20 20 20 20
```
---

## Analysis of the data

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

---

## Complete pooling

Linear model where we ignore information about the countries:

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

```
## 
## Call:
## lm(formula = RunningDist ~ Age, data = Babies)
## 
## Residuals:
## Min 1Q Median 3Q Max 
## -309.65 -127.68 -39.48 76.39 622.61 
## 
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) 
## (Intercept) 387.262 29.575 13.094 < 2e-16 ***
## Age 11.695 1.643 7.116 1.99e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 195.4 on 198 degrees of freedom
## Multiple R-squared: 0.2037,	Adjusted R-squared: 0.1996 
## F-statistic: 50.64 on 1 and 198 DF, p-value: 1.991e-11
```

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

## Complete pooling

Potential problems:

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

---

## No pooling

Linear model where we estimates parameters for each country separately:

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

```
## 
## Call:
## lm(formula = RunningDist ~ Age + factor(Country_id) - 1, data = Babies)
## 
## Residuals:
## Min 1Q Median 3Q Max 
## -282.176 -46.913 -2.057 41.743 259.480 
## 
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) 
## Age 12.0933 0.7851 15.404 <2e-16 ***
## factor(Country_id)1 646.5979 23.0742 28.023 <2e-16 ***
## factor(Country_id)2 335.1276 22.7026 14.762 <2e-16 ***
## factor(Country_id)3 381.3623 25.4322 14.995 <2e-16 ***
## factor(Country_id)4 243.0205 24.8691 9.772 <2e-16 ***
## factor(Country_id)5 738.8460 23.6285 31.269 <2e-16 ***
## factor(Country_id)6 220.8993 24.0903 9.170 <2e-16 ***
## factor(Country_id)7 256.4331 24.0903 10.645 <2e-16 ***
## factor(Country_id)8 469.6384 23.9834 19.582 <2e-16 ***
## factor(Country_id)9 221.2246 22.3018 9.920 <2e-16 ***
## factor(Country_id)10 296.1329 23.4261 12.641 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 90.15 on 189 degrees of freedom
## Multiple R-squared: 0.9796,	Adjusted R-squared: 0.9784 
## F-statistic: 824.4 on 11 and 189 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), ylab='Residuals',xlab='Index')
```

---

## Trade-offs between the 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

---

## Partial pooling in multilevel model

No predictors situation:

The goal is to estimate average running distance for all babies in country `$j$`, for which we have a random sample size of `$n_j$`

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}$`)

---

## Intercept

`$n_j$` - number of babies in a country 
`$\sigma^2_y$` - within country variance in running distance 
`$\sigma^2_\alpha$` - variance among the average running distance of the different countries

---

## 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^2_y)$$` 
 
Soft constraint: `$\alpha_j \sim N(\mu_\alpha, \sigma^2_\alpha)$`

---

## Intercept adjustment: R

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

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

```
## Linear mixed model fit by REML ['lmerMod']
## Formula: RunningDist ~ Age + (1 | Country_id)
##    Data: Babies
## 
## REML criterion at convergence: 2399
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.08265 -0.51390 -0.02426  0.47520  2.92714 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  Country_id (Intercept) 33092    181.91  
##  Residual                8128     90.15  
## Number of obs: 200, groups:  Country_id, 10
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 381.0106    59.2094   6.435
## Age          12.0881     0.7847  15.404
## 
## Correlation of Fixed Effects:
##     (Intr)
## Age -0.211
```

---

## Results: intercept adjustment

`$$y_i \sim N(\alpha_{j[i]}+12.08*x_i, 90.15)$$`

Soft constraint: `$\alpha_j \sim N(381.01, 181.91)$`

Variation between the countries versus within the country: 
`$\frac{\sigma^2_{\alpha}}{\sigma^2_y} = \frac{181.91}{90.15}=2$`

Our variance between the countries is 2x larger (more informative) than variance within countries - our model will always benefit from inclusion of the random adjustments

---

## Multilevel model: intercept and slope adjustments

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

---

## Intercept and slope adjustments: R

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

```
## Linear mixed model fit by REML ['lmerMod']
## Formula: RunningDist ~ Age + (1 + Age | Country_id)
##    Data: Babies
## 
## REML criterion at convergence: 1349.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4137 -0.6028 -0.0925  0.6037  3.3570 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev. Corr
##  Country_id (Intercept) 2010.65  44.840       
##             Age          111.93  10.580   0.33
##  Residual                 23.93   4.892       
## Number of obs: 200, groups:  Country_id, 10
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  399.652     14.202   28.14
## Age           11.008      3.346    3.29
## 
## Correlation of Fixed Effects:
##     (Intr)
## Age 0.330
```

---

## Fixed 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 
##   399.65175    11.00775
```

---

## Random effects

```r
ranef(mult2)
```

```
## $Country_id
##    (Intercept)        Age
## 1     69.57041  13.487071
## 2    -29.63758  -1.536129
## 3    -37.81309   2.075458
## 4     67.50182 -11.000094
## 5     38.32357  20.248860
## 6     17.97654 -10.625830
## 7    -32.89473  -5.480356
## 8    -14.06215   6.164935
## 9    -48.00096  -9.647232
## 10   -30.96383  -3.686683
## 
## with conditional variances for "Country_id"
```

---

## Why multilevel modelling

- Helps us with better inference about fixed effects
- Interest in between group effects
- Better inferences about levels of group not included in the sample

- Deals better with the non-normality of outcomes
- Deals better with outliers (shrinks them towards the mean)
- Deals better with unequal number of observations per group

---

## When multilevel models = linear regression

Very little group-level variation - complete pooling

Too much group-level variation - no pooling

We can build complete or no pooling model (next slide), but multilevel model automatically adjusts estimates for both groups, ones close to complete pooling and ones close to 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. 
---

## No pooling: 2

```r
mod1NP2<-lm(RunningDist~Age*factor(Country_id)-1, data=Babies)
summary(mod1NP2)[4]
```

```
## $coefficients
##                            Estimate Std. Error    t value      Pr(>|t|)
## Age                       24.489579  0.1273112  192.35994 2.623964e-210
## factor(Country_id)1      469.331037  2.1238473  220.98153 4.179553e-221
## factor(Country_id)2      369.926521  2.3538825  157.15590 1.337919e-194
## factor(Country_id)3      361.521752  3.5548706  101.69758 5.775693e-161
## factor(Country_id)4      467.649653  3.0302334  154.32793 3.431913e-193
## factor(Country_id)5      437.964345  2.0929519  209.25677 7.352799e-217
## factor(Country_id)6      417.770998  2.4320250  171.77907 1.648858e-201
## factor(Country_id)7      366.677101  2.5114917  145.99973 6.886336e-189
## factor(Country_id)8      385.492273  2.4490023  157.40789 1.004817e-194
## factor(Country_id)9      351.585467  2.0343006  172.82867 5.544225e-202
## factor(Country_id)10     368.627039  2.1004179  175.50176 3.556781e-203
## Age:factor(Country_id)2  -15.012744  0.2019115  -74.35308 4.665983e-137
## Age:factor(Country_id)3  -11.391698  0.2133982  -53.38236 2.727194e-112
## Age:factor(Country_id)4  -24.505671  0.1985354 -123.43222 6.831701e-176
## Age:factor(Country_id)5    6.768155  0.1706634   39.65791  6.865861e-91
## Age:factor(Country_id)6  -24.114839  0.1814550 -132.89708 1.322168e-181
## Age:factor(Country_id)7  -18.958425  0.1852506 -102.33933 1.897456e-161
## Age:factor(Country_id)8   -7.311926  0.1836773  -39.80855  3.719432e-91
## Age:factor(Country_id)9  -23.125574  0.1900993 -121.64998 9.080362e-175
## Age:factor(Country_id)10 -17.165635  0.1735666  -98.89943 8.025977e-159
```

---

## Residuals (no pooling versus multilevel model)

```r
par(mfrow=c(2,1), bty='n',mar = c(5, 4, .1, .1), cex=1.1,pch=16)
plot(resid(mult2))
plot(resid(mod1NP2))
```

---
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 reach 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

<img src="crossed.png" width="50%" style="display: block; margin: auto;" />
 
 
Random effect specification: (1|Country) + (1|Region) 
Responses have multimembership status, they are nested in a combination of countries

---

## Crossed random effects 2

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

---

## Output of the model

```r
mult2<-lmer(RunningDist~Age+(1+Age|Country_id), data=Babies)
print(summary(mult2), cor=F)
```

???

REML: [LINK](https://towardsdatascience.com/maximum-likelihood-ml-vs-reml-78cf79bef2cf) 
Intercept: Average predicted score for Running distance for Children at birth 
Age: If we compare Babies that differ in their age by 1 month, we expect to see that older babies are run by the size of `$\beta$` on average 
---

## Significance testing: fixed effects

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

```
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RunningDist ~ Age + (1 + Age | Country_id)
##    Data: Babies
## 
## REML criterion at convergence: 1349.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4137 -0.6028 -0.0925  0.6037  3.3570 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev. Corr
##  Country_id (Intercept) 2010.65  44.840       
##             Age          111.93  10.580   0.33
##  Residual                 23.93   4.892       
## Number of obs: 200, groups:  Country_id, 10
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  399.652     14.202   8.998   28.14 4.41e-10 ***
## Age           11.008      3.346   9.001    3.29  0.00938 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

???
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:

```r
anova(mult1, mult2)
```

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

```
## Data: Babies
## Models:
## mult1: RunningDist ~ Age + (1 | Country_id)
## mult2: RunningDist ~ Age + (1 + Age | Country_id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq) 
## mult1 4 2418.2 2431.4 -1205.1 2410.2 
## mult2 6 1372.6 1392.4 -680.3 1360.6 1049.6 2 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

---
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