We did a multilevel modeling using random intercepts and slopes model with unstructured covariances and an interaction term. 1) First, do model-fit analysis using -2LL (Old model=Random Intercepts and Slopes Model) with unstructured covariances, New model = Random Intercepts and Slopes Model with unstructured covariances and an interaction term) (see slides 37 & 44)
# Firstly, we need to prepare the data.
library(haven)
Cosmetic_Surgery <- read_sav("Cosmetic_Surgery.sav")
head(Cosmetic_Surgery)## # A tibble: 6 x 9
## ID Post_QoL Base_QoL Surgery Clinic Age BDI Reason Gender
## <dbl> <dbl> <dbl> <dbl+lbl> <dbl> <dbl> <dbl> <dbl+lbl> <dbl+l>
## 1 1 71.3 73 0 [Waitin… 1 31 12 0 [Change … 0 [Fem…
## 2 2 77 74 0 [Waitin… 1 32 16 0 [Change … 0 [Fem…
## 3 3 73 80 0 [Waitin… 1 33 13 0 [Change … 0 [Fem…
## 4 4 68.9 76 0 [Waitin… 1 59 11 0 [Change … 1 [Mal…
## 5 5 69 71 0 [Waitin… 1 61 11 0 [Change … 1 [Mal…
## 6 6 68.5 72 0 [Waitin… 1 32 10 1 [Physica… 0 [Fem…## Take a glance of the data
summary(Cosmetic_Surgery)## ID Post_QoL Base_QoL Surgery
## Min. : 1.00 Min. :40.00 Min. :43.00 Min. :0.0000
## 1st Qu.: 69.75 1st Qu.:52.30 1st Qu.:57.00 1st Qu.:0.0000
## Median :138.50 Median :58.00 Median :63.00 Median :0.0000
## Mean :138.50 Mean :59.61 Mean :63.56 Mean :0.4746
## 3rd Qu.:207.25 3rd Qu.:67.00 3rd Qu.:71.00 3rd Qu.:1.0000
## Max. :276.00 Max. :88.20 Max. :91.00 Max. :1.0000
## Clinic Age BDI Reason
## Min. : 1.000 Min. :18.00 Min. : 0.00 Min. :0.0000
## 1st Qu.: 3.000 1st Qu.:31.00 1st Qu.:10.00 1st Qu.:0.0000
## Median : 6.000 Median :38.00 Median :20.00 Median :1.0000
## Mean : 5.754 Mean :39.17 Mean :23.05 Mean :0.6449
## 3rd Qu.: 8.000 3rd Qu.:48.00 3rd Qu.:36.00 3rd Qu.:1.0000
## Max. :10.000 Max. :65.00 Max. :63.00 Max. :1.0000
## Gender
## Min. :0.0000
## 1st Qu.:0.0000
## Median :0.0000
## Mean :0.4312
## 3rd Qu.:1.0000
## Max. :1.0000# Load the nlme package for the data analysis
library(nlme)
# Run a Random Intercepts and Slopes Model as the Old model
OldModel <-lme(Post_QoL ~ Surgery + Base_QoL, data = Cosmetic_Surgery, random = ~Surgery|Clinic, method = "ML")
summary(OldModel)## Linear mixed-effects model fit by maximum likelihood
## Data: Cosmetic_Surgery
## 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: 10logLik(OldModel)*-2 # the -2LL for OldModel is 1798.624## 'log Lik.' 1798.624 (df=7)# Then, we add an interaction term into the old model, and name it: New Model
NewModel <- lme(Post_QoL ~ Surgery + Base_QoL + Reason + Surgery:Reason, data = Cosmetic_Surgery, random = ~Surgery|Clinic, method = "ML")
summary(NewModel)## Linear mixed-effects model fit by maximum likelihood
## Data: Cosmetic_Surgery
## 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: 10logLik(NewModel)*-2 # the -2LL for NewModel is 1789.045## 'log Lik.' 1789.045 (df=9)# Use anova function to compare two models
anova(OldModel,NewModel)## Model df AIC BIC logLik Test L.Ratio p-value
## OldModel 1 7 1812.624 1837.966 -899.3119
## NewModel 2 9 1807.045 1839.629 -894.5226 1 vs 2 9.578516 0.0083# By comparing the Old Model with New Model, the delta difference is 9.58, and
# The Change is significant because the P value is less than 0.01 (0.0083). Therefore we choose random intercepts and slopes model with unstructured covariance and interaction term.
# Prepare the function first to automatically calculate the ICC
ICC <- function(out) {
varests <- as.numeric(VarCorr(out)[1:nrow(VarCorr(out))])
return(paste("ICC=",varests[1]/ (varests[1]+varests[-1])))
}
# Use the ICC function to calculate the ICC in New Model
ICC(NewModel) ## [1] "ICC= 0.505950710632734" "ICC= 0.470246407862322"# The ICC for our New Model is 47.02%
# This tell us that about 47% of the variability in outcome difference occurs between level-2 units (clinic), with the remaining 53% occurring within level-2. This represents a high level of dependence. It is good to use a multilevel method. # Check the four parameter estimates of fixed effects
summary(NewModel)## Linear mixed-effects model fit by maximum likelihood
## Data: Cosmetic_Surgery
## 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# According to the outcome:
# γ_00=42.52, t = 10.97, p < 0.001; the average intercept across level-2 unites is 42.52 and significant.
# γ_10=-3.19, t = -1.47, p = 0.145; the average surgery slope across level-2 units is -3.19 and not significant.
# γ_20=0.31, t = 5.74, p < 0.001; the average Base_QoL slope across level-2 units is 0.31 and significant.
# γ_30=-3.52, t = -3.11, p < 0.05; the average Reason slope across level-2 units is -3.52 and significant.
# γ_40=4.22, t = 2.48, p < 0.05; the average Surgery and Reason interaction slope across level-2 units is -4.22 and significant.(QoL After Surgery)_ij= 42.52 - 3.19 * (Surgery)_ij + 0.31*(Base_QoL)_ij - 3.52 * (Reason)_ij+ 4.22 * (ReasonXSurgery)_ij