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/).
library(ggplot2); library(nlme); library(reshape); library(ggpubr)
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.
For a given variable we take each score and subtract from it the mean of all scores (for that variable).
For a given variable we take each score and subtract from it the mean of the scores (for that variable) within a given group.
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:
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
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 ...
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")
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))
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
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).
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
logLik(interceptOnly)*-2
## 'log Lik.' 2013.124 (df=2)
logLik(randomInterceptOnly)*-2
## 'log Lik.' 1905.473 (df=3)
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
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
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
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
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
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
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).
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 |
Field, Andy, Jeremy Miles, and Zoë Field. “Discovering statistics using R.” (2012).