---
title: "Intro to Multilevel Model: Lecture 6 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)
out <- readr::read_rds("lecture6saveresults.rds")
ss21 <- out[[1]]
ss22 <- out[[2]]
ss11 <- out[[3]]
ss12 <- out[[4]]
ss31 <- out[[5]]
ss32 <- out[[6]]
```

## Cross-level Interactions

##### Example 1: Interviewees' Rating Scores within Interviewers

Read the data set from Lecture 4.

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

First, check whether the interviewers' sex predict the random intercept.

```{r checkrandomintercept1}
library(lme4)
out2m1 <- lmer(score ~ 1 + eesex + I(iq - 100) + (1 + eesex|erid), data=dat2, REML=FALSE)
out2m1a <- lmer(score ~ 1 + eesex + I(iq - 100) + ersex + (1 + eesex|erid), data=dat2, REML=FALSE)
anova(out2m1, out2m1a)
```

Even though the interviewers' sex effect on interview scores was not significant, the interviewers' sex will be used to predict random slopes of interviewees' sex. For running likelihood ratio test, in the full model, `out2m1b`, interviewers' sex predicts both random intercept and random slope while, in the restricted model, `out2m1a`, interviewers' sex predicts only random intercept.

```{r checkrandomslope1}
out2m1b <- lmer(score ~ 1 + eesex + I(iq - 100) + ersex + eesex:ersex 
                + (1 + eesex|erid), data=dat2, REML=FALSE)
anova(out2m1a, out2m1b)
```

Run the `summary` function, the fixed effect of interaction was significant.

```{r summaryrandomslope1}
summary(out2m1b)
```

The `effects` package is used to calcuate the expected means of interview ratings in each combination of interviewer's and interviewee's sex. 

```{r effectspackage1}
library(effects)
eff2m1b <- effect(term="eesex:ersex", mod=out2m1b, 
                         xlevels=list(ersex=c(0, 1), eesex=c(0, 1)))
summary(eff2m1b)
```

The expected means are used to make a clustered bar graph using the `ggplot2` package.

```{r bargraph1}
library(ggplot2)
mydat <- as.data.frame(eff2m1b)
mydat$eesex <- factor(mydat$eesex, labels=c("Male", "Female"))
mydat$ersex <- factor(mydat$ersex, labels=c("Male", "Female"))
ggplot(mydat, aes(factor(eesex), fit, fill = ersex)) + 
  geom_bar(stat="identity", position = "dodge") + 
  labs(x = "Interviewee", y = "Interview Ratings", fill = "Interviewer")
