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

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

data <- read.table("Cosmetic Surgery.dat", sep = "\t", header = TRUE)
head(data)
##   particnu Post_QoL Base_QoL Clinic Surgery Reason Age Gender BDI Surgery_Text
## 1        1     71.3       73      1       0      0  31      0  12 Waiting List
## 2        2     77.0       74      1       0      0  32      0  16 Waiting List
## 3        3     73.0       80      1       0      0  33      0  13 Waiting List
## 4        4     68.9       76      1       0      0  59      1  11 Waiting List
## 5        5     69.0       71      1       0      0  61      1  11 Waiting List
## 6        6     68.5       72      1       0      1  32      0  10 Waiting List
##         Reason_Text Gender_Text
## 1 Change Appearance      Female
## 2 Change Appearance      Female
## 3 Change Appearance      Female
## 4 Change Appearance        Male
## 5 Change Appearance        Male
## 6   Physical reason      Female

## Exploring data

dim(data) #number of rows and columns
## [1] 276  12
 str(data) #basic description of the dataframe
## 'data.frame':    276 obs. of  12 variables:
##  $particnu : int 1 2 3 4 5 6 7 8 9 10 ... ##$ Post_QoL    : num  71.3 77 73 68.9 69 68.5 70 75 61.5 68 ...
##  $Base_QoL : int 73 74 80 76 71 72 71 73 80 64 ... ##$ Clinic      : int  1 1 1 1 1 1 1 1 1 1 ...
##  $Surgery : int 0 0 0 0 0 0 0 0 1 0 ... ##$ Reason      : int  0 0 0 0 0 1 1 1 0 0 ...
##  $Age : int 31 32 33 59 61 32 33 35 25 55 ... ##$ Gender      : int  0 0 0 1 1 0 0 0 0 1 ...
##  $BDI : int 12 16 13 11 11 10 11 15 30 36 ... ##$ Surgery_Text: Factor w/ 2 levels "Cosmetic Surgery",..: 2 2 2 2 2 2 2 2 1 2 ...
##  $Reason_Text : Factor w/ 2 levels "Change Appearance",..: 1 1 1 1 1 2 2 2 1 1 ... ##$ Gender_Text : Factor w/ 2 levels "Female","Male": 1 1 1 2 2 1 1 1 1 2 ...

## 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

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

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

surgery_aov <- aov(Post_QoL ~ Surgery, data = data)
summary(surgery_aov)
##              Df Sum Sq Mean Sq F value Pr(>F)
## Surgery       1     29   28.62    0.33  0.566
## Residuals   274  23748   86.67
surgery_linearModel<-lm(Post_QoL ~  Surgery,    data    =   data)
summary(surgery_linearModel)
##
## Call:
## lm(formula = Post_QoL ~ Surgery, data = data)
##
## Residuals:
##     Min      1Q  Median      3Q     Max
## -19.916  -7.271  -1.271   7.084  28.284
##
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  59.9159     0.7731  77.498   <2e-16 ***
## Surgery      -0.6449     1.1222  -0.575    0.566
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.31 on 274 degrees of freedom
## Multiple R-squared:  0.001204,   Adjusted R-squared:  -0.002442
## F-statistic: 0.3302 on 1 and 274 DF,  p-value: 0.566

## ANCOVA using Surgery as a predictor, Post_QoL as the outcome and Pre_QoL as covariate

surgery_linearModel<-lm(Post_QoL ~  Surgery +   Base_QoL,   data = data)
summary(surgery_linearModel)
##
## Call:
## lm(formula = Post_QoL ~ Surgery + Base_QoL, data = data)
##
## Residuals:
##      Min       1Q   Median       3Q      Max
## -13.4142  -5.1326  -0.6495   4.0540  23.5005
##
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 18.14702    2.90767   6.241 1.65e-09 ***
## Surgery     -1.69723    0.84404  -2.011   0.0453 *
## Base_QoL     0.66504    0.04537  14.659  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.977 on 273 degrees of freedom
## Multiple R-squared:  0.4411, Adjusted R-squared:  0.437
## F-statistic: 107.7 on 2 and 273 DF,  p-value: < 2.2e-16

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

