At the time of writing, semhelpinghands
has these groups of functions: manipulating parameter estimates tables,
comparing results from different methods, bootstrapping, and …
others.
This article will only introduce very briefly the groups. Some of them will be introduced in details in forthcoming article dedicated to them.
Let’s load the package first, and also load lavaan.
In using lavaan, I prefer reading the output of
parameterEstimates(), which is compact and clear to me. I
sometimes would like to organize the rows and columns in ways meaningful
to a particular research questions. This can certainly be done using
base R or dplyr. However, I am lazy and want to be able to
do things using just one or two functions, with just one or two
arguments.
This is a sample dataset for illustration, dvs_ivs, with
3 predictors (x1, x2, and x3), 3
outcome variables (y1 y2, and
y3), and a group variable (gp).
First a single sample model:
data(dvs_ivs)
mod <-
"
y1 ~ x1 + x2 + x3
y2 ~ x1 + x3
y3 ~ y2 + x2
"
fit <- sem(model = mod,
data = dvs_ivs,
fixed.x = FALSE)The parameter estimates tables:
est <- parameterEstimates(fit)
est
#> lhs op rhs est se z pvalue ci.lower ci.upper
#> 1 y1 ~ x1 0.206 0.087 2.354 0.019 0.034 0.377
#> 2 y1 ~ x2 0.381 0.093 4.100 0.000 0.199 0.563
#> 3 y1 ~ x3 0.162 0.088 1.852 0.064 -0.009 0.334
#> 4 y2 ~ x1 0.149 0.105 1.425 0.154 -0.056 0.354
#> 5 y2 ~ x3 0.230 0.105 2.187 0.029 0.024 0.437
#> 6 y3 ~ y2 0.044 0.106 0.420 0.675 -0.163 0.251
#> 7 y3 ~ x2 0.295 0.121 2.445 0.014 0.059 0.532
#> 8 y1 ~~ y1 0.695 0.098 7.071 0.000 0.502 0.888
#> 9 y2 ~~ y2 1.013 0.143 7.071 0.000 0.732 1.294
#> 10 y3 ~~ y3 1.206 0.171 7.071 0.000 0.872 1.540
#> 11 y1 ~~ y3 -0.027 0.092 -0.292 0.770 -0.206 0.153
#> 12 x1 ~~ x1 0.926 0.131 7.071 0.000 0.669 1.183
#> 13 x1 ~~ x2 0.118 0.088 1.339 0.180 -0.055 0.291
#> 14 x1 ~~ x3 0.000 0.092 -0.001 0.999 -0.180 0.180
#> 15 x2 ~~ x2 0.828 0.117 7.071 0.000 0.598 1.057
#> 16 x2 ~~ x3 0.073 0.087 0.840 0.401 -0.098 0.244
#> 17 x3 ~~ x3 0.912 0.129 7.071 0.000 0.659 1.165A two-sample model is also fitted to the dataset:
This is the parameter estimates table:
est_gp <- parameterEstimates(fit_gp)
est_gp
#> lhs op rhs block group est se z pvalue ci.lower ci.upper
#> 1 y1 ~ x1 1 1 0.210 0.149 1.404 0.160 -0.083 0.502
#> 2 y1 ~ x2 1 1 0.243 0.138 1.761 0.078 -0.027 0.513
#> 3 y1 ~ x3 1 1 0.246 0.144 1.710 0.087 -0.036 0.527
#> 4 y2 ~ x1 1 1 0.036 0.179 0.199 0.842 -0.315 0.386
#> 5 y2 ~ x3 1 1 0.217 0.176 1.232 0.218 -0.128 0.561
#> 6 y3 ~ y2 1 1 0.096 0.120 0.802 0.423 -0.139 0.332
#> 7 y3 ~ x2 1 1 0.372 0.135 2.754 0.006 0.107 0.637
#> 8 y1 ~~ y1 1 1 0.754 0.157 4.796 0.000 0.446 1.062
#> 9 y2 ~~ y2 1 1 1.148 0.239 4.796 0.000 0.679 1.617
#> 10 y3 ~~ y3 1 1 0.787 0.164 4.796 0.000 0.466 1.109
#> 11 y1 ~~ y3 1 1 0.010 0.114 0.085 0.932 -0.213 0.232
#> 12 x1 ~~ x1 1 1 0.798 0.166 4.796 0.000 0.472 1.125
#> 13 x1 ~~ x2 1 1 0.223 0.132 1.694 0.090 -0.035 0.481
#> 14 x1 ~~ x3 1 1 0.121 0.121 1.001 0.317 -0.116 0.358
#> 15 x2 ~~ x2 1 1 0.938 0.196 4.796 0.000 0.555 1.321
#> 16 x2 ~~ x3 1 1 0.139 0.131 1.061 0.289 -0.118 0.397
#> 17 x3 ~~ x3 1 1 0.825 0.172 4.796 0.000 0.488 1.162
#> 18 y1 ~1 1 1 -0.295 0.137 -2.148 0.032 -0.564 -0.026
#> 19 y2 ~1 1 1 0.027 0.168 0.163 0.870 -0.302 0.357
#> 20 y3 ~1 1 1 0.233 0.131 1.775 0.076 -0.024 0.490
#> 21 x1 ~1 1 1 0.282 0.132 2.141 0.032 0.024 0.540
#> 22 x2 ~1 1 1 -0.067 0.143 -0.469 0.639 -0.347 0.213
#> 23 x3 ~1 1 1 -0.123 0.134 -0.915 0.360 -0.385 0.140
#> 24 y1 ~ x1 2 2 0.258 0.102 2.530 0.011 0.058 0.458
#> 25 y1 ~ x2 2 2 0.484 0.121 3.998 0.000 0.247 0.721
#> 26 y1 ~ x3 2 2 0.121 0.104 1.165 0.244 -0.082 0.324
#> 27 y2 ~ x1 2 2 0.226 0.127 1.774 0.076 -0.024 0.475
#> 28 y2 ~ x3 2 2 0.261 0.129 2.016 0.044 0.007 0.514
#> 29 y3 ~ y2 2 2 -0.003 0.162 -0.020 0.984 -0.322 0.315
#> 30 y3 ~ x2 2 2 0.291 0.190 1.527 0.127 -0.082 0.664
#> 31 y1 ~~ y1 2 2 0.568 0.109 5.196 0.000 0.354 0.783
#> 32 y2 ~~ y2 2 2 0.882 0.170 5.196 0.000 0.549 1.214
#> 33 y3 ~~ y3 2 2 1.412 0.272 5.196 0.000 0.880 1.945
#> 34 y1 ~~ y3 2 2 0.050 0.122 0.410 0.682 -0.189 0.289
#> 35 x1 ~~ x1 2 2 1.020 0.196 5.196 0.000 0.635 1.405
#> 36 x1 ~~ x2 2 2 0.042 0.117 0.361 0.718 -0.187 0.271
#> 37 x1 ~~ x3 2 2 -0.099 0.137 -0.719 0.472 -0.367 0.170
#> 38 x2 ~~ x2 2 2 0.722 0.139 5.196 0.000 0.450 0.994
#> 39 x2 ~~ x3 2 2 0.013 0.115 0.110 0.912 -0.212 0.238
#> 40 x3 ~~ x3 2 2 0.985 0.190 5.196 0.000 0.613 1.356
#> 41 y1 ~1 2 2 0.007 0.104 0.068 0.946 -0.196 0.211
#> 42 y2 ~1 2 2 -0.050 0.129 -0.392 0.695 -0.303 0.202
#> 43 y3 ~1 2 2 -0.333 0.163 -2.042 0.041 -0.652 -0.013
#> 44 x1 ~1 2 2 0.101 0.137 0.738 0.460 -0.168 0.371
#> 45 x2 ~1 2 2 0.092 0.116 0.797 0.425 -0.134 0.319
#> 46 x3 ~1 2 2 -0.065 0.135 -0.482 0.630 -0.330 0.200So, these are what semhelpinghands has for now:
add_sig()Despite the controversy over null hypothesis significance testing, we
still need to report them, for now. They are there, but I have to decode
them mentally. Here comes add_sig():
add_sig(est)
#> lhs op rhs est sig se z pvalue ci.lower ci.upper
#> 1 y1 ~ x1 0.206 * 0.087 2.354 0.019 0.034 0.377
#> 2 y1 ~ x2 0.381 *** 0.093 4.100 0.000 0.199 0.563
#> 3 y1 ~ x3 0.162 0.088 1.852 0.064 -0.009 0.334
#> 4 y2 ~ x1 0.149 0.105 1.425 0.154 -0.056 0.354
#> 5 y2 ~ x3 0.230 * 0.105 2.187 0.029 0.024 0.437
#> 6 y3 ~ y2 0.044 0.106 0.420 0.675 -0.163 0.251
#> 7 y3 ~ x2 0.295 * 0.121 2.445 0.014 0.059 0.532
#> 8 y1 ~~ y1 0.695 *** 0.098 7.071 0.000 0.502 0.888
#> 9 y2 ~~ y2 1.013 *** 0.143 7.071 0.000 0.732 1.294
#> 10 y3 ~~ y3 1.206 *** 0.171 7.071 0.000 0.872 1.540
#> 11 y1 ~~ y3 -0.027 0.092 -0.292 0.770 -0.206 0.153
#> 12 x1 ~~ x1 0.926 *** 0.131 7.071 0.000 0.669 1.183
#> 13 x1 ~~ x2 0.118 0.088 1.339 0.180 -0.055 0.291
#> 14 x1 ~~ x3 0.000 0.092 -0.001 0.999 -0.180 0.180
#> 15 x2 ~~ x2 0.828 *** 0.117 7.071 0.000 0.598 1.057
#> 16 x2 ~~ x3 0.073 0.087 0.840 0.401 -0.098 0.244
#> 17 x3 ~~ x3 0.912 *** 0.129 7.071 0.000 0.659 1.165If bootstrapping was used, it also supports using the confidence intervals, which may yield results different from the p-values when bootstrapping confidence intervals are used (but same results in this example):
add_sig(est,
use = c("pvalue", "ci"))
#> lhs op rhs est sig ci se z pvalue ci.lower ci.upper
#> 1 y1 ~ x1 0.206 * Sig 0.087 2.354 0.019 0.034 0.377
#> 2 y1 ~ x2 0.381 *** Sig 0.093 4.100 0.000 0.199 0.563
#> 3 y1 ~ x3 0.162 0.088 1.852 0.064 -0.009 0.334
#> 4 y2 ~ x1 0.149 0.105 1.425 0.154 -0.056 0.354
#> 5 y2 ~ x3 0.230 * Sig 0.105 2.187 0.029 0.024 0.437
#> 6 y3 ~ y2 0.044 0.106 0.420 0.675 -0.163 0.251
#> 7 y3 ~ x2 0.295 * Sig 0.121 2.445 0.014 0.059 0.532
#> 8 y1 ~~ y1 0.695 *** Sig 0.098 7.071 0.000 0.502 0.888
#> 9 y2 ~~ y2 1.013 *** Sig 0.143 7.071 0.000 0.732 1.294
#> 10 y3 ~~ y3 1.206 *** Sig 0.171 7.071 0.000 0.872 1.540
#> 11 y1 ~~ y3 -0.027 0.092 -0.292 0.770 -0.206 0.153
#> 12 x1 ~~ x1 0.926 *** Sig 0.131 7.071 0.000 0.669 1.183
#> 13 x1 ~~ x2 0.118 0.088 1.339 0.180 -0.055 0.291
#> 14 x1 ~~ x3 0.000 0.092 -0.001 0.999 -0.180 0.180
#> 15 x2 ~~ x2 0.828 *** Sig 0.117 7.071 0.000 0.598 1.057
#> 16 x2 ~~ x3 0.073 0.087 0.840 0.401 -0.098 0.244
#> 17 x3 ~~ x3 0.912 *** Sig 0.129 7.071 0.000 0.659 1.165It also works on the standardized solution:
std <- standardizedSolution(fit)
add_sig(std)
#> lhs op rhs est.std sig se z pvalue ci.lower ci.upper
#> 1 y1 ~ x1 0.208 * 0.087 2.394 0.017 0.038 0.378
#> 2 y1 ~ x2 0.364 *** 0.083 4.364 0.000 0.200 0.527
#> 3 y1 ~ x3 0.163 0.087 1.871 0.061 -0.008 0.333
#> 4 y2 ~ x1 0.138 0.096 1.438 0.150 -0.050 0.326
#> 5 y2 ~ x3 0.212 * 0.095 2.236 0.025 0.026 0.397
#> 6 y3 ~ y2 0.041 0.097 0.420 0.674 -0.149 0.231
#> 7 y3 ~ x2 0.237 * 0.094 2.517 0.012 0.053 0.422
#> 8 y1 ~~ y1 0.767 *** 0.074 10.369 0.000 0.622 0.912
#> 9 y2 ~~ y2 0.936 *** 0.047 19.799 0.000 0.844 1.029
#> 10 y3 ~~ y3 0.941 *** 0.046 20.649 0.000 0.852 1.031
#> 11 y1 ~~ y3 -0.029 0.100 -0.292 0.770 -0.225 0.167
#> 12 x1 ~~ x1 1.000 0.000 NA NA 1.000 1.000
#> 13 x1 ~~ x2 0.135 0.098 1.377 0.169 -0.057 0.328
#> 14 x1 ~~ x3 0.000 0.100 -0.001 0.999 -0.196 0.196
#> 15 x2 ~~ x2 1.000 0.000 NA NA 1.000 1.000
#> 16 x2 ~~ x3 0.084 0.099 0.849 0.396 -0.110 0.279
#> 17 x3 ~~ x3 1.000 0.000 NA NA 1.000 1.000Note: But be careful about the interpretation of the
p-values in the standardized solution, which are based on the
delta method. See
vignette("standardizedSolution_boot_ci").
See the help page of add_sig() for other options.
filter_by()Its purpose is very simple, selecting rows by commonly used columns:
operators (op), “dependent variables” (lhs),
and “independent variables” (rhs).
filter_by(est,
op = "~")
#> lhs op rhs est se z pvalue ci.lower ci.upper
#> 1 y1 ~ x1 0.206 0.087 2.354 0.019 0.034 0.377
#> 2 y1 ~ x2 0.381 0.093 4.100 0.000 0.199 0.563
#> 3 y1 ~ x3 0.162 0.088 1.852 0.064 -0.009 0.334
#> 4 y2 ~ x1 0.149 0.105 1.425 0.154 -0.056 0.354
#> 5 y2 ~ x3 0.230 0.105 2.187 0.029 0.024 0.437
#> 6 y3 ~ y2 0.044 0.106 0.420 0.675 -0.163 0.251
#> 7 y3 ~ x2 0.295 0.121 2.445 0.014 0.059 0.532It also supports filtering by group using group labels:
filter_by(est_gp,
op = "~",
group = "gp1",
fit = fit_gp)
#> lhs op rhs block group est se z pvalue ci.lower ci.upper
#> 1 y1 ~ x1 1 1 0.210 0.149 1.404 0.160 -0.083 0.502
#> 2 y1 ~ x2 1 1 0.243 0.138 1.761 0.078 -0.027 0.513
#> 3 y1 ~ x3 1 1 0.246 0.144 1.710 0.087 -0.036 0.527
#> 4 y2 ~ x1 1 1 0.036 0.179 0.199 0.842 -0.315 0.386
#> 5 y2 ~ x3 1 1 0.217 0.176 1.232 0.218 -0.128 0.561
#> 6 y3 ~ y2 1 1 0.096 0.120 0.802 0.423 -0.139 0.332
#> 7 y3 ~ x2 1 1 0.372 0.135 2.754 0.006 0.107 0.637See the help page of filter_by() for other options.
group_by_dvs() and
group_by_ivs()Sometimes I want the conventional iv-column-dv-row or
dv-column-iv-row format. That’s what group_by_dvs() and
group_by_ivs() are for:
group_by_dvs(est)
#> iv est_y1 est_y2 est_y3
#> x1 x1 0.206 0.149 --
#> x2 x2 0.381 -- 0.295
#> x3 x3 0.162 0.230 --
#> y2 y2 -- -- 0.044
group_by_ivs(est)
#> dv est_x1 est_x2 est_x3 est_y2
#> y1 y1 0.206 0.381 0.162 --
#> y2 y2 0.149 -- 0.230 --
#> y3 y3 -- 0.295 -- 0.044They also supports extracting another column:
group_by_dvs(est,
col_name = "pvalue")
#> iv pvalue_y1 pvalue_y2 pvalue_y3
#> x1 x1 0.019 0.154 --
#> x2 x2 0.000 -- 0.014
#> x3 x3 0.064 0.029 --
#> y2 y2 -- -- 0.675
group_by_ivs(est,
col_name = "pvalue")
#> dv pvalue_x1 pvalue_x2 pvalue_x3 pvalue_y2
#> y1 y1 0.019 0.000 0.064 --
#> y2 y2 0.154 -- 0.029 --
#> y3 y3 -- 0.014 -- 0.675See the help page of group_by_dvs() and
group_by_ivs() for other options.
group_by_groups()In multiple sample models, one common task is comparing results
across groups. I wrote group_by_groups() for this task, to
compare results side-by-side:
group_by_groups(est_gp)
#> lhs op rhs est_1 est_2
#> 1 y1 ~ x1 0.210 0.258
#> 2 y1 ~ x2 0.243 0.484
#> 3 y1 ~ x3 0.246 0.121
#> 4 y2 ~ x1 0.036 0.226
#> 5 y2 ~ x3 0.217 0.261
#> 6 y3 ~ x2 0.372 0.291
#> 7 y3 ~ y2 0.096 -0.003
#> 8 x1 ~~ x1 0.798 1.020
#> 9 x1 ~~ x2 0.223 0.042
#> 10 x1 ~~ x3 0.121 -0.099
#> 11 x2 ~~ x2 0.938 0.722
#> 12 x2 ~~ x3 0.139 0.013
#> 13 x3 ~~ x3 0.825 0.985
#> 14 y1 ~~ y1 0.754 0.568
#> 15 y1 ~~ y3 0.010 0.050
#> 16 y2 ~~ y2 1.148 0.882
#> 17 y3 ~~ y3 0.787 1.412
#> 18 x1 ~1 0.282 0.101
#> 19 x2 ~1 -0.067 0.092
#> 20 x3 ~1 -0.123 -0.065
#> 21 y1 ~1 -0.295 0.007
#> 22 y2 ~1 0.027 -0.050
#> 23 y3 ~1 0.233 -0.333It also supports extracting several columns:
group_by_groups(est_gp,
col_names = c("est", "pvalue"))
#> lhs op rhs est_1 est_2 pvalue_1 pvalue_2
#> 1 y1 ~ x1 0.210 0.258 0.160 0.011
#> 2 y1 ~ x2 0.243 0.484 0.078 0.000
#> 3 y1 ~ x3 0.246 0.121 0.087 0.244
#> 4 y2 ~ x1 0.036 0.226 0.842 0.076
#> 5 y2 ~ x3 0.217 0.261 0.218 0.044
#> 6 y3 ~ x2 0.372 0.291 0.006 0.127
#> 7 y3 ~ y2 0.096 -0.003 0.423 0.984
#> 8 x1 ~~ x1 0.798 1.020 0.000 0.000
#> 9 x1 ~~ x2 0.223 0.042 0.090 0.718
#> 10 x1 ~~ x3 0.121 -0.099 0.317 0.472
#> 11 x2 ~~ x2 0.938 0.722 0.000 0.000
#> 12 x2 ~~ x3 0.139 0.013 0.289 0.912
#> 13 x3 ~~ x3 0.825 0.985 0.000 0.000
#> 14 y1 ~~ y1 0.754 0.568 0.000 0.000
#> 15 y1 ~~ y3 0.010 0.050 0.932 0.682
#> 16 y2 ~~ y2 1.148 0.882 0.000 0.000
#> 17 y3 ~~ y3 0.787 1.412 0.000 0.000
#> 18 x1 ~1 0.282 0.101 0.032 0.460
#> 19 x2 ~1 -0.067 0.092 0.639 0.425
#> 20 x3 ~1 -0.123 -0.065 0.360 0.630
#> 21 y1 ~1 -0.295 0.007 0.032 0.946
#> 22 y2 ~1 0.027 -0.050 0.870 0.695
#> 23 y3 ~1 0.233 -0.333 0.076 0.041If the fit object is used, it can print group labels:
group_by_groups(fit_gp,
col_names = c("est", "pvalue"))
#> lhs op rhs est_gp1 est_gp2 pvalue_gp1 pvalue_gp2
#> 1 y1 ~ x1 0.210 0.258 0.160 0.011
#> 2 y1 ~ x2 0.243 0.484 0.078 0.000
#> 3 y1 ~ x3 0.246 0.121 0.087 0.244
#> 4 y2 ~ x1 0.036 0.226 0.842 0.076
#> 5 y2 ~ x3 0.217 0.261 0.218 0.044
#> 6 y3 ~ x2 0.372 0.291 0.006 0.127
#> 7 y3 ~ y2 0.096 -0.003 0.423 0.984
#> 8 x1 ~~ x1 0.798 1.020 0.000 0.000
#> 9 x1 ~~ x2 0.223 0.042 0.090 0.718
#> 10 x1 ~~ x3 0.121 -0.099 0.317 0.472
#> 11 x2 ~~ x2 0.938 0.722 0.000 0.000
#> 12 x2 ~~ x3 0.139 0.013 0.289 0.912
#> 13 x3 ~~ x3 0.825 0.985 0.000 0.000
#> 14 y1 ~~ y1 0.754 0.568 0.000 0.000
#> 15 y1 ~~ y3 0.010 0.050 0.932 0.682
#> 16 y2 ~~ y2 1.148 0.882 0.000 0.000
#> 17 y3 ~~ y3 0.787 1.412 0.000 0.000
#> 18 x1 ~1 0.282 0.101 0.032 0.460
#> 19 x2 ~1 -0.067 0.092 0.639 0.425
#> 20 x3 ~1 -0.123 -0.065 0.360 0.630
#> 21 y1 ~1 -0.295 0.007 0.032 0.946
#> 22 y2 ~1 0.027 -0.050 0.870 0.695
#> 23 y3 ~1 0.233 -0.333 0.076 0.041See the help page of group_by_groups() for other
options.
group_by_models()This is inspired by the proposal Rönkkö proposed in a GitHub issue
for semTools. I want something simple for a quick overview
and so I wrote group_by_models().
Suppose this is the other model fitted:
mod2 <-
"
y1 ~ x1 + x2 + x3
y2 ~ x1 + x2
y3 ~ y2 + x1
"
fit2 <- sem(model = mod2,
data = dvs_ivs,
fixed.x = FALSE)
est2 <- parameterEstimates(fit2)
est2
#> lhs op rhs est se z pvalue ci.lower ci.upper
#> 1 y1 ~ x1 0.209 0.087 2.394 0.017 0.038 0.381
#> 2 y1 ~ x2 0.387 0.093 4.171 0.000 0.205 0.569
#> 3 y1 ~ x3 0.162 0.088 1.853 0.064 -0.009 0.334
#> 4 y2 ~ x1 0.176 0.106 1.655 0.098 -0.032 0.384
#> 5 y2 ~ x2 -0.209 0.112 -1.864 0.062 -0.430 0.011
#> 6 y3 ~ y2 0.023 0.109 0.210 0.834 -0.190 0.236
#> 7 y3 ~ x1 -0.157 0.117 -1.339 0.181 -0.388 0.073
#> 8 y1 ~~ y1 0.695 0.098 7.071 0.000 0.502 0.888
#> 9 y2 ~~ y2 1.026 0.145 7.071 0.000 0.741 1.310
#> 10 y3 ~~ y3 1.254 0.177 7.071 0.000 0.906 1.601
#> 11 y1 ~~ y3 -0.028 0.093 -0.299 0.765 -0.211 0.155
#> 12 x1 ~~ x1 0.926 0.131 7.071 0.000 0.669 1.183
#> 13 x1 ~~ x2 0.118 0.088 1.339 0.180 -0.055 0.291
#> 14 x1 ~~ x3 0.000 0.092 -0.001 0.999 -0.180 0.180
#> 15 x2 ~~ x2 0.828 0.117 7.071 0.000 0.598 1.057
#> 16 x2 ~~ x3 0.073 0.087 0.840 0.401 -0.098 0.244
#> 17 x3 ~~ x3 0.912 0.129 7.071 0.000 0.659 1.165These two models have no nested relationships. To compare the
estimates, group_by_models() can be used:
group_by_models(list(Model1 = est,
Model2 = est2))
#> lhs op rhs est_Model1 est_Model2
#> 1 y1 ~ x1 0.206 0.209
#> 2 y1 ~ x2 0.381 0.387
#> 3 y1 ~ x3 0.162 0.162
#> 4 y2 ~ x1 0.149 0.176
#> 5 y2 ~ x2 NA -0.209
#> 6 y2 ~ x3 0.230 NA
#> 7 y3 ~ x1 NA -0.157
#> 8 y3 ~ x2 0.295 NA
#> 9 y3 ~ y2 0.044 0.023
#> 10 x1 ~~ x1 0.926 0.926
#> 11 x1 ~~ x2 0.118 0.118
#> 12 x1 ~~ x3 0.000 0.000
#> 13 x2 ~~ x2 0.828 0.828
#> 14 x2 ~~ x3 0.073 0.073
#> 15 x3 ~~ x3 0.912 0.912
#> 16 y1 ~~ y1 0.695 0.695
#> 17 y1 ~~ y3 -0.027 -0.028
#> 18 y2 ~~ y2 1.013 1.026
#> 19 y3 ~~ y3 1.206 1.254It can also compare several columns:
group_by_models(list(Model1 = est,
Model2 = est2),
col_names = c("est", "pvalue"))
#> lhs op rhs est_Model1 est_Model2 pvalue_Model1 pvalue_Model2
#> 1 y1 ~ x1 0.206 0.209 0.019 0.017
#> 2 y1 ~ x2 0.381 0.387 0.000 0.000
#> 3 y1 ~ x3 0.162 0.162 0.064 0.064
#> 4 y2 ~ x1 0.149 0.176 0.154 0.098
#> 5 y2 ~ x2 NA -0.209 NA 0.062
#> 6 y2 ~ x3 0.230 NA 0.029 NA
#> 7 y3 ~ x1 NA -0.157 NA 0.181
#> 8 y3 ~ x2 0.295 NA 0.014 NA
#> 9 y3 ~ y2 0.044 0.023 0.675 0.834
#> 10 x1 ~~ x1 0.926 0.926 0.000 0.000
#> 11 x1 ~~ x2 0.118 0.118 0.180 0.180
#> 12 x1 ~~ x3 0.000 0.000 0.999 0.999
#> 13 x2 ~~ x2 0.828 0.828 0.000 0.000
#> 14 x2 ~~ x3 0.073 0.073 0.401 0.401
#> 15 x3 ~~ x3 0.912 0.912 0.000 0.000
#> 16 y1 ~~ y1 0.695 0.695 0.000 0.000
#> 17 y1 ~~ y3 -0.027 -0.028 0.770 0.765
#> 18 y2 ~~ y2 1.013 1.026 0.000 0.000
#> 19 y3 ~~ y3 1.206 1.254 0.000 0.000See the help page of group_by_models() for other
options.
sort_by()The function sort_by() can be used to sort rows using
(a) common fields such as "op", "lhs", and
"rhs", (b) operators such as "~" and
"~~". It can be used on the output of some other functions
that manipulate a parameter estimates table.
out <- group_by_groups(est_gp,
col_names = c("est", "pvalue"))
out <- filter_by(out,
op = c("~", "~~"))
sort_by(out,
by = c("op", "rhs"))
#> lhs op rhs est_1 est_2 pvalue_1 pvalue_2
#> 1 y1 ~ x1 0.210 0.258 0.160 0.011
#> 2 y2 ~ x1 0.036 0.226 0.842 0.076
#> 3 y1 ~ x2 0.243 0.484 0.078 0.000
#> 4 y3 ~ x2 0.372 0.291 0.006 0.127
#> 5 y1 ~ x3 0.246 0.121 0.087 0.244
#> 6 y2 ~ x3 0.217 0.261 0.218 0.044
#> 7 y3 ~ y2 0.096 -0.003 0.423 0.984
#> 8 x1 ~~ x1 0.798 1.020 0.000 0.000
#> 9 x1 ~~ x2 0.223 0.042 0.090 0.718
#> 10 x2 ~~ x2 0.938 0.722 0.000 0.000
#> 11 x1 ~~ x3 0.121 -0.099 0.317 0.472
#> 12 x2 ~~ x3 0.139 0.013 0.289 0.912
#> 13 x3 ~~ x3 0.825 0.985 0.000 0.000
#> 14 y1 ~~ y1 0.754 0.568 0.000 0.000
#> 15 y2 ~~ y2 1.148 0.882 0.000 0.000
#> 16 y1 ~~ y3 0.010 0.050 0.932 0.682
#> 17 y3 ~~ y3 0.787 1.412 0.000 0.000Some functions, such as group_by_models(), will
automatically call sort_by() to sort the results.
The default order should be acceptable in most cases, but can also be
customized. See the help page of sort_by() on customizing
the order.
Though not officially supported, piping using |> can
be used with most of the functions that manipulate parameter estimates
tables. For example:
est_gp |>
add_sig() |>
group_by_groups(col_names = c("est", "pvalue", "sig"),
group_first = FALSE) |>
filter_by(op = c("~"))
#> lhs op rhs est_1 pvalue_1 sig_1 est_2 pvalue_2 sig_2
#> 1 y1 ~ x1 0.210 0.160 0.258 0.011 *
#> 2 y1 ~ x2 0.243 0.078 0.484 0.000 ***
#> 3 y1 ~ x3 0.246 0.087 0.121 0.244
#> 4 y2 ~ x1 0.036 0.842 0.226 0.076
#> 5 y2 ~ x3 0.217 0.218 0.261 0.044 *
#> 6 y3 ~ x2 0.372 0.006 ** 0.291 0.127
#> 7 y3 ~ y2 0.096 0.423 -0.003 0.984Though I believe the choice of the estimation method should be
justified, suppose we want to assess the sensitivity of the parameter
estimate results to the methods used, compare_estimators()
can be used as a quick way to compare results from different
methods.
out <- compare_estimators(fit,
estimator = c("ML", "GLS", "MLR"))
group_by_models(out,
col_names = c("se", "pvalue"))
#> lhs op rhs se_ML se_GLS se_MLR pvalue_ML pvalue_GLS pvalue_MLR
#> 1 y1 ~ x1 0.087 0.094 0.067 0.019 0.028 0.002
#> 2 y1 ~ x2 0.093 0.100 0.091 0.000 0.000 0.000
#> 3 y1 ~ x3 0.088 0.089 0.084 0.064 0.067 0.053
#> 4 y2 ~ x1 0.105 0.112 0.107 0.154 0.148 0.165
#> 5 y2 ~ x3 0.105 0.107 0.123 0.029 0.028 0.061
#> 6 y3 ~ x2 0.121 0.133 0.111 0.014 0.019 0.008
#> 7 y3 ~ y2 0.106 0.116 0.108 0.675 0.654 0.681
#> 8 x1 ~~ x1 0.131 0.129 0.136 0.000 0.000 0.000
#> 9 x1 ~~ x2 0.088 0.090 0.077 0.180 0.202 0.122
#> 10 x1 ~~ x3 0.092 0.092 0.102 0.999 0.840 0.999
#> 11 x2 ~~ x2 0.117 0.114 0.127 0.000 0.000 0.000
#> 12 x2 ~~ x3 0.087 0.089 0.091 0.401 0.414 0.423
#> 13 x3 ~~ x3 0.129 0.131 0.111 0.000 0.000 0.000
#> 14 y1 ~~ y1 0.098 0.100 0.119 0.000 0.000 0.000
#> 15 y1 ~~ y3 0.092 0.093 0.104 0.770 0.793 0.797
#> 16 y2 ~~ y2 0.143 0.139 0.141 0.000 0.000 0.000
#> 17 y3 ~~ y3 0.171 0.167 0.171 0.000 0.000 0.000It simply refits the models for each estimator and returns the
results. They can then be treated as different “models” and processed by
group_by_models().
se_ratios() is a wrapper of
group_by_models() used to compare the standard errors by
different estimators in the output of
compare_estimators():
se_ratios(out,
reference = "ML")
#> lhs op rhs se_ML se_GLS se_MLR ratio_ML ratio_GLS ratio_MLR
#> 1 y1 ~ x1 0.087 0.094 0.067 1 1.072 0.765
#> 2 y1 ~ x2 0.093 0.100 0.091 1 1.073 0.977
#> 3 y1 ~ x3 0.088 0.089 0.084 1 1.013 0.958
#> 4 y2 ~ x1 0.105 0.112 0.107 1 1.068 1.027
#> 5 y2 ~ x3 0.105 0.107 0.123 1 1.011 1.165
#> 6 y3 ~ x2 0.121 0.133 0.111 1 1.101 0.918
#> 7 y3 ~ y2 0.106 0.116 0.108 1 1.099 1.020
#> 8 x1 ~~ x1 0.131 0.129 0.136 1 0.985 1.036
#> 9 x1 ~~ x2 0.088 0.090 0.077 1 1.015 0.867
#> 10 x1 ~~ x3 0.092 0.092 0.102 1 0.999 1.114
#> 11 x2 ~~ x2 0.117 0.114 0.127 1 0.974 1.084
#> 12 x2 ~~ x3 0.087 0.089 0.091 1 1.015 1.049
#> 13 x3 ~~ x3 0.129 0.131 0.111 1 1.012 0.859
#> 14 y1 ~~ y1 0.098 0.100 0.119 1 1.015 1.210
#> 15 y1 ~~ y3 0.092 0.093 0.104 1 1.015 1.134
#> 16 y2 ~~ y2 0.143 0.139 0.141 1 0.973 0.987
#> 17 y3 ~~ y3 0.171 0.167 0.171 1 0.979 1.005See the help page of compare_estimators() for other
options.
One issue with the standardized solution is the confidence intervals.
They are based on the delta method even if se = "boot" is
used. For indirect effects, for which bootstrap confidence intervals are
commonly used, the confidence intervals for them in the standardized
solution are not what usually reported other tools for mediation. There
are some other powerful tools on the Internet and the [Google Group for
lavaan], see this
thread, to address this problem. I wrote
standardizedSolution_boot_ci() not to replace them (it
obviously can’t), but to address a very specific case I usually
encounter myself: Generating the bootstrap confidence intervals for
standardized estimates based on the bootstrap estimates already
generated by se = "boot".
Please see vignette("standardizedSolution_boot_ci") for
an illustration on standardizedSolution_boot_ci().
Another issue is examine bootstrap estimates, such as the
distribution of the estimates. The function plot_boot() and
related functions can be used to compute bootstrap estimates and plot
them. The bootstrap estimates of free parameters (stored by
lavaan), user-defined parameters (computed by
store_boot_def()), and the standardized solution (computed
by store_boot_est_std()) can be plotted.
Please see the article
on plot_boot() on how to use this function.
lavaan::lavaan() and its wrappers suc has
lavaan::sem() and lavaan::cfa() allow users to
set several options using an estimator: ML,
GLS, WLSMV, and others. However, it is not
easy to remember what the options set for an estimator. Instead of
finding them from the output of summary(), this function
shows some of them in one table for a quick overview. This is an
example:
data(dvs_ivs)
mod <-
"
y1 ~ x1 + x2 + x3
y2 ~ x1 + x3
y3 ~ y2 + x2
"
fit_default <- sem(model = mod,
data = dvs_ivs)
show_more_options(fit_default)
#> Options Call Actual
#> Estimator(s) default ML
#> Standard Error (SE) default standard
#> Model Test Statistic(s) default standard
#> How Missing Data is Handled default listwise
#> Information Matrix (for SE) default expected
#> Information Matrix (for Model Test) default expected
#> Mean Structure default No
#> 'x' Fixed default TRUE
fit_MLR <- sem(model = mod,
data = dvs_ivs,
estimator = "MLR")
show_more_options(fit_MLR)
#> Options Call Actual
#> Estimator(s) MLR ML
#> Standard Error (SE) default robust.huber.white
#> Model Test Statistic(s) default yuan.bentler.mplus
#> How Missing Data is Handled default listwise
#> Information Matrix (for SE) default observed
#> Information Matrix (for Model Test) default observed
#> Mean Structure default No
#> 'x' Fixed default TRUE
fit_MLR_fiml <- sem(model = mod,
data = dvs_ivs,
estimator = "MLR",
missing = "fiml")
show_more_options(fit_MLR_fiml)
#> Options Call Actual
#> Estimator(s) MLR ML
#> Standard Error (SE) default robust.huber.white
#> Model Test Statistic(s) default yuan.bentler.mplus
#> How Missing Data is Handled fiml ml
#> Information Matrix (for SE) default observed
#> Information Matrix (for Model Test) default observed
#> Mean Structure default Yes
#> 'x' Fixed default TRUEIn structural equation modeling, closed-form solution is rare and
optimization (minimization) is used to find the/a solution. Out of
curiosity and teaching, I wrote a function to capture the minimization
history so I can examine it or even plot the process. The function,
record_history(), is still in early development but should
work for now for common simple scenarios. Please refer the help page of
record_history() to learn more about it.
In lavaan, I rarely need to manually add covariances
between exogenous variables (defined in a loose sense: variables appear
on the right-hand side but not on the left-hand side of ~).
However, I came into a situation in which lavaan will not
do this (for good reasons). For example, when a covariance between a
residual term and an exogenous variable is set to free. I wrote two
simple functions, add_exo_cov() and
auto_exo_cov(), for this purpose. Please refer to their
help pages for further information.