---
title: "Intro to Multilevel Model: Lecture 13 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(nlme)
library(psych)
library(interactions)
library(lavaan)
library(MASS)
library(r2mlm)
```

## Three-level Null Model

##### Example 1: Math Achievement Data across Schools and Countries 

Read the new data set where students' math achievement and IQ were collected across schools and countries. The nesting structure has three levels; that is, students:schools:countries.  

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

See descriptive statistics of all variables within the data set.

```{r describe1}
library(psych)
describe(dat1)
```

Let's start with the null model which `math` is the dependent variable. Note that `(1|schoolid)` and `(1|countryid)` means that the random intercepts are varied across schools and countries, respectively.  

```{r nullmodel11}
m10 <- lmer(math ~ 1 + (1|schoolid) + (1|countryid), data=dat1, REML=FALSE)
summary(m10)
```

Next, check the null model where `iq` is the dependent variable.

```{r nullmodel12}
m11 <- lmer(iq ~ 1 + (1|schoolid) + (1|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
summary(m11)
```

For the following analysis, `iq` is subtracted by 100 and divided by 15, which are the mean and standard deviation of the standard IQ score. The transformation will help model easier to converge.

```{r fixedtime1}
dat1$iqc <- (dat1$iq - 100)/15
```

## Centering

##### Example 1: Math Achievement Data across Schools and Countries 

Centering in three-level models is quite complicated. In this example, two predictors will be centered: `iqc` and `private`. `iqc` is a level-1 predictor. First, make the school means of IQ and deviation scores of IQ from school means. This process is similar to those in two-level models.

```{r level2center1}
dat1$aveschooliqc <- ave(dat1$iqc, dat1$schoolid)
dat1$diffschooliqc <- dat1$iqc - dat1$aveschooliqc
```

From the school means of IQ, the country means of IQ and school-level deviation scores of IQ from country means will be calculated. To find the country means, the school-level dataset is calculated.

```{r schooldata1}
dat1school <- dat1[!duplicated(dat1$schoolid),]
```

Next, calculate the country means of IQ and school-level deviation scores of IQ from country means.

```{r level3center11}
dat1school$avecountryiqc <- ave(dat1school$aveschooliqc, dat1school$countryid)
dat1school$diffcountryiqc <- dat1school$aveschooliqc - dat1school$avecountryiqc
```

`private` is a level-2 predictor. Thus, the centering must be implemented in the school-level data.

```{r level3center12}
dat1school$aveprivate <- ave(dat1school$private, dat1school$countryid)
dat1school$diffprivate <- dat1school$private - dat1school$aveprivate
```

Four new transformed variables in the school-level data will be attached in the original data. First, the school ID between both data are matched. The result will show that the school in each row of the original data is in which rows in the school-level data.

```{r matchstudentschooldata1}
matchit <- match(dat1$schoolid, dat1school$schoolid)
```

Then, the result of the `match` function is used to attach transformed variables from the school-level data to the student-level (original) data.

```{r attachdata1}
dat1$avecountryiqc <- dat1school[matchit, "avecountryiqc"]
dat1$diffcountryiqc <- dat1school[matchit, "diffcountryiqc"]
dat1$aveprivate <- dat1school[matchit, "aveprivate"]
dat1$diffprivate <- dat1school[matchit, "diffprivate"]
```

After this, many models centered and uncentered variables are shown to illustrate their differences.

First, only deviation-from-school IQ is used in the model. The slope represents the student-level effect (controlling schools) of IQ on math achievement.

```{r center11}
m1xa1 <- lmer(math ~ 1 + diffschooliqc + (1|schoolid) + (1|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
summary(m1xa1)
```

Next, the predictor is uncentered IQ. The slope represents all student-level (controlling schools), school-level (controlling countries), and country-levle effects. The slopes in three levels are equally constrained. 

```{r uncenter11}
m1xa2 <- lmer(math ~ 1 + iqc + (1|schoolid) + (1|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
summary(m1xa2)
```

Next, deviation-from-school IQ, deviation-from-country IQ, and country IQ means are the predictors. The slope of deviation-from-school IQ represents the student-level effect (controlling schools) of IQ on math achievement. The slope of deviation-from-country IQ represents the school-level effect (controlling countries) of IQ on math achievement. The slope of country IQ means represents the country-level effect of IQ on math achievement.

```{r center12}
m1xb1 <- lmer(math ~ 1 + diffschooliqc + diffcountryiqc + avecountryiqc + (1|schoolid) + (1|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
summary(m1xb1)
```

Next, deviation-from-school IQ, school IQ means, and country IQ means are the predictors. The slope of deviation-from-school IQ represents the student-level effect (controlling schools) of IQ on math achievement. The slope of school IQ means still represents the school-level effect (controlling countries) of IQ on math achievement. However, the slope of country IQ means represents the country-level effect subtracted by the school-level effect.

```{r center13}
m1xb2 <- lmer(math ~ 1 + diffschooliqc + aveschooliqc + avecountryiqc + (1|schoolid) + (1|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
summary(m1xb2)
```

Next, student IQ scores, school IQ means, and country IQ means are the predictors. The slope of student IQ scores represents the student-level effect (controlling schools) of IQ on math achievement. The slope of school IQ means represents the school-level effect subtracted by the student-level effect. The slope of country IQ means represents the country-level effect subtracted by the school-level effect.

```{r uncenter12}
m1xb3 <- lmer(math ~ 1 + iqc + aveschooliqc + avecountryiqc + (1|schoolid) + (1|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
summary(m1xb3)
```

Let's illustrate the impact of centering on level-2 predictors. Initially, the deviation-from-country calculated from `private` is used. The slope represents the school-level effect (controlling countries) of type of schools on math achievement.

```{r center2level11}
m1wa1 <- lmer(math ~ 1 + diffprivate + (1|schoolid) + (1|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
summary(m1wa1)
```

Next, the original `private` variable is the predictor. The slope represents both school-level effect (controlling countries) and country-level effect of type of schools on math achievement. The effects on both levels are equally constrained.

```{r uncenter2level11}
m1wa2 <- lmer(math ~ 1 + private + (1|schoolid) + (1|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
summary(m1wa2)
```

Next, the deviation-from-country calculated from `private` and the country means of `private` are used as predictors. The first slope represents the school-level effect (controlling countries) of type of schools on math achievement. The second slope represents the country-level effect.

```{r center2level12}
m1wb1 <- lmer(math ~ 1 + diffprivate + aveprivate + (1|schoolid) + (1|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
summary(m1wb1)
```

Next, the original `private` variable and the country means of `private` are the predictors. The first slope represents both school-level effect (controlling countries) and country-level effect of type of schools on math achievement. The second slope represent the country-level effect subtracted by the school-level effect.

```{r uncenter2level12}
m1wb2 <- lmer(math ~ 1 + private + aveprivate + (1|schoolid) + (1|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
summary(m1wb2)
```

Let's combine all school-mean and country-mean centered variables. Also, add country's educational quality and country's educational opportunity indices in the model. The full random intercept model is as follows:

```{r randomintercept1}
dat1country <- dat1[!duplicated(dat1$countryid),]
apply(dat1country, 2, mean)

m13 <- lmer(math ~ 1 + diffschooliqc + diffcountryiqc + I(avecountryiqc+0.691) + diffprivate + I(aveprivate-0.506) + I(quality-53.33) + I(opportunity-47.98) + (1|schoolid) + (1|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
summary(m13)
```

## Random Slopes and Cross-Level Interactions

##### Example 1: Math Achievement Data across Schools and Countries 

Next, let's check whether the **student-level** effect (`diffschooliqc`) on math is random across **schools**. See that the only differences between two models in `anova` is `(1 + diffschooliqc|schoolid)`. The effect was not significant.

```{r randomslope11}
m14 <- lmer(math ~ 1 + diffschooliqc + diffcountryiqc + I(avecountryiqc+0.691) + diffprivate + I(aveprivate-0.506) + I(quality-53.33) + I(opportunity-47.98) + (1 + diffschooliqc|schoolid) + (1|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead", calc.derivs=FALSE))
anova(m13, m14)
```

Next, check whether the **student-level** effect (`diffschooliqc`) on math is random across **countries**. See that the only differences between two models in `anova` is `(1 + diffschooliqc|countryid)`. The effect was significant.

```{r randomslope12}
m15 <- lmer(math ~ 1 + diffschooliqc + diffcountryiqc + I(avecountryiqc+0.691) + diffprivate + I(aveprivate-0.506) + I(quality-53.33) + I(opportunity-47.98) + (1|schoolid) + (1 + diffschooliqc|countryid), data=dat1, REML=FALSE)
anova(m13, m15)
```

Next, check whether the **school-level** effect (`diffcountryiqc`) on math is random across **countries**. The last model, `m15`, is used as the restricted model. The full model added the random effect of `diffcountryiqc` across countries, which shown in `(1 + diffschooliqc + diffcountryiqc|countryid)`. The effect was not significant. Then, `m15` is retained.

```{r randomslope13}
m15 <- lmer(math ~ 1 + diffschooliqc + diffcountryiqc + I(avecountryiqc+0.691) + diffprivate + I(aveprivate-0.506) + I(quality-53.33) + I(opportunity-47.98) + (1|schoolid) + (1 + diffschooliqc|countryid), data=dat1, REML=FALSE)
anova(m13, m15)
```

Next, check whether the **school-level** effect of type of schools (`diffprivate`) on math is random across **countries**. `m15` is used as the restricted model. The full model added the random effect of `diffprivate` across countries, which shown in `(1 + diffschooliqc + diffprivate|countryid)`. The effect was significant. Then, the full model, `m17`, is retained.

```{r randomslope14}
m17 <- lmer(math ~ 1 + diffschooliqc + diffcountryiqc + I(avecountryiqc+0.691) + diffprivate + I(aveprivate-0.506) + I(quality-53.33) + I(opportunity-47.98) + (1|schoolid) + (1 + diffschooliqc + diffprivate|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(m15, m17)
```

Check the output of the final random-slope model.

```{r randomslope15}
summary(m17)
```

Only one random effect is in the school level, which is the random intercept. However, three random effects are in the country level: (a) random intercepts, (b) the effect of student-level IQ, and (c) the effect of school-level type of schools. The correlation between random effects are plotted here. Note that the random intercepts are the adjusted means of math achievement when the education quality and opportunity of a country are equal to their grand means.

```{r randomeffectplot1}
ranefm17 <- ranef(m17)
names(ranefm17)
ranefm17country <- ranefm17$countryid
plot(ranefm17country)
```

Because the random slopes are only in the country level. Thus, four country-level predictors are used to predict both random slopes. As a result, there are eight cross-level interactions.

```{r crosslevel11, warn=TRUE}
m18 <- lmer(math ~ 1 + diffschooliqc + diffcountryiqc + I(avecountryiqc+0.691) + diffprivate + I(aveprivate-0.506) + I(quality-53.33) + I(opportunity-47.98) + diffschooliqc*I(avecountryiqc+0.691) + diffschooliqc*I(aveprivate-0.506) + diffschooliqc*I(quality-53.33) + diffschooliqc*I(opportunity-47.98) + diffprivate*I(avecountryiqc+0.691) + diffprivate*I(aveprivate-0.506) + diffprivate*I(quality-53.33) + diffprivate*I(opportunity-47.98) + (1|schoolid) + (1 + diffschooliqc + diffprivate|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
summary(m18)
```

Unfortunately, the warning, `boundary (singular) fit`, is shown. It means that at least a pair of random effects have a perfect correlation. The result is not trustworthy. Usually, a random effect should be dropped. However, in this case, the random effect of student-level IQ is predicted by two school-level predictors. The perfect random effect is not shown even though the two added cross-level interactions were not significant. Thus, this model is used for probing interactions.

```{r crosslevel12}
m19 <- lmer(math ~ 1 + diffschooliqc + diffcountryiqc + I(avecountryiqc+0.691) + diffprivate + I(aveprivate-0.506) + I(quality-53.33) + I(opportunity-47.98) + diffschooliqc*diffcountryiqc + diffschooliqc*I(avecountryiqc+0.691) + diffschooliqc*diffprivate + diffschooliqc*I(aveprivate-0.506) + diffschooliqc*I(quality-53.33) + diffschooliqc*I(opportunity-47.98) + diffprivate*I(avecountryiqc+0.691) + diffprivate*I(aveprivate-0.506) + diffprivate*I(quality-53.33) + diffprivate*I(opportunity-47.98) + (1|schoolid) + (1 + diffschooliqc + diffprivate|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
summary(m19)
```

Before probing interactions, analyze the new model without `I()`.

```{r crosslevel13}
dat1$avecountryiqcc <- dat1$avecountryiqc + 0.691
dat1$aveprivatec <- dat1$aveprivate - 0.506
dat1$qualityc <- dat1$quality - 53.33
dat1$opportunityc <- dat1$opportunity - 47.98

m20 <- lmer(math ~ 1 + diffschooliqc + diffcountryiqc + avecountryiqcc + diffprivate + aveprivatec + qualityc + opportunityc + diffschooliqc*diffcountryiqc + diffschooliqc*avecountryiqcc + diffschooliqc*diffprivate + diffschooliqc*aveprivatec + diffschooliqc*qualityc + diffschooliqc*opportunityc + diffprivate*avecountryiqcc + diffprivate*aveprivatec + diffprivate*qualityc + diffprivate*opportunityc + (1|schoolid) + (1 + diffschooliqc + diffprivate|countryid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
```

Let's pick values for probing interactions. The *M*, *M - SD*, and *M + SD* are used for`diffschooliqc` (using student-level *SD*), `qualityc`, and `opportunityc` (using country-level *SD*s). 

```{r probevalue1}
mdiq <- mean(dat1$diffschooliqc)
sddiq <- sd(dat1$diffschooliqc)

mq <- mean(dat1country$quality)
sdq <- sd(dat1country$quality)

mo <- mean(dat1country$opportunity)
sdo <- sd(dat1country$opportunity)

diffschooliqcval <- c(mdiq - sddiq, mdiq, mdiq + sddiq)
qualitycval <- c(mq - sdq, mq, mq + sdq) - 53.33
opportunitycval <- c(mo - sdo, mo, mo + sdo) - 47.98
```

The `interactions` package is usually used for probing interactions. However, the `sim_slopes` function took very long time to run the analysis (longer than 12 hours). I wrote a function to analyze simple slopes with much shorter time. To be honest, I am still not sure why the `sim_slopes` function took so long. To use my function, the `quicksimslopeslme4.R` file must be placed in the working directory. Then, run the following function to read the function in the file.

```{r readfunction1}
source("quicksimslopeslme4.R")
```

My function is called `quick_sim_slopes`. The features are not as many as the `sim_slopes` function but it saved a lot of time. First, investigate the interaction between the student-level effect of IQ and the country's educational quality. Check the slopes of the student-level effect of IQ at each level of country's educational quality.

```{r simslopes11}
quick_sim_slopes(model=m20, pred="diffschooliqc", modx="qualityc", modx.values=qualitycval)
```

Plot the student-level effect of IQ at each level of country's educational quality using the `interactions` package.

```{r plotslopes11}
library(interactions)
interact_plot(model=m20, pred="diffschooliqc", modx="qualityc", modx.values=qualitycval)
```

Next, check the slopes of the country's educational quality at each level of student's IQ deviation from school means.

```{r simslopes12}
quick_sim_slopes(model=m20, pred="qualityc", modx="diffschooliqc", modx.values=diffschooliqcval)
```

Plot the the country's educational quality effect at each level of student's IQ deviation from school means.

```{r plotslopes12}
interact_plot(model=m20, pred="qualityc", modx="diffschooliqc", modx.values=diffschooliqcval)
```

Next, investigate the interaction between the school-level effect of school types and the country's educational opportunity. Check the slopes of the school-level effect of school type at each level of country's educational opportunity.

```{r simslopes13}
quick_sim_slopes(model=m20, pred="diffprivate", modx="opportunityc", modx.values=opportunitycval)
```

Plot the school-level effect of school type at each level of country's educational opportunity.

```{r plotslopes13}
interact_plot(model=m20, pred="diffprivate", modx="opportunityc", modx.values=opportunitycval)
```

Note that users can run the effect of `opportunityc` at each level of `diffprivate` but it is very hard to interpret. So this setup will be skipped. 

## Three-level models with longitudinal analysis

##### Example 2: Change in schools' standardized scores

Read the new data set where the standardized test data were collected from the same schools three years in a row. However, students who took the tests in each year were different across time points. Thus, the nesting structure has three levels; that is, students:years:schools.  

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

See descriptive statistics of all variables within the data set.

```{r describe2}
describe(dat2)
```

Center the year such that 2017, 2018, and 2019 are 0, 1, 2, respectively.

```{r center2}
dat2$yearc <- dat2$year - 2017
```

Next, check whether schools can significantly explain variance in standardized tests. The random intercept across schools was significant.

```{r schooltest2}
m21 <- lmer(test ~ 1 + yearc + (1|timeid), data=dat2, REML=FALSE)
m22 <- lmer(test ~ 1 + yearc + (1|timeid) + (1|schoolid), data=dat2, REML=FALSE)
anova(m21, m22)
```

Check the changes in standardized tests from 2017 to 2019 where the slopes are equal across schools.  

```{r nullmodel21}
summary(m22)
```

Next, allow the rates of changes to be different across schools. The random slope was significant.

```{r randomslope21}
m23 <- lmer(test ~ 1 + yearc + (1|timeid) + (1 + yearc|schoolid), data=dat2, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(m22, m23)
```

See the output of the random slope model.

```{r randomslope22}
summary(m23)
```

Make the scatterplot to visualize the relationship between the schools' averaged standardized tests scores in 2017 and the rate of changes in each school.

```{r plotrandomslope2}
plot(ranef(m23)$schoolid, xlab="School Averaged Test Score in 2017", ylab="School Change in Test Score per Year")
```

Next, check whether both random intercepts and random slopes can be explained by types of schools (public vs. private).

```{r crosslevel2}
m24 <- lmer(test ~ 1 + yearc*public + (1|timeid) + (1 + yearc|schoolid), data=dat2, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
summary(m24)
```

The interaction was significant. Set the values to probe the interaction.

```{r probeval2}
publicval <- c(0, 1)
yearval <- c(0, 1, 2)
```

Test the rate of changes at each type of school.

```{r simpleslope21}
quick_sim_slopes(model=m24, pred="yearc", modx="public", modx.values=publicval)
```

Plot the rate of changes at each type of school.

```{r plotslope21}
interact_plot(model=m24, pred="yearc", modx="public", modx.values=publicval)
```

Test the differences between each type of schools at each year.

```{r simpleslope22}
quick_sim_slopes(model=m24, pred="public", modx="yearc", modx.values=yearval)
```

Plot the differences between each type of schools at each year.

```{r plotslope22}
interact_plot(model=m24, pred="public", modx="yearc", modx.values=yearval)
```

## R-Squared

##### Example 1: Math Achievement Data across Schools and Countries 

Use the final model with cross-level interactions to find the proportion of variances at each level explained.

```{r model3}
m20 <- lmer(math ~ 1 + diffschooliqc + diffcountryiqc + avecountryiqcc
            + diffprivate + aveprivatec 
            + qualityc + opportunityc
            + diffschooliqc*diffcountryiqc + diffschooliqc*avecountryiqcc
            + diffschooliqc*diffprivate + diffschooliqc*aveprivatec
            + diffschooliqc*qualityc + diffschooliqc*opportunityc
            + diffprivate*avecountryiqcc + diffprivate*aveprivatec
            + diffprivate*qualityc + diffprivate*opportunityc
            + (1|schoolid) + (1 + diffschooliqc + diffprivate|countryid), 
            data=dat1, REML=FALSE, 
            control = lmerControl(optimizer ="Nelder_Mead"))
```

Use the `model.matrix` function to get the `data.frame` containing predictors used in the model. Note that if `I()` is used, the values are already transformed in the result of the `model.matrix` function.

```{r modelmatrix3}
dat1_1 <- model.matrix(m20)
```

Save the results of `summary` function for future use.

```{r summary3}
sumoutm20 <- summary(m20)
```

Next, predictors at each level must be listed and put in the r-squared calculation function. Level-1 predictors are `diffschooliqc` and all interactions involving `diffschooliqc`.

```{r level1predictor3}
lv1_covs <- c("diffschooliqc", "diffschooliqc:diffcountryiqc", 
              "diffschooliqc:avecountryiqcc", 
              "diffschooliqc:diffprivate", 
              "diffschooliqc:aveprivatec", 
              "diffschooliqc:qualityc", 
              "diffschooliqc:opportunityc")
```

Level-2 predictors are `diffschooliqc`, `diffprivate`, and all interactions between level-2 and level-3 predictors.

```{r level2predictor3}
lv2_covs <- c("diffschooliqc", "diffprivate",
              "avecountryiqcc:diffprivate",
              "diffprivate:aveprivatec",
              "diffprivate:qualityc",
              "diffprivate:opportunityc")
```

Level-3 predictors are listed as follows.

```{r level3predictor3}
lv3_covs <- c("avecountryiqcc", "aveprivatec", 
              "qualityc", "opportunityc")
```

Next, the covariance matrix between random effects are extracted.

```{r randomeffectcovariance3}
ranefvar20 <- as.matrix(Matrix::bdiag(VarCorr(m20)))
```

The `r2mlm3_manual` function in the `r2mlm` package is used. The arguments of this function is explained as follows:

- `data` is the data frame from the `model.matrix` function
- `l1_covs` is the vector of level-1 predictors
- `l2_covs` is the vector of level-2 predictors
- `l3_covs` is the vector of level-3 predictors
- `random_covs12` is the vector of level-1 predictors that are varied across level-2 units
- `random_covs13` is the vector of level-1 predictors that are varied across level-3 units
- `random_covs23` is the vector of level-2 predictors that are varied across level-3 units 
- `gamma_1` is the fixed effect estimates of the level-1 predictors
- `gamma_2` is the fixed effect estimates of the level-2 predictors
- `gamma_3` is the fixed effect estimates of the level-3 predictors
- `Tau12` is the covariance matrix of the random intercept at level 2 attached with the random slopes of level-1 predictors across level-2 units
- `Tau13` is the covariance matrix of the random intercept at level 3 attached with the random slopes of level-1 predictors across level-3 units
- `Tau23` is the covariance matrix of the random intercept at level 3 attached with the random slopes of level-1 predictors across level-3 units
- `sigma2` is the level-1 residual variances
- `clustermeancentered` is whether all variables are centered at level-2 and level-3 means

```{r r2mlm3}
library(r2mlm)
r2mlm3_manual(data=dat1_1,
              l1_covs=lv1_covs,
              l2_covs=lv2_covs,
              l3_covs=lv3_covs,
              random_covs12 = NULL,
              random_covs13 = "diffschooliqc",
              random_covs23 = "diffprivate",
              gamma_1 = coef(sumoutm20)[lv1_covs, "Estimate"],
              gamma_2 = coef(sumoutm20)[lv2_covs, "Estimate"],
              gamma_3 = coef(sumoutm20)[lv3_covs, "Estimate"],
              Tau12 = ranefvar20[1, 1, drop=FALSE],
              Tau13 = ranefvar20[2:3, 2:3, drop=FALSE],
              Tau23 = ranefvar20[c(2,4), c(2,4), drop=FALSE],
              sigma2=getME(m20, "sigma")^2,
              clustermeancentered = TRUE)
```

##### Example 2: Change in schools' standardized scores

Use the final model with cross-level interactions to find the proportion of variances at each level explained. Note that this is the example where a predictor is not group-mean centered, which is `yearc`.

```{r model4}
m24 <- lmer(test ~ 1 + yearc*public + (1|timeid) + (1 + yearc|schoolid), 
            data=dat2, REML=FALSE, 
            control = lmerControl(optimizer ="Nelder_Mead"))
```

Use the `model.matrix` function to get the `data.frame` containing predictors used in the model. 

```{r modelmatrix4}
dat2_1 <- model.matrix(m24)
```

Attach level-2 and level-3 ID in the resulting data frame.

```{r modelmatrix42}
varnames <- colnames(dat2_1)
dat2_1 <- data.frame(dat2_1, timeid = dat2$timeid,schoolid = dat2$schoolid)
colnames(dat2_1)[1:length(varnames)] <- varnames
```

Save the results of `summary` function for future use.

```{r summary4}
sumoutm24 <- summary(m24)
```

Next, predictors at each level are listed. Note that this model does not have level-1 predictors.

```{r predictor4}
lv2_covs <- c("yearc", "yearc:public")
lv3_covs <- c("public")
```

Next, the covariance matrix between random effects are extracted.

```{r randomeffectcovariance4}
ranefvar24 <- as.matrix(Matrix::bdiag(VarCorr(m24)))
```

The `r2mlm3_manual` function in the `r2mlm` package is used. Most arguments of this function have explained above. However, for the model without cluster mean centered, some arguments are needed.

- `Tau2_noncmc` (replacing `Tau12`) is the covariance matrix of the random intercept at level 2 attached with the random slopes of level-1 predictors across level-2 units
- `Tau3_noncmc` (replacing `Tau13` and `Tau23`) is the covariance matrix of the random intercept at level 3 attached with the random slopes of level-1 and level-2 predictors across level-3 units
- `l2clusterID_noncmc` is the level-2 ID variable name
- `l3clusterID_noncmc` is the level-3 ID variable name

```{r r2mlm4}
r2mlm3_manual(data=dat2_1,
              l1_covs=NULL,
              l2_covs=lv2_covs,
              l3_covs=lv3_covs,
              random_covs12 = NULL,
              random_covs13 = NULL,
              random_covs23 = "yearc",
              gamma_1 = NULL,
              gamma_2 = coef(sumoutm24)[lv2_covs, "Estimate"],
              gamma_3 = coef(sumoutm24)[lv3_covs, "Estimate"],
              Tau2_noncmc = ranefvar24[1, 1,drop=FALSE],
              Tau3_noncmc = ranefvar24[2:3, 2:3, drop=FALSE],
              sigma2=getME(m24, "sigma")^2,
              l2clusterID_noncmc = "timeid",
              l3clusterID_noncmc = "schoolid",
              clustermeancentered = FALSE)
```