interceptOnly   <-gls(Post_QoL ~    1,  data    =   data,   method  = "ML")
summary(interceptOnly)
## Generalized least squares fit by maximum likelihood
##   Model: Post_QoL ~ 1
##   Data: data
##        AIC      BIC    logLik
##   2017.124 2024.365 -1006.562
##
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 59.60978 0.5596972 106.5036       0
##
## Standardized residuals:
##        Min         Q1        Med         Q3        Max
## -2.1127754 -0.7875625 -0.1734394  0.7962286  3.0803354
##
## Residual standard error: 9.281527
## Degrees of freedom: 276 total; 275 residual

randomInterceptOnly <-lme(Post_QoL  ~   1,  data    =   data,   random
=   ~1|Clinic,  method  =   "ML")
summary(randomInterceptOnly)
## Linear mixed-effects model fit by maximum likelihood
##  Data: data
##        AIC      BIC    logLik
##   1911.473 1922.334 -952.7364
##
## Random effects:
##  Formula: ~1 | Clinic
##         (Intercept) Residual
## StdDev:    5.909691 7.238677
##
## Fixed effects: Post_QoL ~ 1
##                Value Std.Error  DF  t-value p-value
## (Intercept) 60.08377  1.923283 266 31.24022       0
##
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max
## -1.8828507 -0.7606631 -0.1378732  0.7075242  2.8607949
##
## Number of Observations: 276
## Number of Groups: 10

### Computing Log-likelihood

logLik(interceptOnly)*-2
## 'log Lik.' 2013.124 (df=2)
logLik(randomInterceptOnly)*-2
## 'log Lik.' 1905.473 (df=3)

### Testing if the improvement of our new model is significant.

If significant, we can say that intercepts vary significantly acrosss clinics

anova(interceptOnly,    randomInterceptOnly)
##                     Model df      AIC      BIC     logLik   Test  L.Ratio
## interceptOnly           1  2 2017.124 2024.365 -1006.5622
## randomInterceptOnly     2  3 1911.473 1922.334  -952.7364 1 vs 2 107.6517
##                     p-value
## interceptOnly
## randomInterceptOnly  <.0001

randomInterceptSurgery  <-lme(Post_QoL ~    Surgery,    data    =
data,   random  =   ~1|Clinic,  method  =   "ML")
summary(randomInterceptSurgery)
## Linear mixed-effects model fit by maximum likelihood
##  Data: data
##        AIC      BIC    logLik
##   1910.137 1924.619 -951.0686
##
## Random effects:
##  Formula: ~1 | Clinic
##         (Intercept) Residual
## StdDev:    6.099513  7.18542
##
## Fixed effects: Post_QoL ~ Surgery
##                Value Std.Error  DF  t-value p-value
## (Intercept) 59.30517 2.0299632 265 29.21490   0.000
## Surgery      1.66583 0.9091314 265  1.83233   0.068
##  Correlation:
##         (Intr)
## Surgery -0.21
##
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max
## -1.8904290 -0.7191399 -0.1420998  0.7177762  2.8644538
##
## Number of Observations: 276
## Number of Groups: 10

Surgery here does not improve the model. The effect of clinic is significant but not the effect of surgery

randomInterceptSurgeryQoL   <-lme(Post_QoL  ~   Surgery +   Base_QoL,   data
=   data,   random  =   ~1|Clinic,  method  =   "ML")
summary(randomInterceptSurgeryQoL)
## Linear mixed-effects model fit by maximum likelihood
##  Data: data
##       AIC      BIC   logLik
##   1847.49 1865.592 -918.745
##
## Random effects:
##  Formula: ~1 | Clinic
##         (Intercept) Residual
## StdDev:    3.039264 6.518986
##
## Fixed effects: Post_QoL ~ Surgery + Base_QoL
##                 Value Std.Error  DF   t-value p-value
## (Intercept) 29.563601  3.471879 264  8.515160  0.0000
## Surgery     -0.312999  0.843145 264 -0.371228  0.7108
## Base_QoL     0.478630  0.052774 264  9.069465  0.0000
##  Correlation:
##          (Intr) Surgry
## Surgery   0.102
## Base_QoL -0.947 -0.222
##
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max
## -1.8872666 -0.7537675 -0.0954987  0.5657241  3.0020852
##
## Number of Observations: 276
## Number of Groups: 10

### Testing if the improvement of our new model is significant

anova(randomInterceptOnly,
randomInterceptSurgery, randomInterceptSurgeryQoL)
##                           Model df      AIC      BIC    logLik   Test  L.Ratio
## randomInterceptOnly           1  3 1911.473 1922.334 -952.7364
## randomInterceptSurgery        2  4 1910.137 1924.619 -951.0686 1 vs 2  3.33564
## randomInterceptSurgeryQoL     3  5 1847.490 1865.592 -918.7450 2 vs 3 64.64721
##                           p-value
## randomInterceptOnly
## randomInterceptSurgery     0.0678
## randomInterceptSurgeryQoL  <.0001

