Detailed guide to setting priors
James Uanhoro
priors.RmdDefault priors
We begin by loading the package:
## ## ################################################################################# This is minorbsem 0.2.14##
## All users of R (or SEM) are invited to report bugs, submit functions or ideas
## for functions. An efficient way to do this is to open an issue on GitHub
## https://github.com/jamesuanhoro/minorbsem/issues/.## ###############################################################################The default priors are:
priors_default <- new_mbsempriors()
priors_default## An object of class "mbsempriors"
## Slot "lkj_shape":
## [1] 2
##
## Slot "ml_par":
## [1] 0
##
## Slot "sl_par":
## [1] 1
##
## Slot "rs_par":
## [1] 1
##
## Slot "rc_par":
## [1] 2
##
## Slot "sc_par":
## [1] 0.5
##
## Slot "rm_par":
## [1] 0.15
lkj_shape
The lkj_shape parameter is for setting priors on the
interfactor correlation matrix. By default this value is 2, and implies
each correlation in the matrix follows the distribution: \(\frac{r + 1}{2} \sim \text{Beta}(a, a)\),
where \(a = \eta + \frac{d - 2}{2}\),
where \(\eta\) is the
lkj_shape parameter and \(d\) is the number of factors. For example,
for three factors and \(\eta = 5\),
\(a = 5 + (3 - 2) / 2 = 5.5\), so each
correlation follows the distribution: Beta(5.5, 5.5). Assuming \(\eta = 20\), we have Beta(20.5, 20.5). Both
distributions are visualized below.
We see that larger values of lkj_shape can be used to
constrain interfactor correlations.
ml_par and sl_par
ml_par and sl_par are the prior mean and SD
for loadings which are assumed to be normally distributed. Note that
factors are assumed to be standardized.
rs_par
Residual standard deviation parameters are assumed to follow a
location-scale Student-t distribution with degrees of freedom 3
(unchangeable), location 0 (unchangeable) and scale
rs_par.
rc_par
Residual correlations between indicators assumed \(\text{Beta}(a, a)\) where \(a = \mathsf{rc\_par}\).
sc_par
All coefficients in minorbsem are assumed to be
standardized. These coefficients are assumed normal with mean 0
(changeable, see section on latent
regression below) and SD of sc_par.
Simple changes to default priors
Holzinger-Swineford example
This dataset is discussed in the CFA tutorial.
## x1 x2 x3 x4 x5 x6 x7 x8 x9
## 1 3.333333 7.75 0.375 2.333333 5.75 1.2857143 3.391304 5.75 6.361111
## 2 5.333333 5.25 2.125 1.666667 3.00 1.2857143 3.782609 6.25 7.916667
## 3 4.500000 5.25 1.875 1.000000 1.75 0.4285714 3.260870 3.90 4.416667
## 4 5.333333 7.75 3.000 2.666667 4.50 2.4285714 3.000000 5.30 4.861111
## 5 4.833333 4.75 0.875 2.666667 4.00 2.5714286 3.695652 6.30 5.916667
## 6 5.333333 5.00 2.250 1.000000 3.00 0.8571429 4.347826 6.65 7.500000The scale of data is important for setting priors on model
parameters. The default priors for models fit with
minorbsem are reasonable when variables have standard
deviations close to 1. For this reason, we first check the standard
deviations of the relevant variables for this analysis:
apply(item_data, 2, sd) # compute SD of each variable## x1 x2 x3 x4 x5 x6 x7 x8 x9
## 1.167432 1.177451 1.130979 1.164116 1.290472 1.095603 1.089534 1.012615 1.009152All variables have standard deviations close to 1, so we can move forward with the data as they are. Otherwise, we would recommend re-scaling the variables.1
syntax_basic <- "
Visual =~ x1 + x2 + x3
Verbal =~ x4 + x5 + x6
Speed =~ x7 + x8 + x9"Let’s fit the model with modified priors:
priors_modified <- new_mbsempriors(lkj_shape = 10, sl_par = .75)
priors_modified## An object of class "mbsempriors"
## Slot "lkj_shape":
## [1] 10
##
## Slot "ml_par":
## [1] 0
##
## Slot "sl_par":
## [1] 0.75
##
## Slot "rs_par":
## [1] 1
##
## Slot "rc_par":
## [1] 2
##
## Slot "sc_par":
## [1] 0.5
##
## Slot "rm_par":
## [1] 0.15We can fit make the call to the minorbsem() function and
obtain the data list passed to Stan instead of fitting the model by
setting ret_data_list = TRUE:
base_dl <- minorbsem(
syntax_basic, item_data,
priors = priors_modified, ret_data_list = TRUE
)## Processing user input ...
base_dl$shape_phi_c # reflects the shape of interfactor correlation matrix## [1] 10
base_dl$load_est # reflects prior mean for loadings## Visual Verbal Speed
## x1 0 0 0
## x2 0 0 0
## x3 0 0 0
## x4 0 0 0
## x5 0 0 0
## x6 0 0 0
## x7 0 0 0
## x8 0 0 0
## x9 0 0 0
base_dl$load_se # reflects prior SD for loadings## Visual Verbal Speed
## x1 0.75 0.75 0.75
## x2 0.75 0.75 0.75
## x3 0.75 0.75 0.75
## x4 0.75 0.75 0.75
## x5 0.75 0.75 0.75
## x6 0.75 0.75 0.75
## x7 0.75 0.75 0.75
## x8 0.75 0.75 0.75
## x9 0.75 0.75 0.75
base_dl$loading_pattern # shows indices for non-zero loadings## Visual Verbal Speed
## x1 1 0 0
## x2 2 0 0
## x3 3 0 0
## x4 0 4 0
## x5 0 5 0
## x6 0 6 0
## x7 0 0 7
## x8 0 0 8
## x9 0 0 9We can either run the model with:
base_mod <- minorbsem(syntax_basic, item_data, priors = priors_modified)or
base_mod <- minorbsem(data_list = base_dl, priors = priors_modified)The results would be identical. Adding
priors = priors_modified is optional.2
Examples of priors on specific loadings and coefficients
Small variance prior on all cross-loadings
Following Muthén and Asparouhov (2012),
one can estimate all cross-loadings by placing small variance priors on
them. Note that minorbsem contains a cleaner approach to the same task
with models of the form:
minorbsem(..., simple_struc = FALSE). But we demonstrate
the small variance approach to show how to set specific priors on some
loadings.
syntax_cross <- paste(
paste0("Visual =~ ", paste0("x", 1:9, collapse = " + ")),
paste0("Verbal =~ ", paste0("x", 1:9, collapse = " + ")),
paste0("Speed =~ ", paste0("x", 1:9, collapse = " + ")),
sep = "\n"
)
writeLines(syntax_cross)## Visual =~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9
## Verbal =~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9
## Speed =~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9Then we call minorbsem() with
ret_data_list = TRUE to obtain the data list:
small_var_dl <- minorbsem(syntax_cross, item_data, ret_data_list = TRUE)## Processing user input ...Set all cross-loadings to have prior SDs of 0.1:
small_var_dl$load_se## Visual Verbal Speed
## x1 1 1 1
## x2 1 1 1
## x3 1 1 1
## x4 1 1 1
## x5 1 1 1
## x6 1 1 1
## x7 1 1 1
## x8 1 1 1
## x9 1 1 1
small_var_dl$load_se[4:9, 1] <- 0.1
small_var_dl$load_se[-c(4:6), 2] <- 0.1
small_var_dl$load_se[1:6, 3] <- 0.1
small_var_dl$load_se## Visual Verbal Speed
## x1 1.0 0.1 0.1
## x2 1.0 0.1 0.1
## x3 1.0 0.1 0.1
## x4 0.1 1.0 0.1
## x5 0.1 1.0 0.1
## x6 0.1 1.0 0.1
## x7 0.1 0.1 1.0
## x8 0.1 0.1 1.0
## x9 0.1 0.1 1.0Before fitting the model, minorbsem automatically identifies an indicator per factor that it uses to align the factor in the right direction. Given that all indicators reflect all factors, we would need to identify these indicators manually:
small_var_dl$markers## [1] 1 1 1
small_var_dl$markers <- c(1, 4, 7) # for visual, verbal and speed respectivelyThen fit the model:
small_var_mod <- minorbsem(data_list = small_var_dl)## Parameter estimates (method = normal, sample size(s) = 301)
## from op to mean sd 5.000% 95.000% rhat ess_bulk
## ───────────────────────────────────────────────────────────────────────────────
## Goodness of fit
## ───────────────────────────────────────────────────────────────────────────────
## PPP 0.450 1.000 2174
## RMSE 0.027 0.013 0.007 0.050 1.003 563
## ───────────────────────────────────────────────────────────────────────────────
## Factor loadings
## ───────────────────────────────────────────────────────────────────────────────
## Visual =~ x1 0.757 0.123 0.571 0.969 1.001 1547
## Visual =~ x2 0.557 0.102 0.387 0.728 1.001 3157
## Visual =~ x3 0.755 0.113 0.572 0.939 1.003 1369
## Visual =~ x4 0.036 0.073 -0.085 0.154 1.000 3283
## Visual =~ x5 -0.063 0.077 -0.188 0.062 1.002 3344
## Visual =~ x6 0.069 0.070 -0.047 0.178 1.000 3530
## Visual =~ x7 -0.142 0.082 -0.278 -0.008 1.000 3543
## Visual =~ x8 0.005 0.084 -0.135 0.145 1.001 2502
## Visual =~ x9 0.212 0.076 0.084 0.334 1.000 2433
## Verbal =~ x1 0.123 0.081 -0.014 0.253 1.001 2477
## Verbal =~ x2 0.014 0.072 -0.105 0.130 1.000 4317
## Verbal =~ x3 -0.081 0.075 -0.207 0.039 1.005 3759
## Verbal =~ x4 0.980 0.076 0.857 1.109 1.001 2956
## Verbal =~ x5 1.139 0.085 1.005 1.285 1.001 2546
## Verbal =~ x6 0.887 0.072 0.773 1.012 1.001 2783
## Verbal =~ x7 0.024 0.075 -0.099 0.151 1.002 3216
## Verbal =~ x8 -0.035 0.075 -0.159 0.090 1.001 3352
## Verbal =~ x9 0.032 0.065 -0.075 0.141 1.000 3343
## Speed =~ x1 0.049 0.081 -0.088 0.176 1.000 4136
## Speed =~ x2 -0.050 0.073 -0.175 0.068 1.001 4151
## Speed =~ x3 0.040 0.081 -0.097 0.168 1.002 3184
## Speed =~ x4 0.002 0.071 -0.113 0.117 1.000 4014
## Speed =~ x5 0.007 0.076 -0.121 0.133 1.000 3471
## Speed =~ x6 0.003 0.066 -0.104 0.115 1.001 3233
## Speed =~ x7 0.733 0.103 0.562 0.900 1.003 1732
## Speed =~ x8 0.777 0.101 0.620 0.950 1.008 1298
## Speed =~ x9 0.541 0.089 0.397 0.694 1.001 2399
## ───────────────────────────────────────────────────────────────────────────────
## Inter-factor correlations
## ───────────────────────────────────────────────────────────────────────────────
## Verbal ~~ Visual 0.353 0.106 0.176 0.520 1.001 2688
## Speed ~~ Visual 0.327 0.131 0.099 0.530 1.001 2675
## Speed ~~ Verbal 0.242 0.118 0.039 0.425 1.000 2649
## ───────────────────────────────────────────────────────────────────────────────
## Residual variances
## ───────────────────────────────────────────────────────────────────────────────
## x1 ~~ x1 0.662 0.143 0.427 0.872 1.001 1457
## x2 ~~ x2 1.089 0.115 0.905 1.283 1.001 3433
## x3 ~~ x3 0.725 0.133 0.503 0.929 1.004 1447
## x4 ~~ x4 0.375 0.084 0.236 0.510 1.002 2294
## x5 ~~ x5 0.421 0.111 0.237 0.607 1.001 1620
## x6 ~~ x6 0.372 0.072 0.255 0.487 1.001 2424
## x7 ~~ x7 0.698 0.124 0.489 0.892 1.003 1529
## x8 ~~ x8 0.449 0.121 0.233 0.627 1.013 1056
## x9 ~~ x9 0.583 0.079 0.454 0.712 1.003 2497
## ───────────────────────────────────────────────────────────────────────────────
##
##
## Column names: from, op, to, mean, sd, 5%, 95%, rhat, ess_bulkPriors on specific coefficients
Consider the latent regression model:
syntax_reg <- "
Visual =~ x1 + x2 + x3 + x9
Verbal =~ x4 + x5 + x6
Speed =~ x7 + x8 + x9
Verbal ~ Visual + Speed"Assume we believe the most likely estimate of the coefficient from visual to verbal is .5, we can change the prior mean for this coefficient to .5. We first retrieve the data list object.
reg_dl <- minorbsem(
syntax_reg, item_data,
ret_data_list = TRUE
)## Processing user input ...
reg_dl$coef_pattern # coefficient pattern, note non-zero cells## Visual Verbal Speed
## Visual 0 0 0
## Verbal 1 0 2
## Speed 0 0 0
reg_dl$coef_est # default prior means## Visual Verbal Speed
## Visual 0 0 0
## Verbal 0 0 0
## Speed 0 0 0
reg_dl$coef_est["Verbal", "Visual"] <- .5
reg_dl$coef_est # updated## Visual Verbal Speed
## Visual 0.0 0 0
## Verbal 0.5 0 0
## Speed 0.0 0 0
reg_dl$coef_se # default prior SDs## Visual Verbal Speed
## Visual 0.5 0.5 0.5
## Verbal 0.5 0.5 0.5
## Speed 0.5 0.5 0.5Then we fit the model with the updated data list:
reg_mod <- minorbsem(data_list = reg_dl)## Parameter estimates (method = normal, sample size(s) = 301)
## from op to mean sd 5.000% 95.000% rhat ess_bulk
## ──────────────────────────────────────────────────────────────────────────────
## Goodness of fit
## ──────────────────────────────────────────────────────────────────────────────
## PPP 0.386 1.000 1386
## RMSE 0.037 0.014 0.016 0.060 1.015 415
## ──────────────────────────────────────────────────────────────────────────────
## Latent regression coefficients (outcome ~
## predictor)
## ──────────────────────────────────────────────────────────────────────────────
## Verbal ~ Visual 0.419 0.076 0.289 0.540 1.014 989
## Verbal ~ Speed 0.078 0.080 -0.054 0.211 1.003 2289
## ──────────────────────────────────────────────────────────────────────────────
## R square
## ──────────────────────────────────────────────────────────────────────────────
## Visual ~~ Visual 0.000 0.000 0.000 0.000
## Verbal ~~ Verbal 0.209 0.059 0.119 0.310 1.011 1779
## Speed ~~ Speed 0.000 0.000 0.000 0.000
## ──────────────────────────────────────────────────────────────────────────────
## Factor loadings
## ──────────────────────────────────────────────────────────────────────────────
## Visual =~ x1 0.911 0.100 0.747 1.082 1.005 1443
## Visual =~ x2 0.498 0.086 0.357 0.639 1.001 2992
## Visual =~ x3 0.649 0.085 0.513 0.791 1.000 1930
## Visual =~ x9 0.395 0.081 0.263 0.522 1.004 1864
## Verbal =~ x4 0.994 0.073 0.879 1.117 1.001 1788
## Verbal =~ x5 1.081 0.079 0.948 1.206 1.000 2270
## Verbal =~ x6 0.929 0.070 0.819 1.050 1.002 1601
## Speed =~ x9 0.450 0.100 0.312 0.607 1.008 541
## Speed =~ x7 0.644 0.093 0.495 0.802 1.002 1053
## Speed =~ x8 0.842 0.106 0.668 1.018 1.008 576
## ──────────────────────────────────────────────────────────────────────────────
## Inter-factor correlations
## ──────────────────────────────────────────────────────────────────────────────
## Verbal ~~ Visual 0.000 0.000 0.000 0.000
## Speed ~~ Visual 0.275 0.089 0.130 0.417 1.000 2504
## Speed ~~ Verbal 0.000 0.000 0.000 0.000
## ──────────────────────────────────────────────────────────────────────────────
## Residual variances
## ──────────────────────────────────────────────────────────────────────────────
## x1 ~~ x1 0.537 0.157 0.265 0.773 1.002 1160
## x2 ~~ x2 1.143 0.111 0.964 1.328 1.002 3512
## x3 ~~ x3 0.859 0.110 0.683 1.041 1.000 1751
## x9 ~~ x9 0.574 0.101 0.438 0.713 1.004 671
## x4 ~~ x4 0.365 0.100 0.194 0.520 1.001 1144
## x5 ~~ x5 0.495 0.118 0.310 0.692 1.003 1665
## x6 ~~ x6 0.344 0.091 0.189 0.486 1.001 1179
## x7 ~~ x7 0.776 0.116 0.586 0.966 1.003 1161
## x8 ~~ x8 0.324 0.160 0.036 0.573 1.011 472
## ──────────────────────────────────────────────────────────────────────────────
##
##
## Column names: from, op, to, mean, sd, 5%, 95%, rhat, ess_bulkWe also show results from the fit with default priors:
reg_def_priors <- minorbsem(syntax_reg, item_data)## Parameter estimates (method = normal, sample size(s) = 301)
## from op to mean sd 5.000% 95.000% rhat ess_bulk
## ──────────────────────────────────────────────────────────────────────────────
## Goodness of fit
## ──────────────────────────────────────────────────────────────────────────────
## PPP 0.390 1.000 1752
## RMSE 0.037 0.013 0.017 0.058 1.004 641
## ──────────────────────────────────────────────────────────────────────────────
## Latent regression coefficients (outcome ~
## predictor)
## ──────────────────────────────────────────────────────────────────────────────
## Verbal ~ Visual 0.414 0.074 0.293 0.538 1.001 2824
## Verbal ~ Speed 0.083 0.080 -0.051 0.214 1.003 2586
## ──────────────────────────────────────────────────────────────────────────────
## R square
## ──────────────────────────────────────────────────────────────────────────────
## Visual ~~ Visual 0.000 0.000 0.000 0.000
## Verbal ~~ Verbal 0.206 0.058 0.119 0.308 1.000 3445
## Speed ~~ Speed 0.000 0.000 0.000 0.000
## ──────────────────────────────────────────────────────────────────────────────
## Factor loadings
## ──────────────────────────────────────────────────────────────────────────────
## Visual =~ x1 0.914 0.099 0.763 1.088 1.002 1034
## Visual =~ x2 0.498 0.090 0.348 0.646 1.001 2487
## Visual =~ x3 0.646 0.086 0.502 0.789 1.000 2142
## Visual =~ x9 0.397 0.079 0.269 0.526 1.001 2431
## Verbal =~ x4 0.996 0.071 0.882 1.114 1.000 2151
## Verbal =~ x5 1.086 0.077 0.961 1.214 0.999 1944
## Verbal =~ x6 0.926 0.066 0.819 1.038 1.000 2234
## Speed =~ x9 0.444 0.086 0.309 0.592 1.000 1642
## Speed =~ x7 0.646 0.092 0.497 0.804 1.000 1718
## Speed =~ x8 0.844 0.102 0.681 1.016 1.001 1078
## ──────────────────────────────────────────────────────────────────────────────
## Inter-factor correlations
## ──────────────────────────────────────────────────────────────────────────────
## Verbal ~~ Visual 0.000 0.000 0.000 0.000
## Speed ~~ Visual 0.271 0.088 0.125 0.413 1.000 3243
## Speed ~~ Verbal 0.000 0.000 0.000 0.000
## ──────────────────────────────────────────────────────────────────────────────
## Residual variances
## ──────────────────────────────────────────────────────────────────────────────
## x1 ~~ x1 0.531 0.161 0.241 0.772 1.003 934
## x2 ~~ x2 1.140 0.115 0.962 1.337 1.003 3288
## x3 ~~ x3 0.861 0.105 0.692 1.038 1.001 2177
## x9 ~~ x9 0.578 0.079 0.452 0.708 1.000 2327
## x4 ~~ x4 0.361 0.098 0.201 0.516 1.003 1817
## x5 ~~ x5 0.491 0.119 0.299 0.684 1.000 1571
## x6 ~~ x6 0.347 0.087 0.200 0.488 1.002 1668
## x7 ~~ x7 0.774 0.118 0.574 0.962 1.001 1531
## x8 ~~ x8 0.322 0.155 0.037 0.561 1.001 965
## ──────────────────────────────────────────────────────────────────────────────
##
##
## Column names: from, op, to, mean, sd, 5%, 95%, rhat, ess_bulkWorks Cited
Muthén, Bengt, and Tihomir Asparouhov. 2012. “Bayesian Structural
Equation Modeling: A More Flexible Representation of
Substantive Theory.” Psychological Methods 17 (3):
313–35. https://doi.org/10.1037/a0026802.