---
title: "Intro to Multilevel Model: Lecture 7 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)
out <- readr::read_rds("lecture7saveresults.rds")
ss41 <- out[[1]]
ss42 <- out[[2]]
ss31 <- out[[3]]
```

## No-Centering vs. Group-Mean-Centering

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

Read the data set from Lecture 4.

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

Calculate group means by the `ave` function. The result has the same length as the original score and provides the group means that each data point belongs to. The group-mean centered variable is the difference between the original score and the group mean.

```{r calcgroupmean1}
dat3$aveage <- ave(dat3$age, dat3$tableid)
dat3$diffage <- dat3$age - dat3$aveage
```

Calculate the output when `age` is not centered. The regression coefficient of `age` equals both the between-table and within-table effects.

```{r nocentering1}
library(lme4)
out3m1 <- lmer(sat ~ 1 + age + (1|tableid), data=dat3, REML=FALSE)
summary(out3m1)
```

Check when `age` is group-mean centered. The regression coefficient represents the within-table effect.

```{r groupcenter1}
out3m2 <- lmer(sat ~ 1 + diffage + (1|tableid), data=dat3, REML=FALSE)
summary(out3m2)
```

Let's predict the customer satisfaction with group means of `age` too. Check the result when both the original scale of `age` and the group means of `age` are the predictors. The regression coefficient of the original scale of `age` is the within-table effect. However, the regression coefficient of the `age` table averages is the difference between within- and between-table effects.

```{r nocenterwithmean1}
out3m3 <- lmer(sat ~ 1 + age + aveage + (1|tableid), data=dat3, REML=FALSE)
summary(out3m3)
```

Check the result when both the group-centered `age` and the group means of `age` are the predictors. The regression coefficient of the group-centered `age` is the within-table effect and the regression coefficient of the `age` table averages is the between-table effect.

```{r groupcenterwithmean1}
out3m4 <- lmer(sat ~ 1 + diffage + aveage + (1|tableid), data=dat3, REML=FALSE)
summary(out3m4)
```

Before trying to find the cross-level interaction of `age`, the `age` needs to be rescaled to make the multilevel analysis convergent. The scaled age, `sage`, is the original scale of `age` divided by 10.

```{r rescale1}
dat3$sage <- dat3$age / 10
dat3$avesage <- ave(dat3$sage, dat3$tableid)
dat3$diffsage <- dat3$sage - dat3$avesage
```

Let's rerun the previous analyses. 

```{r rerunwithscaleage1}
out3m1 <- lmer(sat ~ 1 + sage + (1|tableid), data=dat3, REML=FALSE)
summary(out3m1)

out3m2 <- lmer(sat ~ 1 + diffsage + (1|tableid), data=dat3, REML=FALSE)
summary(out3m2)

out3m3 <- lmer(sat ~ 1 + sage + avesage + (1|tableid), 
               data=dat3, REML=FALSE)
summary(out3m3)

out3m4 <- lmer(sat ~ 1 + diffsage + avesage + (1|tableid), 
               data=dat3, REML=FALSE)
summary(out3m4)
```

The fifth model has `sage` (the original scale divided by 10), `avesage` (the group means), and their interaction. As mentioned in the lecture, the regression coefficients are not intuitive. The effect of `sage` represents the within-table effect at `avesage` of 0. The effect of `avesage` represents the difference between the within-table and between-table effect at `avesage` of 0. The interaction effect represents the increase in both within-table and between-table effects. 

```{r nocenteringinteraction1}
out3m5 <- lmer(sat ~ 1 + sage + avesage + sage:avesage 
               + (1 + sage|tableid), data=dat3, REML=FALSE)
summary(out3m5)
```

The sixth model has `diffsage` (the group-centered age), `avesage` (the group means), and their interaction. As mentioned in the lecture, the regression coefficients are easier to understand. The effect of `diffsage` represents the within-table effect at `avesage` of 0. The effect of `avesage` represents the between-table effect. The interaction effect represents the increase in within-table effects as `avesage` increased by 1. 

```{r groupcenteringinteraction1}
out3m6 <- lmer(sat ~ 1 + diffsage + avesage + diffsage:avesage 
               + (1 + diffsage|tableid), data=dat3, REML=FALSE)
