Multilevel analysis
Sunny Avry
November 28th, 2020
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:
- Group mean clustering should be used if the primary interest is in an association between variables measured at level 1
- 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
- Both types of centring can be used to look at the differential influence of a variable at level 1 and 2
- 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 FemaleExploring 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.67surgery_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.566ANCOVA 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-16At 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 residualAdding Random intercepts
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: 10Computing 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 <.0001Adding the fixed effect
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: 10Surgery here does not improve the model. The effect of clinic is significant but not the effect of surgery
Adding the covariate
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: 10Testing 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 <.0001Adding random slopes
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: 10Testing 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 <.0001Grouped 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")
Adding an interaction term
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: 10We 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: 10Testing 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.0162Computing 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.351214Final 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 reasonphysicalModel<-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: 10cosmeticModel<-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: 9People 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).