addRandomSlope<-lme(Post_QoL    ~   Surgery +   Base_QoL,   data    =
data,   random  =   ~Surgery|Clinic,    method  =   "ML")
summary(addRandomSlope)
## Linear mixed-effects model fit by maximum likelihood
##  Data: data
##        AIC      BIC    logLik
##   1812.624 1837.967 -899.3119
##
## Random effects:
##  Formula: ~Surgery | Clinic
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev   Corr
## (Intercept) 6.132655 (Intr)
## Surgery     6.197489 -0.965
## Residual    5.912335
##
## Fixed effects: Post_QoL ~ Surgery + Base_QoL
##                Value Std.Error  DF   t-value p-value
## (Intercept) 40.10253  3.892945 264 10.301334  0.0000
## Surgery     -0.65453  2.110917 264 -0.310069  0.7568
## Base_QoL     0.31022  0.053506 264  5.797812  0.0000
##  Correlation:
##          (Intr) Surgry
## Surgery  -0.430
## Base_QoL -0.855 -0.063
##
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max
## -2.4114778 -0.6628574 -0.1138411  0.6833110  2.8334730
##
## Number of Observations: 276
## Number of Groups: 10

### Testing the model with random intercepts and slopes

anova(randomInterceptSurgeryQoL,addRandomSlope)
##                           Model df      AIC      BIC    logLik   Test  L.Ratio
## randomInterceptSurgeryQoL     1  5 1847.490 1865.592 -918.7450
## addRandomSlope                2  7 1812.624 1837.966 -899.3119 1 vs 2 38.86626
##                           p-value
## randomInterceptSurgeryQoL
## addRandomSlope             <.0001

### Grouped scatterplot

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

We add the reason for the person having cosmetic surgery as a fixed effect

addReason<-lme(Post_QoL ~   Surgery +   Base_QoL    +   Reason, data    =
data,   random  =   ~Surgery|Clinic,    method  =   "ML")
summary(addReason)
## Linear mixed-effects model fit by maximum likelihood
##  Data: data
##        AIC      BIC    logLik
##   1810.825 1839.788 -897.4124
##
## Random effects:
##  Formula: ~Surgery | Clinic
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev   Corr
## (Intercept) 5.838006 (Intr)
## Surgery     6.171602 -0.969
## Residual    5.886964
##
## Fixed effects: Post_QoL ~ Surgery + Base_QoL + Reason
##                Value Std.Error  DF   t-value p-value
## (Intercept) 41.43182  3.929923 263 10.542653  0.0000
## Surgery     -0.56817  2.106684 263 -0.269697  0.7876
## Base_QoL     0.30535  0.053417 263  5.716335  0.0000
## Reason      -1.69137  0.854411 263 -1.979571  0.0488
##  Correlation:
##          (Intr) Surgry Bas_QL
## Surgery  -0.400
## Base_QoL -0.862 -0.064
## Reason   -0.231 -0.042  0.118
##
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max
## -2.2006350 -0.6830482 -0.1256909  0.6730347  2.9817687
##
## Number of Observations: 276
## Number of Groups: 10

We add the interaction between the reason and the surgery as a fixed effect

finalModel<-lme(Post_QoL    ~   Surgery +   Base_QoL    +   Reason  +
Surgery:Reason, data    =   data,   random  =   ~Surgery|Clinic,
method  =   "ML")
summary(finalModel)
## Linear mixed-effects model fit by maximum likelihood
##  Data: data
##        AIC      BIC    logLik
##   1807.045 1839.629 -894.5226
##
## Random effects:
##  Formula: ~Surgery | Clinic
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev   Corr
## (Intercept) 5.482359 (Intr)
## Surgery     5.417495 -0.946
## Residual    5.818911
##
## Fixed effects: Post_QoL ~ Surgery + Base_QoL + Reason + Surgery:Reason
##                   Value Std.Error  DF   t-value p-value
## (Intercept)    42.51781  3.875317 262 10.971441  0.0000
## Surgery        -3.18768  2.185367 262 -1.458646  0.1459
## Base_QoL        0.30536  0.053125 262  5.747835  0.0000
## Reason         -3.51515  1.140934 262 -3.080939  0.0023
## Surgery:Reason  4.22129  1.700269 262  2.482718  0.0137
##  Correlation:
##                (Intr) Surgry Bas_QL Reason
## Surgery        -0.356
## Base_QoL       -0.865 -0.078
## Reason         -0.233  0.306  0.065
## Surgery:Reason  0.096 -0.505  0.024 -0.661
##
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max
## -2.2331483 -0.6972193 -0.1541073  0.6326387  3.1641799
##
## Number of Observations: 276
## Number of Groups: 10