summary(out3m6)
```

## Probing Cross-Level Interactions between the Same Predictor

##### Example 2: The Effect of Math Achievement on Perceived Self-Efficacy

Read the new data set and check the variables within data.

```{r readdata2}
dat4 <- read.table("lecture7ex1.csv", sep=",", header=TRUE)
library(psych)
describe(dat4)
```

Calculate the group means and group-mean-centered math achievement scores.

```{r calcgroupmean2}
dat4$aveach <- ave(dat4$ach, dat4$schoolid)
dat4$diffach <- dat4$ach - dat4$aveach
```

Predict self-efficacy with the group-mean-centered achievement, `diffach`, group means of achievement, `aveach`, its interaction, and whether the school is private, `private`.

```{r crosslevel2}
out4m1a <- lmer(efficacy ~ 1 + diffach + aveach + private
               + diffach:aveach
               + (1 + diffach|schoolid), data=dat4, REML=FALSE)
summary(out4m1a)
```

Find values for probing interactions. Let's check the within-school and between-school standard deviations of math achievement.

```{r nullmodel2}
out4m0 <- lmer(ach ~ 1 + (1|schoolid), data=dat4, REML=FALSE)
summary(out4m0)
```

The result of the null model showed that the school-level standard deviation was 9.637 and student-level standard deviation was 9.745. The grand mean of math achievement was 54.215. Thus, `aveach` values to examine the within-group effect is 55, $55 - 10 = 45$, $55 + 10 = 65$. Thus, `diffach` values to examine the between-group effect is -10, 0, 10.

```{r modvalues2}
aveachval <- c(45, 55, 65)
diffachval <- c(-10, 0, 10)
```

Use the `interactions` package to calculate simple slopes. First, check the within-school effect of math achievement at each level of schools' math achievement averages. 

```{r simpleslope21, eval=FALSE}
library(interactions)
ss41 <- sim_slopes(model=out4m1a, pred=diffach, modx=aveach, modx.values=aveachval)
```
```{r simpleslope211}
ss41
```

Use the `interactions` package can visualize the interaction. 

```{r interactionplot21}
interact_plot(model=out4m1a, pred=diffach, modx=aveach, modx.values=aveachval)
```

Next, check the simple slope of between-school effect of math achievement at each level of within-school deviation of math achievement. 

```{r simpleslope22, eval=FALSE}
ss42 <- sim_slopes(model=out4m1a, pred=aveach, modx=diffach, modx.values=diffachval)
```
```{r simpleslope222}
ss42
```
  
The interaction can be plotted as well.

```{r interactionplot22}
interact_plot(model=out4m1a, pred=aveach, modx=diffach, modx.values=diffachval)
```

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

Let's examine another example of cross-level interactions. Instead of `age`, the cross-level interaction of `female` is used. Calculate the group means and group-mean-centered `female`.

```{r calcgroupmean3}
dat3$avefemale <- ave(dat3$female, dat3$tableid)
dat3$difffemale <- dat3$female - dat3$avefemale
```

First, check whether the effect of `female` is varied across tables.

```{r randomslope3}
out3g1 <- lmer(sat ~ 1 + female + avefemale 
               + (1|tableid), data=dat3, REML=FALSE)
out3g2 <- lmer(sat ~ 1 + female + avefemale 
               + (1 + female|tableid), data=dat3, REML=FALSE)
anova(out3g1, out3g2)
```

Next, check whether the table proportions of females can predict the varying effect of `female`.

```{r crosslevel3}
out3g3 <- lmer(sat ~ 1 + female + avefemale + female:avefemale 
                + (1 + female|tableid), data=dat3, REML=FALSE)
