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

## Multivariate Latent Covariate Model

##### Example 1: Group Flow on Counseling Satisfaction

Read the new data set about clients nested in counselors in group therapy.

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

See descriptive statistics of all variables within the dataset.

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

Next, run the multilevel latent covariate (MLC) model. In this model, each variable will be separated into two variables: between-group deviation (Level 2) and within-group deviation (Level 1). MLC can be run by any programs that analyze multilevel structural equation modeling (MSEM). The `lavaan` package is able to handle it. To run `lavaan`, users specify the relationship among variables via a text object. In the text, `~` represents regression where the variable on the left is outcome and the variables on the right are predictors. Multiple predictors can be combined by `+` like the formula in `lm` function. `~~` represents covariances, as shown at the end of this page. `level: 1` and `level: 2` tell `lavaan` that the following relationships belong to Level 1 or 2. In this example, the clients' satisfaction on group therapy is predicted by the perception of group flows in the group therapy.

When the script is created, the script is run via the `sem` function. Put the script in the `model` argument. Put the data object in the `data` argument. `cluster` means the variable name of the group ID. After running the result, users can use `summary` function to see the results.

```{r mlc1}
library(lavaan)
model1 <- '
level: 1
sat ~ groupflow
level: 2
sat ~ groupflow
'
fit1 <- sem(model = model1, data = dat, cluster = "groupid")
summary(fit1)
```

The similar model in multilevel analysis is Multivariate Manifest Covariate Model (MMC). The outcome is predicted by the group means of predictors and the group-mean centered predictors, as studied in Lecture 7. First, the `groupflow` is used to create group means and group-mean-centered variable.

```{r calcgroupmean1}
dat$aveflow <- ave(dat$groupflow, dat$groupid)
dat$diffflow <- dat$groupflow - dat$aveflow
```

Next, use the `lme4` package to run MMC.

```{r mmcmlm1}
library(lme4)
out2 <- lmer(sat ~ 1 + diffflow + aveflow + (1|groupid), data=dat, REML=FALSE)
summary(out2)
```

Users can use `lavaan` to run MMC too. Note that `lavaan` cannot run restricted maximum likelihood (`REML`). Full information likelihood is used in both `lavaan` and `lme4` packages.

```{r mmclavaan1}
model2 <- '
level: 1
sat ~ diffflow
level: 2
sat ~ aveflow
'
fit2 <- sem(model = model2, data = dat, cluster = "groupid")
summary(fit2)
```

The regression coefficients from MLC and MMC were different. In this case, group flow is the property of groups so MLC is more appropriate. 

##### Example 2: Emotional Stability on Group Counseling Satisfaction

Emotional stability is the property of individuals so it represents the formative measurement model. MMC is more appropriate, especially, when the sampling ratio is close to 100% like in this example. The differences between MLC and MMC are shown here. First, MLC is analyzed. 

```{r mlcemo1}
model3 <- '
level: 1
sat ~ emostability
level: 2
sat ~ emostability
'
fit3 <- sem(model = model3, data = dat, cluster = "groupid")
summary(fit3)
```

Next, MMC is analyzed by the `lme4` package. 

```{r mmcemo1}
dat$aveemo <- ave(dat$emostability, dat$groupid)
dat$diffemo <- dat$emostability - dat$aveemo
out4 <- lmer(sat ~ 1 + diffemo + aveemo + (1|groupid), data=dat, REML=FALSE)
summary(out4)
```

MMC can be analyzed by the `lavaan` package as well. 

```{r mmcemolavaan1}
model4 <- '
level: 1
sat ~ diffemo
level: 2
sat ~ aveemo
'
fit4 <- sem(model = model4, data = dat, cluster = "groupid")
summary(fit4)
```

Similary, the results from MLC and MMC were different but, for emotional stability, MMC results are more accurate.

## Multilevel Mediation

##### Example 3: Perceived Counseling Proficiency on Counseling Satisfaction

Because `lavaan` is the package for MSEM, users can run path analysis in each level. In this example, the effects of counselor's perceived counseling proficiency on counseling satisfaction is investigated. Perceived group flow could be the mediator of the effect between the perceived counseling proficiency and counseling satisfaction. 

Similary, `~` represents the regression effect. `*` in front of each predictor indicates that the names in front of `*` are the labels of each regression coefficient. The labels will be used at the end of script where labels are multiplied by each other. The multiplication is used to create the products of regression coefficients, i.e., indirect effects. The results of the `summary` function show the significance of each regression coefficients and the indirect effects. Unfortunately, nonparametric bootstrap cannot be run in multilevel models; therefore, the delta method, which assume symmetric sampling distribution, is used here. 

```{r multilevelmediation}
model5 <- '
level: 1
groupflow ~ aw*perceivedability
sat ~ bw*groupflow + perceivedability
level: 2
groupflow ~ ab*perceivedability + experience
sat ~ bb*groupflow + perceivedability + experience
abw := aw*bw
abb := ab*bb
'
fit5 <- sem(model = model5, data = dat, cluster = "groupid")
summary(fit5)
```

See the lecture for the interpretation of the results.

## Multilevel Correlation Matrix

##### Example 4: Within- and Between-Group Relationship among Variables in Counseling Study

`lavaan` is the MSEM package so it can analyze covariances among variables as well. `~~` represents relationships among variables. For example, `emostability ~~ groupflow + perceivedability` means that the covariance between emotional stability and group flows and the covariance between emotional stability and perceived proficiency are estimated. In this example, level-2 variable (`experience`) is included too. This variable will show in `level: 2` only.

```{r multilevelcovariance}
model6 <- '
level: 1
sat ~~ emostability + groupflow + perceivedability
emostability ~~ groupflow + perceivedability
groupflow ~~ perceivedability
level: 2
sat ~~ emostability + groupflow + perceivedability + experience
emostability ~~ groupflow + perceivedability + experience
groupflow ~~ perceivedability + experience
perceivedability ~~ experience
'
fit6 <- sem(model = model6, data = dat, cluster = "groupid")
summary(fit6, standardize=TRUE)
```

See the right column for the correlation (standardized covariances). The correlation can be tested whether they are significantly different from 0 by the `standardizedsolution` function.

```{r multilevelcorrelation}
standardizedsolution(fit6)
```