Read the new data set about clients nested in counselors in group therapy.
dat <- read.table("lecture10ex1.csv", sep=",", header=TRUE)
See descriptive statistics of all variables within the dataset.
library(psych)
describe(dat)
## vars n mean sd median trimmed mad min max
## clientid 1 1000 500.50 288.82 500.50 500.50 370.65 1.00 1000.00
## groupid 2 1000 100.50 57.76 100.50 100.50 74.13 1.00 200.00
## sat 3 1000 0.05 1.03 0.03 0.04 1.01 -2.89 4.26
## emostability 4 1000 -0.02 1.01 -0.03 -0.02 1.05 -3.65 3.04
## groupflow 5 1000 -0.08 0.99 -0.11 -0.09 0.99 -2.79 3.15
## perceivedability 6 1000 -0.08 0.98 -0.12 -0.08 1.00 -3.03 2.82
## experience 7 1000 -0.05 1.05 -0.01 -0.08 1.06 -2.76 3.33
## range skew kurtosis se
## clientid 999.00 0.00 -1.20 9.13
## groupid 199.00 0.00 -1.20 1.83
## sat 7.15 0.17 0.16 0.03
## emostability 6.69 -0.05 -0.01 0.03
## groupflow 5.94 0.16 0.02 0.03
## perceivedability 5.85 0.06 -0.13 0.03
## experience 6.09 0.28 0.15 0.03
Next, run the multilevel latent covariate (MLC) model. In this model,
each variable will be separated into two variables: between-group
deviation (Level 2) and within-group deviation (Level 1). MLC can be run
by any programs that analyze multilevel structural equation modeling
(MSEM). The lavaan
package is able to handle it. To run
lavaan
, users specify the relationship among variables via
a text object. In the text, ~
represents regression where
the variable on the left is outcome and the variables on the right are
predictors. Multiple predictors can be combined by +
like
the formula in lm
function. ~~
represents
covariances, as shown at the end of this page. level: 1
and
level: 2
tell lavaan
that the following
relationships belong to Level 1 or 2. In this example, the clients’
satisfaction on group therapy is predicted by the perception of group
flows in the group therapy.
When the script is created, the script is run via the
sem
function. Put the script in the model
argument. Put the data object in the data
argument.
cluster
means the variable name of the group ID. After
running the result, users can use summary
function to see
the results.
library(lavaan)
model1 <- '
level: 1
sat ~ groupflow
level: 2
sat ~ groupflow
'
fit1 <- sem(model = model1, data = dat, cluster = "groupid")
summary(fit1)
## lavaan 0.6.17 ended normally after 14 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 5
##
## Number of observations 1000
## Number of clusters [groupid] 200
##
## Model Test User Model:
##
## Test statistic 0.073
## Degrees of freedom 0
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
##
## Level 1 [within]:
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## sat ~
## groupflow 0.288 0.040 7.203 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .sat 0.817 0.041 20.000 0.000
##
##
## Level 2 [groupid]:
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## sat ~
## groupflow 0.528 0.072 7.367 0.000
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|)
## .sat 0.094 0.036 2.606 0.009
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .sat 0.084 0.027 3.137 0.002
The similar model in multilevel analysis is Multivariate Manifest
Covariate Model (MMC). The outcome is predicted by the group means of
predictors and the group-mean centered predictors, as studied in Lecture
7. First, the groupflow
is used to create group means and
group-mean-centered variable.
dat$aveflow <- ave(dat$groupflow, dat$groupid)
dat$diffflow <- dat$groupflow - dat$aveflow
Next, use the lme4
package to run MMC.
library(lme4)
out2 <- lmer(sat ~ 1 + diffflow + aveflow + (1|groupid), data=dat, REML=FALSE)
summary(out2)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sat ~ 1 + diffflow + aveflow + (1 | groupid)
## Data: dat
##
## AIC BIC logLik deviance df.resid
## 2732.9 2757.5 -1361.5 2722.9 995
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5325 -0.6194 0.0081 0.6332 3.5781
##
## Random effects:
## Groups Name Variance Std.Dev.
## groupid (Intercept) 0.08968 0.2995
## Residual 0.81670 0.9037
## Number of obs: 1000, groups: groupid, 200
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.08953 0.03580 2.501
## diffflow 0.28837 0.04004 7.203
## aveflow 0.46384 0.05144 9.018
##
## Correlation of Fixed Effects:
## (Intr) dffflw
## diffflow 0.000
## aveflow 0.115 0.000
Users can use lavaan
to run MMC too. Note that
lavaan
cannot run restricted maximum likelihood
(REML
). Full information likelihood is used in both
lavaan
and lme4
packages.
model2 <- '
level: 1
sat ~ diffflow
level: 2
sat ~ aveflow
'
fit2 <- sem(model = model2, data = dat, cluster = "groupid")
summary(fit2)
## lavaan 0.6.17 ended normally after 12 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 5
##
## Number of observations 1000
## Number of clusters [groupid] 200
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
##
## Level 1 [within]:
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## sat ~
## diffflow 0.288 0.040 7.203 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .sat 0.817 0.041 20.000 0.000
##
##
## Level 2 [groupid]:
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## sat ~
## aveflow 0.464 0.051 9.018 0.000
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|)
## .sat 0.090 0.036 2.501 0.012
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .sat 0.090 0.027 3.373 0.001
The regression coefficients from MLC and MMC were different. In this case, group flow is the property of groups so MLC is more appropriate.
Emotional stability is the property of individuals so it represents the formative measurement model. MMC is more appropriate, especially, when the sampling ratio is close to 100% like in this example. The differences between MLC and MMC are shown here. First, MLC is analyzed.
model3 <- '
level: 1
sat ~ emostability
level: 2
sat ~ emostability
'
fit3 <- sem(model = model3, data = dat, cluster = "groupid")
summary(fit3)
## lavaan 0.6.17 ended normally after 21 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 5
##
## Number of observations 1000
## Number of clusters [groupid] 200
##
## Model Test User Model:
##
## Test statistic 0.009
## Degrees of freedom 0
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
##
## Level 1 [within]:
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## sat ~
## emostability 0.238 0.032 7.444 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .sat 0.813 0.041 20.000 0.000
##
##
## Level 2 [groupid]:
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## sat ~
## emostability -0.238 0.787 -0.302 0.762
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|)
## .sat 0.049 0.044 1.133 0.257
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .sat 0.181 0.040 4.518 0.000
Next, MMC is analyzed by the lme4
package.
dat$aveemo <- ave(dat$emostability, dat$groupid)
dat$diffemo <- dat$emostability - dat$aveemo
out4 <- lmer(sat ~ 1 + diffemo + aveemo + (1|groupid), data=dat, REML=FALSE)
summary(out4)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sat ~ 1 + diffemo + aveemo + (1 | groupid)
## Data: dat
##
## AIC BIC logLik deviance df.resid
## 2793.6 2818.2 -1391.8 2783.6 995
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4474 -0.6601 -0.0129 0.6222 3.2054
##
## Random effects:
## Groups Name Variance Std.Dev.
## groupid (Intercept) 0.1858 0.4311
## Residual 0.8133 0.9018
## Number of obs: 1000, groups: groupid, 200
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.05597 0.04178 1.340
## diffemo 0.23776 0.03194 7.444
## aveemo 0.18174 0.08829 2.058
##
## Correlation of Fixed Effects:
## (Intr) diffem
## diffemo 0.000
## aveemo 0.040 0.000
MMC can be analyzed by the lavaan
package as well.
model4 <- '
level: 1
sat ~ diffemo
level: 2
sat ~ aveemo
'
fit4 <- sem(model = model4, data = dat, cluster = "groupid")
summary(fit4)
## lavaan 0.6.17 ended normally after 14 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 5
##
## Number of observations 1000
## Number of clusters [groupid] 200
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
##
## Level 1 [within]:
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## sat ~
## diffemo 0.238 0.032 7.444 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .sat 0.813 0.041 20.000 0.000
##
##
## Level 2 [groupid]:
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## sat ~
## aveemo 0.182 0.088 2.058 0.040
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|)
## .sat 0.056 0.042 1.340 0.180
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .sat 0.186 0.036 5.193 0.000
Similary, the results from MLC and MMC were different but, for emotional stability, MMC results are more accurate.
Because lavaan
is the package for MSEM, users can run
path analysis in each level. In this example, the effects of counselor’s
perceived counseling proficiency on counseling satisfaction is
investigated. Perceived group flow could be the mediator of the effect
between the perceived counseling proficiency and counseling
satisfaction.
Similary, ~
represents the regression effect.
*
in front of each predictor indicates that the names in
front of *
are the labels of each regression coefficient.
The labels will be used at the end of script where labels are multiplied
by each other. The multiplication is used to create the products of
regression coefficients, i.e., indirect effects. The results of the
summary
function show the significance of each regression
coefficients and the indirect effects. Unfortunately, nonparametric
bootstrap cannot be run in multilevel models; therefore, the delta
method, which assume symmetric sampling distribution, is used here.
model5 <- '
level: 1
groupflow ~ aw*perceivedability
sat ~ bw*groupflow + perceivedability
level: 2
groupflow ~ ab*perceivedability + experience
sat ~ bb*groupflow + perceivedability + experience
abw := aw*bw
abb := ab*bb
'
fit5 <- sem(model = model5, data = dat, cluster = "groupid")
summary(fit5)
## lavaan 0.6.17 ended normally after 23 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 14
##
## Number of observations 1000
## Number of clusters [groupid] 200
##
## Model Test User Model:
##
## Test statistic 0.081
## Degrees of freedom 0
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
##
## Level 1 [within]:
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## groupflow ~
## prcvdblty (aw) 0.525 0.036 14.591 0.000
## sat ~
## groupflow (bw) 0.145 0.044 3.323 0.001
## prcvdblty 0.358 0.050 7.153 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .groupflow 0.503 0.025 20.000 0.000
## .sat 0.768 0.038 20.000 0.000
##
##
## Level 2 [groupid]:
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## groupflow ~
## prcvdblty (ab) 0.235 0.089 2.635 0.008
## experienc 0.175 0.054 3.253 0.001
## sat ~
## groupflow (bb) 0.450 0.084 5.327 0.000
## prcvdblty -0.103 0.074 -1.380 0.168
## experienc 0.158 0.047 3.366 0.001
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|)
## .groupflow -0.052 0.043 -1.207 0.227
## .sat 0.090 0.035 2.556 0.011
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .groupflow 0.255 0.037 6.929 0.000
## .sat 0.068 0.025 2.675 0.007
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|)
## abw 0.076 0.024 3.240 0.001
## abb 0.106 0.046 2.283 0.022
See the lecture for the interpretation of the results.
lavaan
is the MSEM package so it can analyze covariances
among variables as well. ~~
represents relationships among
variables. For example,
emostability ~~ groupflow + perceivedability
means that the
covariance between emotional stability and group flows and the
covariance between emotional stability and perceived proficiency are
estimated. In this example, level-2 variable (experience
)
is included too. This variable will show in level: 2
only.
model6 <- '
level: 1
sat ~~ emostability + groupflow + perceivedability
emostability ~~ groupflow + perceivedability
groupflow ~~ perceivedability
level: 2
sat ~~ emostability + groupflow + perceivedability + experience
emostability ~~ groupflow + perceivedability + experience
groupflow ~~ perceivedability + experience
perceivedability ~~ experience
'
fit6 <- sem(model = model6, data = dat, cluster = "groupid")
summary(fit6, standardize=TRUE)
## lavaan 0.6.17 ended normally after 32 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 30
##
## Number of observations 1000
## Number of clusters [groupid] 200
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
##
## Level 1 [within]:
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## sat ~~
## emostability 0.237 0.034 6.976 0.000 0.237 0.255
## groupflow 0.184 0.027 6.777 0.000 0.184 0.247
## perceivedablty 0.211 0.024 8.725 0.000 0.211 0.324
## emostability ~~
## groupflow 0.212 0.029 7.282 0.000 0.212 0.266
## perceivedablty 0.073 0.025 2.948 0.003 0.073 0.105
## groupflow ~~
## perceivedablty 0.255 0.022 11.788 0.000 0.255 0.458
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## sat 0.870 0.043 20.000 0.000 0.870 1.000
## emostability 0.997 0.050 20.000 0.000 0.997 1.000
## groupflow 0.637 0.032 20.000 0.000 0.637 1.000
## perceivedablty 0.486 0.024 20.000 0.000 0.486 1.000
##
##
## Level 2 [groupid]:
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## sat ~~
## emostability -0.007 0.021 -0.318 0.750 -0.007 -0.102
## groupflow 0.185 0.034 5.514 0.000 0.185 0.732
## perceivedablty 0.105 0.034 3.090 0.002 0.105 0.356
## experience 0.263 0.048 5.461 0.000 0.263 0.586
## emostability ~~
## groupflow 0.008 0.024 0.328 0.743 0.008 0.086
## perceivedablty -0.017 0.026 -0.642 0.521 -0.017 -0.154
## experience -0.012 0.035 -0.344 0.731 -0.012 -0.074
## groupflow ~~
## perceivedablty 0.189 0.041 4.609 0.000 0.189 0.462
## experience 0.296 0.056 5.337 0.000 0.296 0.476
## perceivedability ~~
## experience 0.438 0.064 6.813 0.000 0.438 0.603
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## sat 0.053 0.042 1.246 0.213 0.053 0.123
## emostability -0.019 0.033 -0.562 0.574 -0.019 -0.121
## groupflow -0.080 0.049 -1.630 0.103 -0.080 -0.135
## perceivedablty -0.077 0.054 -1.437 0.151 -0.077 -0.111
## experience -0.052 0.074 -0.693 0.488 -0.052 -0.049
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## sat 0.182 0.037 4.967 0.000 0.182 1.000
## emostability 0.024 0.024 0.990 0.322 0.024 1.000
## groupflow 0.351 0.048 7.272 0.000 0.351 1.000
## perceivedablty 0.478 0.058 8.282 0.000 0.478 1.000
## experience 1.105 0.111 10.000 0.000 1.105 1.000
See the right column for the correlation (standardized covariances).
The correlation can be tested whether they are significantly different
from 0 by the standardizedsolution
function.
standardizedsolution(fit6)
## lhs op rhs est.std se z pvalue ci.lower
## 1 sat ~~ emostability 0.255 0.033 7.697 0.000 0.190
## 2 sat ~~ groupflow 0.247 0.033 7.433 0.000 0.182
## 3 sat ~~ perceivedability 0.324 0.032 10.251 0.000 0.262
## 4 emostability ~~ groupflow 0.266 0.033 8.111 0.000 0.202
## 5 emostability ~~ perceivedability 0.105 0.035 2.998 0.003 0.036
## 6 groupflow ~~ perceivedability 0.458 0.028 16.419 0.000 0.404
## 7 sat ~~ sat 1.000 0.000 NA NA 1.000
## 8 emostability ~~ emostability 1.000 0.000 NA NA 1.000
## 9 groupflow ~~ groupflow 1.000 0.000 NA NA 1.000
## 10 perceivedability ~~ perceivedability 1.000 0.000 NA NA 1.000
## 11 sat ~1 0.000 0.000 NA NA 0.000
## 12 emostability ~1 0.000 0.000 NA NA 0.000
## 13 groupflow ~1 0.000 0.000 NA NA 0.000
## 14 perceivedability ~1 0.000 0.000 NA NA 0.000
## 15 sat ~~ emostability -0.102 0.338 -0.302 0.763 -0.764
## 16 sat ~~ groupflow 0.732 0.078 9.451 0.000 0.581
## 17 sat ~~ perceivedability 0.356 0.097 3.668 0.000 0.166
## 18 sat ~~ experience 0.586 0.079 7.389 0.000 0.430
## 19 emostability ~~ groupflow 0.086 0.254 0.338 0.735 -0.412
## 20 emostability ~~ perceivedability -0.154 0.253 -0.610 0.542 -0.650
## 21 emostability ~~ experience -0.074 0.217 -0.341 0.733 -0.499
## 22 groupflow ~~ perceivedability 0.462 0.072 6.427 0.000 0.321
## 23 groupflow ~~ experience 0.476 0.067 7.100 0.000 0.344
## 24 perceivedability ~~ experience 0.603 0.052 11.564 0.000 0.501
## 25 sat ~~ sat 1.000 0.000 NA NA 1.000
## 26 emostability ~~ emostability 1.000 0.000 NA NA 1.000
## 27 groupflow ~~ groupflow 1.000 0.000 NA NA 1.000
## 28 perceivedability ~~ perceivedability 1.000 0.000 NA NA 1.000
## 29 experience ~~ experience 1.000 0.000 NA NA 1.000
## 30 sat ~1 0.123 0.100 1.236 0.216 -0.072
## 31 emostability ~1 -0.121 0.223 -0.540 0.589 -0.558
## 32 groupflow ~1 -0.135 0.083 -1.620 0.105 -0.297
## 33 perceivedability ~1 -0.111 0.078 -1.431 0.152 -0.264
## 34 experience ~1 -0.049 0.071 -0.692 0.489 -0.188
## ci.upper
## 1 0.319
## 2 0.312
## 3 0.386
## 4 0.331
## 5 0.173
## 6 0.513
## 7 1.000
## 8 1.000
## 9 1.000
## 10 1.000
## 11 0.000
## 12 0.000
## 13 0.000
## 14 0.000
## 15 0.560
## 16 0.884
## 17 0.546
## 18 0.741
## 19 0.584
## 20 0.341
## 21 0.351
## 22 0.602
## 23 0.607
## 24 0.705
## 25 1.000
## 26 1.000
## 27 1.000
## 28 1.000
## 29 1.000
## 30 0.319
## 31 0.317
## 32 0.028
## 33 0.041
## 34 0.090