anova(out3g2, out3g3)
```

The interaction effect was significant. Check the output.

```{r crosslevel32}
summary(out3g3)
```

As explained in the lecture, the between-table effect of females is actually curvilinear. Let's check the graph.

```{r curveplot3}
avefemale <- seq(0, 1, 0.01)
est <- coef(summary(out3g3))[,"Estimate"]
cust <- est[1] + ((est[2] + est[3]) * avefemale) + (est[4]*avefemale*avefemale)
plot(avefemale, cust, type="l", xlim=c(0, 1), ylim=c(60, 70), 
     xlab="Female Proportion in a Table", ylab="Customer Satisfaction")
```

The values of 20%, 50%, and 80% females in each table are used to probe the interaction.

```{r modvalues3}
avefemaleval <- c(0.2, 0.5, 0.8)
```

Check the simple slopes of the within-table effect of females at each level of the table proportions of females. 

```{r simpleslope31, eval=FALSE}
ss31 <- sim_slopes(model=out3g3, pred=female, modx=avefemale, modx.values=avefemaleval)
```
```{r simpleslope311}
ss31
```

Check the interaction plot.

```{r interactionplot31}
interact_plot(model=out3g3, pred=female, modx=avefemale, modx.values=avefemaleval)
```

Please not that the effect of the interaction of no-centered females and the table proportions of females is not only represent the moderating within-table effect but also the between-table effect too. The model is equally constrained the within- and between-table moderating effects. 

Thus, this model is not easy to understand. The better model in this case is to use group-mean-centered females. First, check the random intercept model.

```{r groupcenterrandomintercept3}
out3f1 <- lmer(sat ~ 1 + difffemale + avefemale + (1|tableid), 
               data=dat3, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))
summary(out3f1)
```

Next, check whether the within-table effect of females is varied across tables.

```{r groupcenterrandomslope3}
out3f2 <- lmer(sat ~ 1 + difffemale + avefemale 
               + (1 + difffemale|tableid), data=dat3, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))
anova(out3f1, out3f2)
```

Next, check whether the cross-level interaction is significant.

```{r groupcentercrosslevel3}
out3f3 <- lmer(sat ~ 1 + difffemale + avefemale + difffemale:avefemale 
               + (1 + difffemale|tableid), data=dat3, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))
anova(out3f2, out3f3)
```

The cross-level interaction was not significant. Thus, the variation of within-table effect of females between tables was not significantly explained by the proportion of females in each table. However, the cross-level interaction was significant in the previous analysis using no-centered females. It is possible that the interaction between levels was not equal. The within-table effect did not have interaction whereas the between-table effect has the interaction. Thus, add the quadratic term of tables' female proportion in the random slope model.

```{r groupcenterquadratic3}
out3f4 <- lmer(sat ~ 1 + difffemale + avefemale + I(avefemale^2)
               + (1 + difffemale|tableid), data=dat3, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))
anova(out3f2, out3f4)
```

The quadratic term was significant. Make sure whether the quadratic term of tables' female proportion predict the random slopes of the within-table effect of females.

```{r groupcenterquadraticinteraction3}
out3f5 <- lmer(sat ~ 1 + difffemale + avefemale + I(avefemale^2)
               + avefemale:difffemale + I(avefemale^2):difffemale
               + (1 + difffemale|tableid), data=dat3, REML=FALSE, 
               control = lmerControl(optimizer ="Nelder_Mead"))
anova(out3f4, out3f5)
```

The interaction effect was not significant. Let's check the output of the random slope model with quadratic terms.

```{r finalfemaleeffect3}
summary(out3f4)
```

Let's visualize the quadratic effect of the tables' female proportion.

```{r quadraticplot3}
est <- coef(summary(out3f4))
avefemale <- seq(0, 1, 0.01)
cust <- est[1] + (est[3] * avefemale) + (est[4]*avefemale*avefemale)
plot(avefemale, cust, type="l", xlim=c(0, 1), ylim=c(60, 70), 
     xlab="Female Proportion in a Table", ylab="Customer Satisfaction")

```