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.03Next, 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.002The 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$aveflowNext, 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.000Users 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.001The 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.000Next, 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.000MMC 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.000Similary, 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.022See 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.000See 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