# Read Dataset
therapy <- read.csv(file = "D:/UNY/MySta/SEM 5/AMS/dataset_ams/therapy.csv")
str(therapy)
## 'data.frame': 112 obs. of 6 variables:
## $ obs: int 1 2 3 4 5 6 7 8 9 10 ...
## $ X : int 1 1 1 1 1 1 1 0 0 0 ...
## $ M : int 23 32 32 35 42 38 38 54 32 54 ...
## $ Y : int 31 41 41 50 41 44 44 61 42 61 ...
## $ i : int 1 2 3 4 5 6 7 1 2 3 ...
## $ j : int 1 1 1 1 1 1 1 2 2 2 ...
summary(therapy)
## obs X M Y i
## Min. : 1.00 Min. :0.0 Min. : 4.00 Min. : 2.00 Min. :1
## 1st Qu.: 28.75 1st Qu.:0.0 1st Qu.:27.00 1st Qu.:34.75 1st Qu.:2
## Median : 56.50 Median :0.5 Median :35.00 Median :43.00 Median :4
## Mean : 56.50 Mean :0.5 Mean :35.30 Mean :41.86 Mean :4
## 3rd Qu.: 84.25 3rd Qu.:1.0 3rd Qu.:42.25 3rd Qu.:47.75 3rd Qu.:6
## Max. :112.00 Max. :1.0 Max. :66.00 Max. :82.00 Max. :7
## j
## Min. : 1.00
## 1st Qu.: 4.75
## Median : 8.50
## Mean : 8.50
## 3rd Qu.:12.25
## Max. :16.00
# Grand Mead Clustering
therapy$M.gc <- therapy$M - mean(therapy$M)
# Mengecek nilai rata-rata dan varians
mean(therapy$M); var(therapy$M)
## [1] 35.30357
## [1] 136.88
mean(therapy$Y); var(therapy$Y)
## [1] 41.85714
## [1] 210.9344
# Model 0: Random intercept (null model)
library(nlme)
## Warning: package 'nlme' was built under R version 4.4.3
M0 <- lme(fixed = Y ~ 1, random = ~1|j, data = therapy)
summary(M0)
## Linear mixed-effects model fit by REML
## Data: therapy
## AIC BIC logLik
## 887.4176 895.5462 -440.7088
##
## Random effects:
## Formula: ~1 | j
## (Intercept) Residual
## StdDev: 9.645123 11.08758
##
## Fixed effects: Y ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 41.85714 2.62905 96 15.92101 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.382402691 -0.506970362 -0.001438152 0.300396467 3.414594811
##
## Number of Observations: 112
## Number of Groups: 16
library(lmeInfo)
## Warning: package 'lmeInfo' was built under R version 4.4.3
vc0 <- sqrt(diag(varcomp_vcov(M0)))
vc0
## Tau.j.var((Intercept)) sigma_sq
## 40.46141 17.74407
# Model 1: 𝑌diprediksi oleh X
library(nlme)
M1 <- lme(fixed = Y ~ X, random = ~1|j, data = therapy)
summary(M1)
## Linear mixed-effects model fit by REML
## Data: therapy
## AIC BIC logLik
## 874.396 885.1979 -433.198
##
## Random effects:
## Formula: ~1 | j
## (Intercept) Residual
## StdDev: 6.414986 11.08758
##
## Fixed effects: Y ~ X
## Value Std.Error DF t-value p-value
## (Intercept) 34.62500 2.709108 96 12.780962 0.000
## X 14.46429 3.831257 14 3.775337 0.002
## Correlation:
## (Intr)
## X -0.707
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.208994103 -0.597222780 0.007022532 0.384760090 3.644046482
##
## Number of Observations: 112
## Number of Groups: 16
library(lmeInfo)
vc1 <- sqrt(diag(varcomp_vcov(M1)))
vc1
## Tau.j.var((Intercept)) sigma_sq
## 22.33615 17.74406
# Model 2: 𝑌diprediksi oleh X dan M.gc
library(nlme)
M2 <- lme(fixed = Y ~ X + M.gc, random = ~1|j, data = therapy)
summary(M2)
## Linear mixed-effects model fit by REML
## Data: therapy
## AIC BIC logLik
## 811.639 825.0958 -400.8195
##
## Random effects:
## Formula: ~1 | j
## (Intercept) Residual
## StdDev: 6.357745 7.92601
##
## Fixed effects: Y ~ X + M.gc
## Value Std.Error DF t-value p-value
## (Intercept) 39.42724 2.532964 95 15.565654 0.0000
## X 4.85980 3.648941 14 1.331838 0.2042
## M.gc 0.82240 0.084155 95 9.772501 0.0000
## Correlation:
## (Intr) X
## X -0.720
## M.gc 0.194 -0.269
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.94069659 -0.45669625 -0.05496793 0.61276725 2.58566146
##
## Number of Observations: 112
## Number of Groups: 16
library(lmeInfo)
vc2 <- sqrt(diag(varcomp_vcov(M2)))
vc2
## Tau.j.var((Intercept)) sigma_sq
## 18.798501 9.114945
# Model 3: M diprediksi oleh X
library(nlme)
M3 <- lme(fixed = M ~ X, random = ~1|j, data = therapy)
summary(M1)
## Linear mixed-effects model fit by REML
## Data: therapy
## AIC BIC logLik
## 874.396 885.1979 -433.198
##
## Random effects:
## Formula: ~1 | j
## (Intercept) Residual
## StdDev: 6.414986 11.08758
##
## Fixed effects: Y ~ X
## Value Std.Error DF t-value p-value
## (Intercept) 34.62500 2.709108 96 12.780962 0.000
## X 14.46429 3.831257 14 3.775337 0.002
## Correlation:
## (Intr)
## X -0.707
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.208994103 -0.597222780 0.007022532 0.384760090 3.644046482
##
## Number of Observations: 112
## Number of Groups: 16
library(lmeInfo)
vc3 <- sqrt(diag(varcomp_vcov(M3)))
vc3
## Tau.j.var((Intercept)) sigma_sq
## 11.93692 12.50109
Group mean centering
mean.Mg <- with(therapy, tapply(M, j, FUN = mean))
mean.Mg
## 1 2 3 4 5 6 7 8
## 34.28571 35.28571 23.42857 41.00000 36.00000 33.57143 42.42857 26.28571
## 9 10 11 12 13 14 15 16
## 23.14286 48.42857 25.28571 37.57143 41.00000 29.85714 48.42857 38.85714
gm <- data.frame(j = 1:16, mean.Mg)
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.4.3
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:nlme':
##
## collapse
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
therapy2 <- full_join(therapy, gm, by = "j")
# group-mean centering
therapy2$M.gm <- therapy2$M - therapy2$mean.Mg
Model dengan𝑌diprediksi oleh 𝑀yang dipusatkan pada rata-rata grup
(group mean centered) dan X
library(nlme)
M4 <- lme(fixed = Y ~ X + M.gm, random = ~1 | j, data = therapy2)
summary(M4)
## Linear mixed-effects model fit by REML
## Data: therapy2
## AIC BIC logLik
## 814.0263 827.4831 -402.0132
##
## Random effects:
## Formula: ~1 | j
## (Intercept) Residual
## StdDev: 7.052398 7.92737
##
## Fixed effects: Y ~ X + M.gm
## Value Std.Error DF t-value p-value
## (Intercept) 34.62500 2.709103 95 12.780982 0.000
## X 14.46429 3.831250 14 3.775343 0.002
## M.gm 0.83748 0.086938 95 9.633080 0.000
## Correlation:
## (Intr) X
## X -0.707
## M.gm 0.000 0.000
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.84425171 -0.49850512 -0.04053074 0.61707575 2.56901768
##
## Number of Observations: 112
## Number of Groups: 16
vc4 <- sqrt(diag(varcomp_vcov(M4)))
vc4
## Tau.j.var((Intercept)) sigma_sq
## 22.229973 9.118248
Model dengan𝑌diprediksi oleh 𝑀yang dipusatkan pada rata-rata grup
(group mean centered) dan X dengan slope acak dalam grup
library(nlme)
M5 <- lme(fixed = Y ~ X + M.gm, random = ~M.gm | j, data = therapy2, control = lmeControl(opt = 'optim'))
summary(M5)
## Linear mixed-effects model fit by REML
## Data: therapy2
## AIC BIC logLik
## 812.7799 831.6193 -399.3899
##
## Random effects:
## Formula: ~M.gm | j
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 7.0241631 (Intr)
## M.gm 0.2732578 0.903
## Residual 7.6039708
##
## Fixed effects: Y ~ X + M.gm
## Value Std.Error DF t-value p-value
## (Intercept) 34.50009 2.497661 95 13.812959 0e+00
## X 14.71411 3.248616 14 4.529347 5e-04
## M.gm 0.73972 0.110076 95 6.720130 0e+00
## Correlation:
## (Intr) X
## X -0.650
## M.gm 0.427 -0.050
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.940076433 -0.406048901 -0.006356577 0.478167129 2.721989124
##
## Number of Observations: 112
## Number of Groups: 16
vc5 <- sqrt(diag(varcomp_vcov(M5)))
vc5
## Tau.j.var((Intercept)) Tau.j.cov(M.gm,(Intercept))
## 21.75681968 0.99774107
## Tau.j.var(M.gm) sigma_sq
## 0.06901692 8.96957176
Model dengan𝑌diprediksi oleh 𝑀yang dipusatkan pada rata-rata grup
serta rata-rata M di setiap grup (group mean centered)
library(nlme)
M6 <- lme(fixed = Y ~ X + M.gm + mean.Mg, random = ~1 | j, data = therapy2)
summary(M6)
## Linear mixed-effects model fit by REML
## Data: therapy2
## AIC BIC logLik
## 813.428 829.5207 -400.714
##
## Random effects:
## Formula: ~1 | j
## (Intercept) Residual
## StdDev: 6.501522 7.92737
##
## Fixed effects: Y ~ X + M.gm + mean.Mg
## Value Std.Error DF t-value p-value
## (Intercept) 17.036727 10.400653 95 1.638044 0.1047
## X 7.492934 5.366560 13 1.396227 0.1860
## M.gm 0.837480 0.086938 95 9.633080 0.0000
## mean.Mg 0.596935 0.342380 13 1.743485 0.1048
## Correlation:
## (Intr) X M.gm
## X 0.608
## M.gm 0.000 0.000
## mean.Mg -0.970 -0.745 0.000
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.90553748 -0.47310930 -0.01052624 0.62723230 2.57988475
##
## Number of Observations: 112
## Number of Groups: 16
vc6 <- sqrt(diag(varcomp_vcov(M6)))
vc6
## Tau.j.var((Intercept)) sigma_sq
## 20.143044 9.118248
Model multilevel dengan 𝑴.𝒈𝒎sebagai variabel dependen
library(nlme)
M7 <- lme(fixed = M.gm ~ X, random = ~1 | j, data = therapy2)
summary(M7)
## Linear mixed-effects model fit by REML
## Data: therapy2
## AIC BIC logLik
## 803.9985 814.8004 -397.9992
##
## Random effects:
## Formula: ~1 | j
## (Intercept) Residual
## StdDev: 0.0002492873 8.694079
##
## Fixed effects: M.gm ~ X
## Value Std.Error DF t-value p-value
## (Intercept) 1.209044e-15 1.161795 96 1.040669e-15 1
## X -2.378748e-15 1.643026 14 -1.447785e-15 1
## Correlation:
## (Intr)
## X -0.707
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.39900539 -0.71888004 0.07394195 0.57921192 2.69477318
##
## Number of Observations: 112
## Number of Groups: 16
vc7 <- sqrt(diag(varcomp_vcov(M7)))
vc7
## Tau.j.var((Intercept)) sigma_sq
## 4.368787 10.910046
Model multilevel dengan 𝑿 berkaitan dengan mean.Mg
library(nlme)
M8 <- lme(fixed = mean.Mg ~ X, random = ~1 | j, data = therapy2)
summary(M8)
## Linear mixed-effects model fit by REML
## Data: therapy2
## AIC BIC logLik
## -5774.142 -5763.34 2891.071
##
## Random effects:
## Formula: ~1 | j
## (Intercept) Residual
## StdDev: 5.588083 1.062107e-14
##
## Fixed effects: mean.Mg ~ X
## Value Std.Error DF t-value p-value
## (Intercept) 29.46429 1.975686 96 14.913448 0e+00
## X 11.67857 2.794042 14 4.179813 9e-04
## Correlation:
## (Intr)
## X -0.707
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -0.6689936 0.6689936 1.3379872 2.0069808 2.6759744
##
## Number of Observations: 112
## Number of Groups: 16