### Testing the model with the new fixed effects

anova(addRandomSlope,addReason,finalModel)
##                Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## addRandomSlope     1  7 1812.624 1837.966 -899.3119
## addReason          2  8 1810.825 1839.788 -897.4124 1 vs 2 3.798961  0.0513
## finalModel         3  9 1807.045 1839.629 -894.5226 2 vs 3 5.779555  0.0162

### Computing confidence intervals for the estimates

intervals(finalModel,   0.95)
## Approximate 95% confidence intervals
##
##  Fixed effects:
##                     lower       est.      upper
## (Intercept)    34.9565174 42.5178133 50.0791092
## Surgery        -7.4516406 -3.1876776  1.0762854
## Base_QoL        0.2017009  0.3053562  0.4090115
## Reason         -5.7412741 -3.5151488 -1.2890236
## Surgery:Reason  0.9038215  4.2212894  7.5387574
## attr(,"label")
## [1] "Fixed effects:"
##
##  Random Effects:
##   Level: Clinic
##                               lower       est.     upper
## sd((Intercept))           3.3130516  5.4823585  9.072076
## sd(Surgery)               3.1323585  5.4174947  9.369697
## cor((Intercept),Surgery) -0.9937882 -0.9455544 -0.598258
##
##  Within-group standard error:
##    lower     est.    upper
## 5.331220 5.818911 6.351214

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

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
physicalModel<-lme(Post_QoL ~   Surgery +   Base_QoL,   data    =
data,   random  =   ~Surgery|Clinic,    subset= physicalSubset,
method  =   "ML")
summary(physicalModel)
## Linear mixed-effects model fit by maximum likelihood
##  Data: data
##   Subset: physicalSubset
##       AIC      BIC    logLik
##   1172.56 1194.832 -579.2798
##
## Random effects:
##  Formula: ~Surgery | Clinic
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev   Corr
## (Intercept) 5.773827 (Intr)
## Surgery     5.804865 -0.948
## Residual    5.798764
##
## Fixed effects: Post_QoL ~ Surgery + Base_QoL
##                Value Std.Error  DF  t-value p-value
## (Intercept) 38.02079  4.705980 166 8.079250  0.0000
## Surgery      1.19655  2.099769 166 0.569848  0.5696
## Base_QoL     0.31771  0.069471 166 4.573271  0.0000
##  Correlation:
##          (Intr) Surgry
## Surgery  -0.306
## Base_QoL -0.908 -0.078
##
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max
## -2.2447342 -0.6505340 -0.1264188  0.6111506  2.9472101
##
## Number of Observations: 178
## Number of Groups: 10
cosmeticModel<-lme(Post_QoL ~   Surgery +   Base_QoL,   data    =
data,   random  =   ~Surgery|Clinic,    subset= cosmeticSubset,
method  =   "ML")
summary(cosmeticModel)
## Linear mixed-effects model fit by maximum likelihood
##  Data: data
##   Subset: cosmeticSubset
##        AIC      BIC    logLik
##   650.9469 669.0417 -318.4734
##
## Random effects:
##  Formula: ~Surgery | Clinic
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev   Corr
## (Intercept) 5.006026 (Intr)
## Surgery     5.292027 -0.969
## Residual    5.738551
##
## Fixed effects: Post_QoL ~ Surgery + Base_QoL
##                Value Std.Error DF   t-value p-value
## (Intercept) 41.78605  5.573849 87  7.496802  0.0000
## Surgery     -4.30702  2.275002 87 -1.893193  0.0617
## Base_QoL     0.33849  0.080274 87  4.216720  0.0001
##  Correlation:
##          (Intr) Surgry
## Surgery  -0.252
## Base_QoL -0.937 -0.058
##
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max
## -1.8945645 -0.6616222 -0.1461451  0.6460834  2.6741347
##
## Number of Observations: 98
## Number of Groups: 9

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