---
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).
