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

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

## Likelihood Ratio Test

Read the data set and run two models which are different in only the fixed effect of `goldcolor`. The restricted model, `out1r`, does not have the fixed effect of `goldcolor` and the full model, `out1f`, have the fixed effect of `goldcolor`.

```{r lrt1}
dat1 <-read.table("lecture4ex1.csv", sep=",", header=TRUE)
library(lme4)
out1r <- lmer(consume ~ 1 + I(length - 15) + (1|groupid), data=dat1, REML=FALSE)
out1f <- lmer(consume ~ 1 + I(length - 15) + goldcolor + (1|groupid), data=dat1, REML=FALSE)
```

To run the likelihood ratio test, the `anova` function is used. If significant, the full model is preferred. If not significant, the restricted model is preferred.

```{r lrt2}
anova(out1r, out1f)
```

## Random Slope Models

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

Let's start with the random analysis of covariance model (which is a submodel of the random intercept model). `length` and `goldcolor` are both goldfish-level predictors. In this model, the effects of both `length` and `goldcolor` are fixed across tanks. Note that `control = lmerControl(optimizer ="Nelder_Mead")` is used here to help the analysis to get a convergent result.

```{r rsm1a}
library(lme4)
out1m0 <- lmer(consume ~ 1 + I(length - 15) + goldcolor + (1|groupid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
```

Next, the effect of `length` is varied across tanks. `(1 + I(length - 15)|groupid)` means that both `1` and `I(length - 15)` effects are varied across `groupid`. `1` means that the intercepts are varied across `groupid` and `I(length - 15)` means that the slopes of this variable is varied across `groupid`. The `anova` function is used to run likelihood ratio test to check whether allowing random slopes is significant.

```{r rsm1b}
out1m0a <- lmer(consume ~ 1 + I(length - 15) + goldcolor + (1 + I(length - 15)|groupid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1m0, out1m0a)
```

Next, the effect of `goldcolor` is varied across tanks. The `anova` function is used to run likelihood ratio test to check whether allowing random slopes is significant.

```{r rsm1c}
out1m0b <- lmer(consume ~ 1 + I(length - 15) + goldcolor + (1 + goldcolor|groupid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1m0, out1m0b)
```

The random slope of `length` is significantly better than the fixed slope of `length`. However, the random slope of `goldcolor` was not significant. Thus, the model with only random slope of `length` is used.

```{r rsm1d}
summary(out1m0a)
```

Check the scatterplot between random effects in this model by the `ranef` and `plot` functions. The random intercepts (on Y axis) is positively correlated with the random slopes of `length` (on X axis). That is, the stronger the effect of `length` in a tank, the more food fish eats in that tank.

```{r rsm1e}
ranef1m0a <- ranef(out1m0a)
plot(ranef1m0a)
```

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

Again, the random ANCOVA model with two interviewee-level predictors is used as the baseline model to compare with the random slope models. Note that the `iq` is scaled because the model with original-scaled `iq` was not convergent.

```{r rsm2a}
dat2 <-read.table("lecture4ex2.csv", sep=",", header=TRUE)
out2m0 <- lmer(score ~ 1 + eesex + I(scale(iq)) + (1|erid), data=dat2, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
```

Next, the effect of interviewees' sex, `eesex`,  is random. The `anova` function is used to run likelihood ratio test to check whether allowing random slopes is significant.

```{r rsm2b}
out2m0a <- lmer(score ~ 1 + eesex + I(scale(iq)) + (1 + eesex|erid), data=dat2, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out2m0, out2m0a)
```

Next, the effect of interviewees' IQ, `iq`,  is random. The `anova` function is used to run likelihood ratio test to check whether allowing random slopes is significant. 

```{r rsm2c}
out2m0b <- lmer(score ~ 1 + eesex + I(scale(iq)) + (1 + I(scale(iq))|erid), data=dat2, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out2m0, out2m0b)
```

The random slope of interviewees' sex was significant so this random slope is retained. Before checking the results by the `summary` function, analyzing `iq` without random slopes does not require the standard scores of `iq` and the `control` argument. Thus, the original score of `iq` centered at 100 is used and the `control` argument is removed. 

```{r rsm2d}
out2m0a <- lmer(score ~ 1 + eesex + I(iq - 100) + (1 + eesex|erid), data=dat2, REML=FALSE)
summary(out2m0a)
```

Plot the scatterplot between random intercepts and random slopes.

```{r rsm2e}
ranef2m0a <- ranef(out2m0a)
plot(ranef2m0a)
```

From the random slope model, the interviewer-level predictor can be added to predict the ratings. This model is similar to the random intercept model from the previous chapter with the random slope of the effect of interviewees' sex.

```{r rsm2f}
out2m0c <- lmer(score ~ 1 + eesex + I(iq - 100) + ersex + (1 + eesex|erid), data=dat2, REML=FALSE)
summary(out2m0c)
```

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

The random ANCOVA model is used as the baseline model. Two customer-level predictors are used to test the random slope: age and sex.

```{r rsm3a}
dat3 <-read.table("lecture4ex3.csv", sep=",", header=TRUE)
out3m0 <- lmer(sat ~ 1 + I(age - 40) + female + (1|tableid), data=dat3, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
```

Check whether the age effects on customer satisfaction are different across tables.

```{r rsm3b}
out3m0a <- lmer(sat ~ 1 + I(age - 40) + female + (1 + I(age - 40)|tableid), data=dat3, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out3m0, out3m0a)
```

Check whether the sex effects on customer satisfaction are different across tables.

```{r rsm3c}
out3m0b <- lmer(sat ~ 1 + I(age - 40) + female + (1 + female|tableid), data=dat3, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out3m0, out3m0b)
```