---
title: "Intro to Multilevel Model: Lecture 9 Example Code"
author: "Sunthud Pornprasertmanit"
output: html_document
date: "`r format(Sys.time(), '%d/%m/%Y')`"
---

```{css style settings, echo = FALSE}
blockquote {
    padding: 10px 20px;
    margin: 0 0 20px;
    font-size: 14px;
    border-left: 5px solid #eee;
    font-family: "Courier New";
}
```

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r xx, include=FALSE}
library(lme4)
library(effects)
library(interactions)
library(ggplot2)
library(psych)
library(r2mlm)
# out <- readr::read_rds("lecture7saveresults.rds")
# ss41 <- out[[1]]
# ss42 <- out[[2]]
# ss31 <- out[[3]]
```

## Build-up Strategy

##### Example 1: Goldfish Food Consumption within Fish Tanks

Read the data set from Lecture 4.

```{r readdata1}
dat1 <- read.table("lecture4ex1.csv", sep=",", header=TRUE)
```

Transform the fish length variable by 10 to make the multilevel analysis easier to converge. Next, the fish length and fish color are group mean centered.

```{r calcgroupmean1}
dat1$lengthx <- dat1$length/10
dat1$avelengthx <- ave(dat1$lengthx, dat1$groupid)
dat1$difflengthx <- dat1$lengthx - dat1$avelengthx
dat1$avegold <- ave(dat1$goldcolor, dat1$groupid)
dat1$diffgold <- dat1$goldcolor - dat1$avegold
```

Let's start using build-up strategy to examine the relationship between predictors and criterion. Start with the null model (Model 0).

```{r null1}
library(lme4)
library(r2mlm)
out1b0 <- lmer(consume ~ 1 + (1|groupid), data=dat1, REML=FALSE)
r2mlm(out1b0, bargraph=FALSE)
```

Next, add level-1 predictors: `difflengthx` (Model 1) and `diffgold` (Model 2).

```{r addl1predictor1}
out1b1 <- lmer(consume ~ 1 + difflengthx + (1|groupid), data=dat1, REML=FALSE)
anova(out1b0, out1b1)
r2mlm(out1b1, bargraph=FALSE)
out1b2 <- lmer(consume ~ 1 + diffgold + (1|groupid), data=dat1, REML=FALSE)
anova(out1b0, out1b2)
r2mlm(out1b2, bargraph=FALSE)
```

Only fish length was significant. Thus Model 1 is retained. The interaction effect between level-1 predictors is investigated. The full model is the model with interaction (Model 4) and the restricted model is the model with both level-1 predictors without interaction.

```{r l1interaction1}
out1b3 <- lmer(consume ~ 1 + difflengthx + diffgold + (1|groupid), data=dat1, REML=FALSE)
out1b4 <- lmer(consume ~ 1 + difflengthx + diffgold + difflengthx:diffgold + (1|groupid), data=dat1, REML=FALSE)
anova(out1b3, out1b4)
```

The interaction was not significant. Thus, Model 1 is still retained. Next, check the effect of level-2 predictors: `volume` (Model 5), `nfish` (Model 6), `avelengthx` (Model 7), `avegold` (Model 8). 

```{r addl2predictor12}
out1b5 <- lmer(consume ~ 1 + difflengthx + volume + (1|groupid), data=dat1, REML=FALSE)
anova(out1b1, out1b5)
out1b6 <- lmer(consume ~ 1 + difflengthx + nfish + (1|groupid), data=dat1, REML=FALSE)
anova(out1b1, out1b6)
out1b7 <- lmer(consume ~ 1 + difflengthx + avelengthx + (1|groupid), data=dat1, REML=FALSE)
anova(out1b1, out1b7)
out1b8 <- lmer(consume ~ 1 + difflengthx + avegold + (1|groupid), data=dat1, REML=FALSE)
anova(out1b1, out1b8)
```

The effects of `avelengthx` and `volume` were significant. Model 9 is then created to include only both significant effects. 

```{r completel2predictor1}
out1b9 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + (1|groupid), data=dat1, REML=FALSE)
```

Next, the interactions between level 2 predictors are investigated. Four level-2 predictors can make six two-way interactions. 

1. `avelengthx*volume` (Model 10) was not significant. The restricted model is Model 9.

```{r l2interaction11}
out1b10 <- lmer(consume ~ 1 + difflengthx + avelengthx*volume + (1|groupid), data=dat1, REML=FALSE)
anova(out1b9, out1b10)
```

2. `avelengthx*nfish` (Model 12) was not significant. The restricted model is Model 9 with `nfish` (Model 11).
3. `volume*nfish` (Model 13) was not significant. The restricted model is Model 9 with `nfish` (Model 11).

```{r l2interaction12}
out1b11 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + nfish + (1|groupid), data=dat1, REML=FALSE)
out1b12 <- lmer(consume ~ 1 + difflengthx + volume + avelengthx*nfish + (1|groupid), data=dat1, REML=FALSE)
out1b13 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume*nfish + (1|groupid), data=dat1, REML=FALSE)
anova(out1b11, out1b12)
anova(out1b11, out1b13)
```

4. `avelengthx*avegold` (Model 15) was not significant. The restricted model is Model 9 with `avegold` (Model 14).
5. `volume*avegold` (Model 16) was not significant. The restricted model is Model 9 with `avegold` (Model 14).

```{r l2interaction13}
out1b14 <- lmer(consume ~ 1 + difflengthx + volume + avelengthx + avegold + (1|groupid), data=dat1, REML=FALSE)
out1b15 <- lmer(consume ~ 1 + difflengthx + volume + avelengthx*avegold + (1|groupid), data=dat1, REML=FALSE)
out1b16 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume*avegold + (1|groupid), data=dat1, REML=FALSE)
anova(out1b14, out1b15)
anova(out1b14, out1b16)
```

6. `avegold*nfish` (Model 18) was not significant. The restricted model is Model 9 with `avegold` and `nfish` (Model 17).

```{r l2interaction14}
out1b17 <- lmer(consume ~ 1 + difflengthx + diffgold + volume + avelengthx + avegold + nfish + (1|groupid), data=dat1, REML=FALSE)
out1b18 <- lmer(consume ~ 1 + difflengthx + diffgold + volume + avelengthx + avegold*nfish + (1|groupid), data=dat1, REML=FALSE)
anova(out1b17, out1b18)
```

All two-way interactions between level-2 predictors are not significant. Thus, Model 9 is still retained. Next, the random effects of two level-1 predictors are investigated. First, the difference slopes of `difflengthx` (Model 19) is compared with Model 9.

```{r randomslope11}
out1b19 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
anova(out1b9, out1b19)
```

Second, the difference slopes of `diffgold` (Model 21) is compared with Model 9 with the fixed effect of `diffgold` (Model 20). 

```{r randomslope12}
out1b20 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + volume + (1|groupid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
out1b21 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + volume + (1 + goldcolor|groupid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1b20, out1b21)
```

The random slope of `difflengthx` was the only significant effect. Thus, Model 19 is retained.

Next, the cross-level interactions were investigated. There are four level-2 predictors that can predict the random slopes of `difflengthx`.

1. `difflengthx:avelengthx` (Model 22) was not significant. The restricted model is Model 19.
2. `difflengthx:volume` (Model 23) was significant. The restricted model is Model 19.

```{r crosslevel11}
out1b22 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + difflengthx:avelengthx + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
anova(out1b19, out1b22)
out1b23 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + difflengthx:volume + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
anova(out1b19, out1b23)
```

3. `difflengthx:nfish` (Model 25) was not significant. The restricted model is Model 19 with `nfish` (Model 24).

```{r crosslevel12}
out1b24 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + nfish + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
out1b25 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + nfish + difflengthx:nfish + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
anova(out1b24, out1b25)
```

4. `difflengthx:avegold` (Model 27) was not significant. The restricted model is Model 19 with `avegold` (Model 26).

```{r crosslevel13}
out1b26 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + avegold + (1 + difflengthx|groupid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
out1b27 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + avegold + difflengthx:avegold + (1 + difflengthx|groupid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1b26, out1b27)
```

`difflengthx:volume` was the only significant cross-level interaction. Therefore, Model 23 is the final model from the build-up strategy.

## Tear-down strategy.

##### Example 1: Goldfish Food Consumption within Fish Tanks

The starting model (Model 0) for the tear-down strategy has all level-1 (fixed and random effects) and level-2 (fixed effects) predictors as well as all cross-level interactions (fixed effects). Note that the interactions among level-1 predictors and interactions among level-2 predictors are not included to simplify the example.

```{r startingmodel2}
out1t0 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold|groupid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
```

First, the cross-level interactions will be investigated. The significant interactions will be retained and the nonsignificant interactions will be dropped. Let's check each interaction as follows:

1. `difflengthx:avelengthx` was not significant. However, the model that drop this interaction was not convergent. The problem is fixed by multiplying `difflengthx` by 10. Then, the model with the interaction (Model 0a) and without the interaction (Model 1a) was compared.

```{r crosslevel21}
dat1a <- dat1
dat1a$difflengthx <- dat1a$difflengthx * 10
out1t1a <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avegold + difflengthx:volume + difflengthx:nfish + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold|groupid), data=dat1a, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
out1t0a <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold|groupid), data=dat1a, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1t1a, out1t0a)
```

2. `difflengthx:avegold` (Model 2) was not significant when compared with Model 0.

```{r crosslevel22}
out1t2 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:volume + difflengthx:nfish + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold|groupid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1t2, out1t0)
```

3. `difflengthx:volume` (Model 3) was not significant when compared with Model 0.

```{r crosslevel23}
out1t3 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:nfish + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold|groupid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1t3, out1t0)
```

4. `difflengthx:nfish` (Model 4) was not significant when compared with Model 0.

```{r crosslevel24}
out1t4 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold|groupid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1t4, out1t0)
```

5. `diffgold:avelengthx` (Model 5) was not significant when compared with Model 0.

```{r crosslevel25}
out1t5 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish + diffgold:avegold + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold|groupid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1t5, out1t0)
```

6. `diffgold:avegold` (Model 6) was not significant when compared with Model 0.

```{r crosslevel26}
out1t6 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish + diffgold:avelengthx + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold|groupid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1t6, out1t0)
```

7. `diffgold:volume` (Model 7) was not significant when compared with Model 0.

```{r crosslevel27}
out1t7 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish + diffgold:avelengthx + diffgold:avegold + diffgold:nfish + (1 + difflengthx + diffgold|groupid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1t7, out1t0)
```

8. `diffgold:nfish` (Model 8) was not significant when compared with Model 0.

```{r crosslevel28}
out1t8 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish + diffgold:avelengthx + diffgold:avegold + diffgold:volume + (1 + difflengthx + diffgold|groupid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1t8, out1t0)
```

All cross-level interactions were not significant so they are dropped from the model as Model 9.

```{r crosslevel29}
out1t9 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + (1 + difflengthx + diffgold|groupid), data=dat1, REML=FALSE)
```

Next, the random effects of both level-1 predictors were investigated. 

```{r randomeff2}
out1t10 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + (1 + diffgold|groupid), data=dat1, REML=FALSE)
anova(out1t10, out1t9)
out1t11 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
anova(out1t11, out1t9)
```

The random effect of `difflengthx` was significant but the random effect of `diffgold` was not significant. Therefore, Model 11 was retained.

Next, all remaining fixed effects were tested except `difflengthx` (because the random effect was significant). 

1. Dropping `diffgold` (Model 12) was significant so the effect will be retained.

```{r fixedeff21}
out1t12 <- lmer(consume ~ 1 + difflengthx + avelengthx + avegold + volume + nfish + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
anova(out1t12, out1t11)
```

2. Dropping `avelengthx` (Model 13) was significant so the effect will be retained.

```{r fixedeff22}
out1t13 <- lmer(consume ~ 1 + difflengthx + diffgold + avegold + volume + nfish + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
anova(out1t13, out1t11)
```

3. Dropping `avegold` (Model 14) was significant so the effect will be retained.

```{r fixedeff23}
out1t14 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + volume + nfish + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
anova(out1t14, out1t11)
```

4. Dropping `volume` (Model 15) was significant so the effect will be retained.

```{r fixedeff24}
out1t15 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + nfish + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
anova(out1t15, out1t11)
```

5. Dropping `nfish` (Model 16) was significant so the effect will be retained.

```{r fixedeff25}
out1t16 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
anova(out1t16, out1t11)
```

All fixed effects were significant. Thus, the final model is Model 11.

## Mixed Strategy with a priori Hypotheses

##### Example 2: The Interaction Effect between Interviewers' and Interviewees' Sex on Interview Ratings

Read the new data set which is the original interview dataset in Lecture 4 with additional covariates.

```{r readdata3}
dat2 <- read.table("lecture9ex1.csv", sep=",", header=TRUE)
```

Based on the research hypothesis, the interaction of interviewers' and interviewees' sex, the random effect of the interviewee's sex, and all fixed effects of all covariates are included in the model (Model 0).

```{r hypothesis3}
library(lme4)
out2m0 <- lmer(score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + erexp + eesex:ersex + (1 + eesex|erid), data=dat2, REML=FALSE)
summary(out2m0)
```

First, all covariates are checked whether they are needed in the model. There are four fixed effects that can be dropped (excluding `eesex` and `ersex`):

1. To check whether `eesalary` can be dropped. First, check whether `eesalary` has a random effect (Model 1). The random effect was significant so the random effect is retained. 

```{r dropfixed31}
out2m1 <- lmer(score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + erexp + eesex:ersex + (1 + eesex + eesalary|erid), data=dat2, REML=FALSE)
anova(out2m1, out2m0)
```

2. To check whether `eeworkexp` can be dropped. First, check whether `eeworkexp` has a random effect (Model 2). The random effect was not significant. Next, dropping `eeworkexp` (Model 3) from the model was significant so this fixed effect will be retained.

```{r dropfixed32}
out2m2 <- lmer(score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + erexp + eesex:ersex + (1 + eesex + eeworkexp|erid), data=dat2, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out2m2, out2m0)
out2m3 <- lmer(score ~ 1 + eesex + eesalary + I(eeiq - 100) + ersex + erexp + eesex:ersex + (1 + eesex|erid), data=dat2, REML=FALSE)
anova(out2m3, out2m0)
```

3. To check whether `eeiq` can be dropped. First, check whether `eeworkexp` has a random effect (Model 4). After some transformations to make the model convergent, the random effect was not significant. Next, dropping `eeiq` (Model 5) from the model was significant so this fixed effect will be retained.
  

```{r dropfixed33}
out2m0a <- lmer(score ~ 1 + eesex + eesalary + I(eeworkexp/10) + I(eeiq/10 - 10) + ersex + erexp + eesex:ersex + (1 + eesex|erid), data=dat2, REML=FALSE)
out2m4a <- lmer(score ~ 1 + eesex + eesalary + I(eeworkexp/10) + I(eeiq/10 - 10) + ersex + erexp + eesex:ersex + (1 + eesex + I(eeiq/10 - 10)|erid), data=dat2, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out2m4a, out2m0a)
out2m5 <- lmer(score ~ 1 + eesex + eesalary + eeworkexp + ersex + erexp + eesex:ersex + (1 + eesex|erid), data=dat2, REML=FALSE)
anova(out2m5, out2m0)
```

4. `erexp` is the level-2 predictor so the random effect is not tested. Dropping `erexp` (Model 6) was not significant so this fixed effect will be dropped. 

```{r dropfixed34}
out2m6 <- lmer(score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + eesex:ersex + (1 + eesex|erid), data=dat2, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out2m6, out2m0)
```

From the results, as shown in Model 7, the random effect of `eesalary` will be added and the fixed effect of `erexp` is dropped.

```{r raneff3}
out2m7 <- lmer(score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + eesex:ersex + (1 + eesex + eesalary|erid), data=dat2, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
```

Next, the cross-level interaction will be examined. There are three combination that can be added. 

1. Model 8 added `eesalary:ersex`. Model 8 was not significantly better than Model 7. Thus, this interaction will be dropped.

```{r crosslevel31}
out2m8 <- lmer(score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + eesex:ersex + eesalary:ersex + (1 + eesex + eesalary|erid), data=dat2, REML=FALSE)
anova(out2m8, out2m7)
```

2. Model 10 added `eesalary:erexp` (and `erexp`) in the model. Model 10 is compared with Model 9, which is Model 7 with `erexp`. The models were not significantly different so `eesalary:erexp` is dropped.

```{r crosslevel32}
out2m9 <- lmer(score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + erexp + eesex:ersex + (1 + eesex + eesalary|erid), data=dat2, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
out2m10 <- lmer(score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + erexp + eesex:ersex + eesalary:erexp + (1 + eesex + eesalary|erid), data=dat2, REML=FALSE)
anova(out2m10, out2m9)
```

3. Model 11 added `eesex:erexp` (and `erexp`) in the model. Model 11 was not significantly better than Model 9, which is Model 7 with `erexp`. Then, `eesalary:erexp` is dropped.

```{r crosslevel33}
out2m11 <- lmer(score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + erexp + eesex:ersex + eesex:erexp + (1 + eesex + eesalary|erid), data=dat2, REML=FALSE)
anova(out2m11, out2m9)
```

All cross-level interactions were not significant. Therefore, Model 7 is the final model.