---
title: "Multilevel analysis"
author: "Sunny Avry"
date: "November 28th, 2020"
output:
html_document:
toc: false
theme: united
---
This R notebook can be used as a template to carry out basic multilevel analyses. It is usable in Rstudio. It is directly adapted from the chapter 19 of *Discovering statistics using R* (https://www.discoveringstatistics.com/books/discovering-statistics-using-r/).
## Packages
```{r message=FALSE, warning=FALSE}
library(ggplot2); library(nlme); library(reshape); library(ggpubr)
```
## Assumptions
Mutlilevel linear models are an extension of regression, so all the assumptions for regression apply to multilevel models. However, a lack of independence caused by higher-level variable can be solved using a multilevel model.
There are two additional assumptions in multilevel models. Random coefficients (intercepts and slopes) are assumed to be normally distributed around the overall model.
## Centring variables
### Grand mean centring
For a given variable we take each score and subtract from it the mean of all scores (for that variable).
### Group mean centring
For a given variable we take each score and subtract from it the mean of the scores (for that variable) within a given group.
### Advantages and consequences of centring variables
If group mean centring is used then a level 1 variable is typically centred around means of a level 2 variable
Sometimes, the intercept of outcome when all the predictors take a value of 0 is **meaningless** (e.g., heart rate at 0). Centring can also be a useful way to combat multicollinearity between predictors. Centring a predictor around its mean changes the meaning of the intercept, which becomes the value of the outcome when the predictor is its average value.
Fitting a multilevel model using the raw score predictors or using grand mean centred predictors are equivalent. The parameters $\beta_s$ will be different but you can find them again by transformation. Grand mean centring does not change the model but change the interpretation of the parameters (intercepts are not raw scores anymore).
Fitting a multilevel model using the raw score predictors gives a different model than using group mean centred predictors. The raw score model is not equivalent to the centred model.
Basically:
1) Group mean clustering should be used if the primary interest is in an association between variables measured at level 1
2) Grand mean centring is appropriate when the primary interest is in the level 2 variable but there is a neeed to control for a level 1 covariate
3) Both types of centring can be used to look at the differential influence of a variable at level 1 and 2
4) Group mean centring is preferable for examining cross-level interactions
## Importing data
```{r}
data <- read.table("Cosmetic Surgery.dat", sep = "\t", header = TRUE)
head(data)
```
## Exploring data
```{r}
dim(data) #number of rows and columns
```
```{r}
str(data) #basic description of the dataframe
```
## Variables
* **Post_QoL**: This is a measure of quality of life after the cosmetic surgery.
This is our outcome variable.
* **Base_QoL**: We need to adjust our outcome for quality of life before the
surgery.
* **Surgery**: This variable is a dummy variable that specifies whether the
person has undergone cosmetic surgery (1) or whether they are on the
waiting list (0), which acts as our control group.
* **Surgery_Text**: This variable is the same as above but specifies group
membership as text (we will use this variable when we create graphs but
not for the main analysis).
* **Clinic**: This variable specifies which of 10 clinics the person attended to
have their surgery.
* **Age**: This variable tells us the person’s age in years.
* **BDI**: It is becoming increasingly apparent that people volunteering for
cosmetic surgery (especially when the surgery is purely for vanity) might
have very different personality profiles than the general public (Cook,
Rossera, Toone, James, & Salmon, 2006). In particular, these people might
have low self-esteem or be depressed. When looking at quality of life it is
important to assess natural levels of depression, and this variable used the
Beck Depression Inventory (BDI) to do just that.
* **Reason**: This dummy variable specifies whether the person had/is waiting
to have surgery purely to change their appearance (0), or because of a
physical reason (1).
# Surgery and baseline quality of life (Base_QoL) as a predictors of quality of life after surgery (Post_QoL)
## Plotting the data
```{r}
scatter1 <-ggplot(data[which(data$Clinic==1),], aes(Base_QoL, Post_QoL , colour = Surgery_Text))
clinic1 <- scatter1 + geom_point() + rremove("x.title") + geom_smooth(method=lm, se = T, alpha = 0.1) + labs(x = "Base_QoL", y = "Post_QoL", colour = "Surgery")
scatter2 <-ggplot(data[which(data$Clinic==2),], aes(Base_QoL, Post_QoL , colour = Surgery_Text))
clinic2 <- scatter2 + geom_point() + rremove("x.title") + geom_smooth(method=lm, se = T, alpha = 0.1) + labs(x = "Base_QoL", y = "Post_QoL", colour = "Surgery")
scatter3 <-ggplot(data[which(data$Clinic==3),], aes(Base_QoL, Post_QoL , colour = Surgery_Text))
clinic3 <- scatter3 + geom_point() + rremove("x.title") + geom_smooth(method=lm, se = T, alpha = 0.1) + labs(x = "Base_QoL", y = "Post_QoL", colour = "Surgery")
scatter4 <-ggplot(data[which(data$Clinic==4),], aes(Base_QoL, Post_QoL , colour = Surgery_Text))
clinic4 <- scatter4 + geom_point() + rremove("x.title") + rremove("y.title") + geom_smooth(method=lm, se = T, alpha = 0.1) + labs(x = "Base_QoL", y = "Post_QoL", colour = "Surgery")
scatter5 <-ggplot(data[which(data$Clinic==5),], aes(Base_QoL, Post_QoL , colour = Surgery_Text))
clinic5 <- scatter5 + geom_point() + rremove("x.title") + rremove("y.title") + geom_smooth(method=lm, se = T, alpha = 0.1) + labs(x = "Base_QoL", y = "Post_QoL", colour = "Surgery")
scatter6 <-ggplot(data[which(data$Clinic==6),], aes(Base_QoL, Post_QoL , colour = Surgery_Text))
clinic6 <- scatter6 + geom_point() + rremove("x.title") + rremove("y.title") + geom_smooth(method=lm, se = T, alpha = 0.1) + labs(x = "Base_QoL", y = "Post_QoL", colour = "Surgery")
scatter7 <-ggplot(data[which(data$Clinic==7),], aes(Base_QoL, Post_QoL , colour = Surgery_Text))
clinic7 <- scatter7 + geom_point() + rremove("x.title") + rremove("y.title") + geom_smooth(method=lm, se = T, alpha = 0.1) + labs(x = "Base_QoL", y = "Post_QoL", colour = "Surgery")
scatter8 <-ggplot(data[which(data$Clinic==8),], aes(Base_QoL, Post_QoL , colour = Surgery_Text))
clinic8 <- scatter8 + geom_point() + rremove("x.title") + rremove("y.title") + geom_smooth(method=lm, se = T, alpha = 0.1) + labs(x = "Base_QoL", y = "Post_QoL", colour = "Surgery")
scatter9 <-ggplot(data[which(data$Clinic==9),], aes(Base_QoL, Post_QoL , colour = Surgery_Text))
clinic9 <- scatter9 + geom_point() + rremove("x.title") + rremove("y.title") + geom_smooth(method=lm, se = T, alpha = 0.1) + labs(x = "Base_QoL", y = "Post_QoL", colour = "Surgery")
scatter10 <-ggplot(data[which(data$Clinic==10),], aes(Base_QoL, Post_QoL , colour = Surgery_Text))
clinic10 <- scatter10 + geom_point() + rremove("x.title") + rremove("y.title") + geom_smooth(method=lm, se = T, alpha = 0.1) + labs(x = "Base_QoL", y = "Post_QoL", colour = "Surgery")
```
## Grid of plots
```{r message=FALSE, warning=FALSE}
figure <- ggarrange(clinic1, clinic2, clinic3, clinic4, clinic5, clinic6, clinic7, clinic8, clinic9, clinic10, common.legend = TRUE, ncol = 5, nrow = 2)
annotate_figure(figure,
top = text_grob("Quality of Life Pre-Post Surgery at 10 Clinics", face = "bold", size = 10),
bottom = text_grob("Pre_QoL", size = 10),
left = text_grob("Post_QoL", size = 10, rot = 90))
```
## ANOVA using Surgery as a predictor and Post_QoL as the outcome
```{r}
surgery_aov <- aov(Post_QoL ~ Surgery, data = data)
summary(surgery_aov)
surgery_linearModel<-lm(Post_QoL ~ Surgery, data = data)
summary(surgery_linearModel)
```
## ANCOVA using Surgery as a predictor, Post_QoL as the outcome and Pre_QoL as covariate
```{r}
surgery_linearModel<-lm(Post_QoL ~ Surgery + Base_QoL, data = data)
summary(surgery_linearModel)
```
At this point, we do not have taken into account the dependence that exists between participants of the same clinic. So we need to check evidence of variation across contexts, i.e., dependence across participants from the same clinic. This is done using the gls() function (generalized least squares).
## Multilevel modelling
### Null model
```{r}
interceptOnly <-gls(Post_QoL ~ 1, data = data, method = "ML")
summary(interceptOnly)
```
### Adding Random intercepts
```{r}
randomInterceptOnly <-lme(Post_QoL ~ 1, data = data, random
= ~1|Clinic, method = "ML")
summary(randomInterceptOnly)
```
### Computing Log-likelihood
```{r}
logLik(interceptOnly)*-2
logLik(randomInterceptOnly)*-2
```
### Testing if the improvement of our new model is significant.
If significant, we can say that intercepts vary significantly acrosss clinics
```{r}
anova(interceptOnly, randomInterceptOnly)
```
### Adding the fixed effect
```{r}
randomInterceptSurgery <-lme(Post_QoL ~ Surgery, data =
data, random = ~1|Clinic, method = "ML")
summary(randomInterceptSurgery)
```
Surgery here does not improve the model. The effect of clinic is significant but not the effect of surgery
### Adding the covariate
```{r}
randomInterceptSurgeryQoL <-lme(Post_QoL ~ Surgery + Base_QoL, data
= data, random = ~1|Clinic, method = "ML")
summary(randomInterceptSurgeryQoL)
```
### Testing if the improvement of our new model is significant
```{r}
anova(randomInterceptOnly,
randomInterceptSurgery, randomInterceptSurgeryQoL)
```
### Adding random slopes
```{r}
addRandomSlope<-lme(Post_QoL ~ Surgery + Base_QoL, data =
data, random = ~Surgery|Clinic, method = "ML")
summary(addRandomSlope)
```
### Testing the model with random intercepts and slopes
```{r}
anova(randomInterceptSurgeryQoL,addRandomSlope)
```
### Grouped scatterplot
```{r message=FALSE, warning=FALSE}
predscore <- fitted(addRandomSlope)
#clinic = as.factor(data$Clinic)
datapred <- data.frame(cbind(predscore=predscore,
Base_QoL=data$Base_QoL, Post_QoL=data$Post_QoL, clinic=data$Clinic))
datapred <- transform(datapred, clinic = as.factor(clinic))
clinics <- ggplot(datapred, aes(Base_QoL,predscore, colour = clinic)) + geom_point() + geom_smooth(method=lm, se = F, alpha = 0.1) + labs(title = "Predicted values from the model plotted againt the observed values (in grey)",x = "Base_QoL", y = "Post_QoL", colour = "Clinic")
clinics + geom_point(aes(y = Post_QoL), colour = "grey")
```
### Adding an interaction term
We add the reason for the person having cosmetic surgery as a fixed effect
```{r}
addReason<-lme(Post_QoL ~ Surgery + Base_QoL + Reason, data =
data, random = ~Surgery|Clinic, method = "ML")
summary(addReason)
```
We add the interaction between the reason and the surgery as a fixed effect
```{r}
finalModel<-lme(Post_QoL ~ Surgery + Base_QoL + Reason +
Surgery:Reason, data = data, random = ~Surgery|Clinic,
method = "ML")
summary(finalModel)
```
### Testing the model with the new fixed effects
```{r}
anova(addRandomSlope,addReason,finalModel)
```
### Computing confidence intervals for the estimates
```{r}
intervals(finalModel, 0.95)
```
### Final model interpretation
Physical reason = 1 and Cosmetic reason = 0
A surgery for medical reason is related to decrease of quality of life. However, this effect includes both people waiting for a surgery and who already had a surgery. So it is not very informative. However, the interaction between reason and surgery allow us to differentiate people who had a surgery from people who had not and their reason for having a surgery. We can also break down this interaction to better understand what is going on.
```{r}
cosmeticSubset<-data$Reason==0 #subset of pople having a surgery for a cosmetic reason
physicalSubset<-data$Reason==1 #subset of pople having a surgery for a physical reason
```
```{r}
physicalModel<-lme(Post_QoL ~ Surgery + Base_QoL, data =
data, random = ~Surgery|Clinic, subset= physicalSubset,
method = "ML")
summary(physicalModel)
```
```{r}
cosmeticModel<-lme(Post_QoL ~ Surgery + Base_QoL, data =
data, random = ~Surgery|Clinic, subset= cosmeticSubset,
method = "ML")
summary(cosmeticModel)
```
People who had a surgery for physical reason do not increase their quality of life (p-value is non significative). People who had a surgery for cosmetic reason decrease their quality of life (p-value is marginally significative).
### Reporting multilevel models
If you have built up your model from one with only fixed parameters to one with a random intercept, and then random slope, it is advisable to report all stages of this process (or at the very least report the fixed-effects-only model and the final model).
Example:
The relationship between surgery and quality of life showed significant variance in intercepts across participants, SD = 5.48 (95% CI: 3.31, 9.07), χ2(1) = 107.65, p < .0001. In addition, the slopes varied across participants, SD = 5.42 (3.13, 9.37), χ2 (2) = 38.87, p < .0001, and the slopes intercepts were negatively and significantly correlated, cor = –.95 (-.99, -.60).
Quality of life before surgery significantly predicted quality of life a surgery, b = 0.31, t(262) = 5.75, p < .001, surgery did not significantly predict quality of life, b = –3.19, t(262) = –1.46, p = .15, but the reason for surgery, b = –3.52, t(262) = –3.08, p < .01, and the interaction of the reason for surgery and surgery, b = 4.22, t(262) = 2.48, p < .05, both significantly predict quality of life. This interaction was broken down conducting separate multilevel models on the ‘physical reason’ and ‘attractiveness reason’. The models specified were the same as the main model but excluded the main effect and interaction term involving the reason for surgery. These analyses showed that for those operated on only to chatheir appearance, surgery almost significantly predicted quality of life a surgery, b = –4.31, t(87) = –1.89, p = .06: quality of life was lower a surgery compared to the control group. However, for those who had surgery to solve a physical problem, surgery did not significantly predict quality of life, b = 1.20, t(166) = 0.57, p = .57. The interaction effect, therefore, reflects the difference in slopes for surgery as a predictor of quality of life in those who had surgery for physical problems (slight positive slope) and those who had surgery purely for vanity (a negative slope).
| | *b* | *SE b* | *95% CI* |
|---|---|---|---|
| Baseline QoL | 0.31 | 0.05 | 0.20, 0.41 |
| Surgery | -3.19 | 2.19 | -7.45, 1.08 |
| Reason | -3.52 | 1.14 | -5.74, -1.29 |
| Surgery x Reason | 4.22 | 1.70 | 0.90,7.54 |
# Sources
Field, Andy, Jeremy Miles, and Zoë Field. "Discovering statistics using R." (2012).