```

Use the `interactions` package to calculate simple slopes. First, check the effect of interviewees' sex at each sex of interviewers.

```{r simpleslope11, eval=FALSE}
library(interactions)
ss21 <- sim_slopes(model=out2m1b, pred=eesex, modx=ersex)
```
```{r simpleslope111}
ss21
```


The `interactions` package can visualize the interaction as well. However, it provides the line plot which might not be appropriate when either the predictor or the moderator is categorical. 

```{r interactionplot11}
interact_plot(model=out2m1b, pred=eesex, modx=ersex)
```

Next, check the simple slope effect of interviewers' sex at each sex of interviewees.

```{r simpleslope12, eval=FALSE}
ss22 <- sim_slopes(model=out2m1b, pred=ersex, modx=eesex)
```
```{r simpleslope122}
ss22
```
  
The interaction can be plotted as well.

```{r interactionplot12}
interact_plot(model=out2m1b, pred=ersex, modx=eesex)
```

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

Read the data set from Lecture 4.

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

First, check whether the tank volume and the number of goldfishes predict the random intercept.

```{r checkrandomintercept2}
out1m1 <- lmer(consume ~ 1 + I(length - 15) + goldcolor + (1 + I(length - 15)|groupid), 
                data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
out1m1intcept <- lmer(consume ~ 1 + I(length - 15) + goldcolor + I(volume - 1) + I(nfish - 8) 
                      + (1 + I(length - 15)|groupid), data=dat1, REML=FALSE, 
                      control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1m1, out1m1intcept)
```

Use the tank volume and the number of goldfishes to predict random slopes of the fish lengths. For running likelihood ratio test, in the full model, `out1m1slope`, the tank-level predictors predict both random intercept and random slope while, in the restricted model, `out1m1intcept`, the tank-level predictors predict only random intercept.

```{r checkrandomslope2}
out1m1slope <- lmer(consume ~ 1 + I(length - 15) + goldcolor + I(volume - 1) + I(nfish - 8) 
                    + I(length - 15):I(volume - 1) + I(length - 15):I(nfish - 8) 
                    + (1 + I(length - 15)|groupid), 
                    data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1m1intcept, out1m1slope)
```

Even though the likelihood ratio test was not significant, let's check each interaction term. 

```{r summaryrandomslope2}
summary(out1m1slope)
```

The the tank volume alone significantly predicted the random slope but the number of fishes could not significantly predict the random slope. Then, the tank volume effect was retained and the effect of the number of fishes was dropped. 

```{r effectspackage2}
out1m1trimmed <- lmer(consume ~ 1 + I(length - 15) + goldcolor + I(volume - 1) + I(nfish - 8) 
                      + I(volume - 1):I(length - 15) + (1 + I(length - 15)|groupid), 
                      data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1m1intcept, out1m1trimmed)
```

Check the output of the cross-level interaction model.

```{r summaryrandomslope22}
summary(out1m1trimmed)
```

Next, find the values of fish lengthes to probe the interaction. The values of 9, 14, and 19 are used because they are close to *M - SD*, *M*, and *M + SD*.

```{r probevalue21}
mean(dat1$length) - sd(dat1$length)
mean(dat1$length)
mean(dat1$length) + sd(dat1$length)
lengthval <- c(9, 14, 19)
```

Using the same logic, the tank volume values of 0.5, 0.8, and 1.1 are used.

```{r probevalue22}
mean(dat1$volume) - sd(dat1$volume)
mean(dat1$volume)
mean(dat1$volume) + sd(dat1$volume)
volumeval <- c(0.5, 0.8, 1.1)
```

Because the `interactions` package does not allow predictors or moderators in `I()`, users need to transform the values in the data frame or simply use the original scale. In this case, the original scales of `length` and `volume` are used. Note that `nfish` is still centered to make the expected means meaningful.

```{r newoutput2}
o1m1 <- lmer(consume ~ 1 + goldcolor + I(nfish - 8) 
                      + length*volume + (1 + length|groupid), 
                      data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
```

Use the `sim_slopes` package to calculate simple slopes. First, check the effect of fish tank volumes at each level of fish lengths.

```{r simpleslope21, eval=FALSE}
ss11 <- sim_slopes(model=o1m1, pred=volume, modx=length, modx.values = lengthval)
```
```{r simpleslope211}
ss11
```

Check the interaction plot. 

```{r interactionplot21}
interact_plot(model=o1m1, pred=volume, modx=length, modx.values = lengthval)
```

Next, check the simple slope effect of fish lengths at each level of fish tank volume.

```{r simpleslope22, eval=FALSE}
ss12 <- sim_slopes(model=o1m1, pred=length, modx=volume, modx.values = volumeval)
```
```{r simpleslope222}
ss12
```

The interaction can be plotted as well.

```{r interactionplot22}
interact_plot(model=out2m1b, pred=ersex, modx=eesex)
```

## Other Interaction Types

##### Example 3: Customers Satisfaction within Each Table in Dining Services

Read the data set from Lecture 4.

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

Run the model with the interaction between age and sex.

```{r interactionoutput3}
out3m1 <- lmer(sat ~ 1 + I(age - 40)*female + I(numperson-4) + outdoor 
               + (1|tableid), data=dat3, REML=FALSE)
summary(out3m1)
```

The interaction was significant. Next, the `interactions` package is used. Remove `I()` and rerun the analysis.

```{r rerunoutput3}
out3m1a <- lmer(sat ~ 1 + age*female + I(numperson-4) + outdoor 
               + (1|tableid), data=dat3, REML=FALSE)
```

Use the `sim_slopes` to calculate simple slopes. First, check the effect of age at each sex.

```{r simpleslope31, eval=FALSE}
ss31 <- sim_slopes(model=out3m1a, pred=age, modx=female, modx.values = c(0, 1))
```

```{r simpleslope312}
ss31
```

Check the interaction plot. 

```{r interactionplot31}
interact_plot(model=out3m1a, pred=age, modx=female, modx.values = c(0, 1))
```

Next, check the sex differences at each level of age. Before probing the interaction, find the appropriate values of age to probe. Let's try the `quantile` function to get 25^th^, 50^th^, and 75^th^ percentiles. The values of 25, 35, and 45 are picked because they are close to the `quantile` result. 

```{r probevalue31}
quantile(dat3$age, c(0.25, 0.50, 0.75))
ageval <- c(25, 35, 45)
```

Next, check the simple effect of sex difference at each level of age.

```{r simpleslope32, eval=FALSE}
ss32 <- sim_slopes(model=out3m1a, pred=female, modx=age, modx.values = ageval)
```

```{r simpleslope322}
ss32
```

The interaction can be plotted as well.

```{r interactionplot32}
interact_plot(model=out3m1a, pred=female, modx=age, modx.values = ageval)
```

Let's make the clustered bar graphs. First, the expected means of males and females at the ages of 25, 35, and 45 are calculated. First, the simple intercepts and slopes are needed.

```{r outputsimslopes32}
ss32$ints
ss32$slopes
```

The simple intercepts are the expected male satisfaction at each age. The sums of simple intercepts and simple slopes are the expected female satisfaction at each age. Next, a data frame containing expected means is created.

```{r dataframe3}
mydat <- data.frame(sex=c(rep("Male", 3), rep("Female", 3)), 
                    age=rep(ageval, 2), 
                    sat=c(ss32$ints[,"Est."], ss32$ints[,"Est."] + ss32$slopes[,"Est."]))
mydat
```

Similar to the first example, the clustered bar graph is created by the `ggplot2` package.

```{r ggplot3}
library(ggplot2)
ggplot(mydat, aes(factor(age), sat, fill = sex)) + 
  geom_bar(stat="identity", position = "dodge") + 
  labs(x = "Age", y = "Customer Satisfaction", fill = "Sex")
```