# use groundhog to make code maximally reproducible
if (!require("groundhog", quietly = TRUE)) {
install.packages("groundhog")
}
## Attached: 'Groundhog' (Version: 3.1.2)
## Tips and troubleshooting: https://groundhogR.com
library("groundhog")
# use groundhog to install and load packages - we may not use all of them and some consider it bad practise
# to still load them, but all of them are nice to have and could be useful
pkgs <- c("here", # System path management
"tidyverse", # ggplot, dplyr, %>%, and friends
"easystats", # bayestestR, performance, effectsize and friends
"modelsummary", # Data and model summaries with tables and plots
"psych", # Descriptive statistics and other tools
"brms", # Bayesian modeling through Stan
"blavaan", # Bayesian SEM
"lavaan", # Frequentist SEM
"semPlot", # Plots for SEM
"tidybayes", # Integration of Bayesian models in tidy + ggplot workflow
"broom", # Convert model objects to data frames
"broom.mixed", # Convert brms model objects to data frames
"tinytable", # Lightweight package to create tables
"hrbrthemes", # Additional ggplot themes
"extrafont", # Additional fonts for plots etc
"ggdist", # Special geoms for posterior distributions
"ggrepel", # Automatically position labels
"gghalves", # Special half geoms
"patchwork" # Combine ggplot objects
)
groundhog.library(pkgs, "2024-03-01") # change to a recent date
## here() starts at C:/Users/Olaf/Documents/olafborghi.github.io/blog/2024/04/16
## Warning: package 'ggplot2' was built under R version 4.3.3
## ── Attaching core tidyverse packages ────────────────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.0 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ──────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Warning: package 'easystats' was built under R version 4.3.3
## # Attaching packages: easystats 0.7.0 (red = needs update)
## ✔ bayestestR 0.13.2 ✔ correlation 0.8.4
## ✖ datawizard 0.9.1 ✖ effectsize 0.8.6
## ✖ insight 0.19.8 ✔ modelbased 0.8.7
## ✖ performance 0.10.9 ✖ parameters 0.21.5
## ✔ report 0.5.8 ✖ see 0.8.2
##
## Restart the R-Session and update packages with `easystats::easystats_update()`.
##
## Version 2.0.0 of `modelsummary`, to be released soon, will introduce a breaking
## change: The default table-drawing package will be `tinytable` instead of
## `kableExtra`. All currently supported table-drawing packages will continue to be
## supported for the foreseeable future, including `kableExtra`, `gt`, `huxtable`,
## `flextable, and `DT`.
##
## You can always call the `config_modelsummary()` function to change the default
## table-drawing package in persistent fashion. To try `tinytable` now:
##
## config_modelsummary(factory_default = 'tinytable')
##
## To set the default back to `kableExtra`:
##
## config_modelsummary(factory_default = 'kableExtra')
##
## Attaching package: 'modelsummary'
##
## The following object is masked from 'package:parameters':
##
## supported_models
##
## The following object is masked from 'package:insight':
##
## supported_models
##
##
## Attaching package: 'psych'
##
## The following object is masked from 'package:modelsummary':
##
## SD
##
## The following object is masked from 'package:effectsize':
##
## phi
##
## The following object is masked from 'package:datawizard':
##
## rescale
##
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
##
## Loading required package: Rcpp
## Loading 'brms' package (version 2.20.4). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
##
## Attaching package: 'brms'
##
## The following object is masked from 'package:psych':
##
## cs
##
## The following object is masked from 'package:stats':
##
## ar
## Warning: package 'blavaan' was built under R version 4.3.3
## This is blavaan 0.5-3
## On multicore systems, we suggest use of future::plan("multicore") or
## future::plan("multisession") for faster post-MCMC computations.
## This is lavaan 0.6-17
## lavaan is FREE software! Please report any bugs.
##
## Attaching package: 'lavaan'
##
## The following object is masked from 'package:psych':
##
## cor2cov
##
##
## Attaching package: 'tidybayes'
##
## The following objects are masked from 'package:brms':
##
## dstudent_t, pstudent_t, qstudent_t, rstudent_t
##
## The following object is masked from 'package:parameters':
##
## parameters
##
## The following object is masked from 'package:bayestestR':
##
## hdi
## Warning: package 'broom.mixed' was built under R version 4.3.3
## NOTE: Either Arial Narrow or Roboto Condensed fonts are required to use these themes.
## Please use hrbrthemes::import_roboto_condensed() to install Roboto Condensed and
## if Arial Narrow is not on your system, please see https://bit.ly/arialnarrow
## Registering fonts with R
##
## Attaching package: 'ggdist'
##
## The following objects are masked from 'package:brms':
##
## dstudent_t, pstudent_t, qstudent_t, rstudent_t
##
## The following object is masked from 'package:bayestestR':
##
## hdi
## Warning: package 'ggrepel' was built under R version 4.3.3
## Warning: package 'gghalves' was built under R version 4.3.3
## Warning: package 'patchwork' was built under R version 4.3.3
## [36mSuccessfully attached 'here_1.0.1'[0m
## [36mSuccessfully attached 'tidyverse_2.0.0'[0m
## [36mSuccessfully attached 'easystats_0.7.0'[0m
## [36mSuccessfully attached 'modelsummary_1.4.5'[0m
## [36mSuccessfully attached 'psych_2.4.1'[0m
## [36mSuccessfully attached 'brms_2.20.4'[0m
## [36mSuccessfully attached 'blavaan_0.5-3'[0m
## [36mSuccessfully attached 'lavaan_0.6-17'[0m
## [36mSuccessfully attached 'semPlot_1.1.6'[0m
## [36mSuccessfully attached 'tidybayes_3.0.6'[0m
## [36mSuccessfully attached 'broom_1.0.5'[0m
## [36mSuccessfully attached 'broom.mixed_0.2.9.4'[0m
## [36mSuccessfully attached 'tinytable_0.0.5'[0m
## [36mSuccessfully attached 'hrbrthemes_0.8.0'[0m
## [36mSuccessfully attached 'extrafont_0.19'[0m
## [36mSuccessfully attached 'ggdist_3.3.1'[0m
## [36mSuccessfully attached 'ggrepel_0.9.5'[0m
## [36mSuccessfully attached 'gghalves_0.1.4'[0m
## [36mSuccessfully attached 'patchwork_1.2.0'[0mTeaching Myself: Bayesian Path Modeling
Using the same dataset as in previous blog posts, this time I will look into path models within a Bayesian statistical framework. Path models allow to directly check complex models, including associations between multiple predictors, multiple outcomes, and indirect associations.
I will not touch upon the topics of causality etc. here. However, I want to point out that before running a path model, one should carefully think about the relationships of interest. In addition, analyses of indirect effects are often called “mediation analysis”. Some use mediation in a merely causal context, and as I don’t touch causality here, I will stick to the term “indirect association”. Path models on cross-sectional data are also criticized from time to time, so here a few suggested reads before you think of doing this in your own research:
- That’s a Lot to Process! Pitfalls of Popular Path Models by Rohrer et al. (2022)
- Thinking Clearly About Correlations and Causation: Graphical Causal Models for Observational Data, again by Rohrer (2018)
- And as many people criticize statistical models that test for indirect associations in cross-sectional data, here a bit of a different view on it.
So beware: This blog post is mainly about messing around with some data and trying out different functions!
Specifically, we want to investigate whether there is an indirect association between Age and Conservatism, Dogmatism, and Religiosity via five cognitive factors Caution, Temporal Discounting, Perceptual Processing Time, Speed of Evidence Accumulation, and Strategic Information Processing (Self-regulation ontology; see Eisenberg et al., 2018; 2019)
Meet the data
I will again use open data from Zmigrod et al. (2021). This dataset includes cognitive and personality factors derived from exploratory factor analysis (EFA) of a large number of behavioral cognitive tasks and personality questionnaires from the Self-Regulation Ontology and in addition to that, several measures of ideology.
Some preparations and explorations
Exploratory factor analysis
I will directly re-use some of the code that was shared together with the manuscript (Zmigrod et al., 2021; see above) to obtain three ideology factors using EFA on the scores of several questionnaires.
# create df with relevant data for EFA
efa_data <- data.frame(
secs_social = data$SocialConservatism,
secs_econ = data$EconConservatism,
nation = data$Nationalism,
auth = data$Authoritarianism,
SJ = data$SJ,
sdo = data$SDO,
patriot = data$Patriotism,
progrSum = data$IngroupAttitude, #
dogma = data$Dogmatism,
IH1 = -1*data$IH1, # IH factors reversed because conceptually inverse of dogmatism
IH2 = -1*data$IH2,
IH3 = -1*data$IH3,
IH4 = -1*data$IH4,
rel_atten = data$ReligiousAttendance,
rel_pray = data$ReligiousPrayer,
rel_imp = data$ReligionImportance
)
# run EFA
scree(efa_data) # scree plot
fa.parallel(efa_data, fm = "ml") # parallel analysis
## Parallel analysis suggests that the number of factors = 3 and the number of components = 3
Ideol_EFA <- fa(efa_data,
nfactors=3,
rotate="oblimin") # EFA oblimin for 3 factors
## Loading required namespace: GPArotation
print.psych(Ideol_EFA)
## Factor Analysis using method = minres
## Call: fa(r = efa_data, nfactors = 3, rotate = "oblimin")
## Standardized loadings (pattern matrix) based upon correlation matrix
## MR1 MR2 MR3 h2 u2 com
## secs_social 0.72 -0.06 0.32 0.76 0.245 1.4
## secs_econ 0.75 -0.08 -0.06 0.51 0.488 1.0
## nation 0.79 0.18 -0.02 0.70 0.299 1.1
## auth 0.44 0.11 0.10 0.28 0.725 1.2
## SJ 0.73 -0.02 -0.08 0.49 0.507 1.0
## sdo 0.59 0.25 -0.13 0.44 0.564 1.5
## patriot 0.76 -0.14 0.06 0.59 0.414 1.1
## progrSum 0.53 0.14 -0.02 0.33 0.668 1.1
## dogma 0.23 0.75 0.05 0.71 0.294 1.2
## IH1 -0.17 0.54 0.13 0.29 0.705 1.3
## IH2 0.01 0.62 0.08 0.40 0.604 1.0
## IH3 -0.12 0.83 -0.07 0.66 0.342 1.1
## IH4 0.07 0.63 0.06 0.44 0.561 1.0
## rel_atten 0.07 0.04 0.72 0.56 0.440 1.0
## rel_pray 0.02 -0.02 0.88 0.78 0.215 1.0
## rel_imp -0.04 0.03 0.99 0.95 0.047 1.0
##
## MR1 MR2 MR3
## SS loadings 3.84 2.54 2.50
## Proportion Var 0.24 0.16 0.16
## Cumulative Var 0.24 0.40 0.56
## Proportion Explained 0.43 0.29 0.28
## Cumulative Proportion 0.43 0.72 1.00
##
## With factor correlations of
## MR1 MR2 MR3
## MR1 1.00 0.23 0.34
## MR2 0.23 1.00 0.08
## MR3 0.34 0.08 1.00
##
## Mean item complexity = 1.1
## Test of the hypothesis that 3 factors are sufficient.
##
## df null model = 120 with the objective function = 8.8 with Chi Square = 2876
## df of the model are 75 and the objective function was 0.92
##
## The root mean square of the residuals (RMSR) is 0.04
## The df corrected root mean square of the residuals is 0.05
##
## The harmonic n.obs is 331 with the empirical chi square 117 with prob < 0.0013
## The total n.obs was 334 with Likelihood Chi Square = 300 with prob < 1.1e-28
##
## Tucker Lewis Index of factoring reliability = 0.868
## RMSEA index = 0.095 and the 90 % confidence intervals are 0.084 0.106
## BIC = -136
## Fit based upon off diagonal values = 0.99
## Measures of factor score adequacy
## MR1 MR2 MR3
## Correlation of (regression) scores with factors 0.95 0.93 0.98
## Multiple R square of scores with factors 0.91 0.86 0.97
## Minimum correlation of possible factor scores 0.82 0.73 0.93Alright, just as in the manuscript the EFA indicates three factors, the first appears to be a factor of political conservatism, as there are high loadings of social and economic conservatism, nationalism, religiosity on this factor. The second appears to be a factor of dogmatism, with high loadings of the dogmatism and intellectual humility measure. The third appears to be a religiosity factor.
If you want to explore EFA a bit more, here again a suggested read Factor Analysis with the psych package by Micheal Clark. For now we don’t really care that much about the details, all we want are the three factor scores so that we can explore a few path models.
# obtain factor loadings and add them to data
factor_scores <- as.data.frame(Ideol_EFA$scores)
data <- add_column(
data,
conservatism_factor = factor_scores$MR1,
dogmatism_factor = factor_scores$MR2,
religiosity_factor = factor_scores$MR3
)Standardizing the data
Again, I will standardize the data. The scales are arbitrary anyway in this case, and this simplifies specifying the priors in this very basic example.
df <- data %>%
select(1, 3:11, 59:61) %>%
rename(ID = SubjectN) %>%
mutate(Edu = as.factor(Edu),
Gender = as.factor(Gender),
ID = as.factor(ID)) %>%
mutate(across(where(is.numeric), scale))
str(df)
## tibble [334 × 13] (S3: tbl_df/tbl/data.frame)
## $ ID : Factor w/ 334 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ Age : num [1:334, 1] -1.776 -0.479 -1.068 -0.715 -0.951 ...
## ..- attr(*, "scaled:center")= num 37.1
## ..- attr(*, "scaled:scale")= num 8.49
## $ Edu : Factor w/ 7 levels "1","2","3","4",..: 4 3 3 1 5 7 5 5 4 4 ...
## $ Gender : Factor w/ 2 levels "0","1": 1 1 1 1 2 1 2 2 2 2 ...
## $ Income : num [1:334, 1] 0.259 0.259 -0.666 -0.666 1.184 ...
## ..- attr(*, "scaled:center")= num 2.72
## ..- attr(*, "scaled:scale")= num 1.08
## $ SpeededIP : num [1:334, 1] -1.844 -2.424 -1.974 -1.407 -0.603 ...
## ..- attr(*, "scaled:center")= num 0.00825
## ..- attr(*, "scaled:scale")= num 0.969
## $ StrategicIP : num [1:334, 1] -1.469 -1.599 0.751 -1.028 -0.522 ...
## ..- attr(*, "scaled:center")= num -0.0448
## ..- attr(*, "scaled:scale")= num 1.04
## $ Percep : num [1:334, 1] -0.0923 -0.6635 1.0988 -0.6263 0.2969 ...
## ..- attr(*, "scaled:center")= num 0.0484
## ..- attr(*, "scaled:scale")= num 1.02
## $ Caution : num [1:334, 1] 0.194 -0.592 -1.213 -0.057 1.131 ...
## ..- attr(*, "scaled:center")= num 0.0285
## ..- attr(*, "scaled:scale")= num 1.03
## $ Discount : num [1:334, 1] 0.492 1.608 1.569 1.196 -0.186 ...
## ..- attr(*, "scaled:center")= num -0.0189
## ..- attr(*, "scaled:scale")= num 1.03
## $ conservatism_factor: num [1:334, 1] 0.359 1.793 -0.143 0.491 0.307 ...
## ..- attr(*, "scaled:center")= num -0.0198
## ..- attr(*, "scaled:scale")= num 0.937
## $ dogmatism_factor : num [1:334, 1] 2.025 0.584 -1.742 1.494 -0.424 ...
## ..- attr(*, "scaled:center")= num -0.0148
## ..- attr(*, "scaled:scale")= num 0.927
## $ religiosity_factor : num [1:334, 1] -0.72 1.81 -0.582 0.061 1.824 ...
## ..- attr(*, "scaled:center")= num 0.00144
## ..- attr(*, "scaled:scale")= num 0.986Descriptive analyses
I recently learned about the modelsummary package, and it looks very very nice. So let’s try to create a descriptive table using the package.
df %>%
select(Age, Income, SpeededIP, StrategicIP, Percep, Caution, Discount, conservatism_factor, dogmatism_factor, religiosity_factor) %>%
datasummary_skim(.)| Unique | Missing Pct. | Mean | SD | Min | Median | Max | ||
|---|---|---|---|---|---|---|---|---|
| Age | 42 | 1 | −0.0 | 1.0 | −1.8 | −0.1 | 3.1 | |
| Income | 6 | 2 | −0.0 | 1.0 | −1.6 | 0.3 | 2.1 | |
| SpeededIP | 321 | 0 | 0.0 | 1.0 | −2.4 | −0.0 | 2.7 | |
| StrategicIP | 325 | 0 | −0.0 | 1.0 | −2.4 | −0.1 | 2.6 | |
| Percep | 315 | 0 | −0.0 | 1.0 | −2.6 | −0.0 | 3.1 | |
| Caution | 316 | 0 | −0.0 | 1.0 | −2.4 | 0.0 | 2.6 | |
| Discount | 319 | 0 | 0.0 | 1.0 | −2.4 | 0.1 | 1.9 | |
| conservatism_factor | 317 | 5 | −0.0 | 1.0 | −2.0 | −0.0 | 2.6 | |
| dogmatism_factor | 317 | 5 | 0.0 | 1.0 | −2.1 | −0.1 | 4.6 | |
| religiosity_factor | 317 | 5 | −0.0 | 1.0 | −0.9 | −0.6 | 2.8 |
Ok, this is actually super super cool! We can also summarize the categorical variables. By default datasummary_skim will only summarize numerical variables (above I selected a few specific variables of interest on purpose). We can also tell it to just summarize categorical variables:
datasummary_skim(df, type = "categorical")
## Warning: These variables were omitted because they include more than 50 levels: ID.| N | % | ||
|---|---|---|---|
| Edu | 1 | 2 | 0.6 |
| 2 | 56 | 16.8 | |
| 3 | 69 | 20.7 | |
| 4 | 48 | 14.4 | |
| 5 | 125 | 37.4 | |
| 6 | 24 | 7.2 | |
| 7 | 8 | 2.4 | |
| Gender | 0 | 163 | 48.8 |
| 1 | 167 | 50.0 |
Again, this is so so cool! I am a big fan of the package already! Let’s see if we can find some more use for it a bit later.
Ready for the Path Models?
Using lavaan for frequentist path models
Path models can also be seen as structural equation models with only observed and no latent variables. They do not model measurement error in a way as latent variable models do, but are thus usually also a bit less complex. We will first run a simple path model in lavaan, and then try to replicate the findings in brms and blavaan. Also, for now we will not control for demographic variables.
Let’s specify our lavaan model.
model1 <-
'
# direct associations (outcome ~ direct predictor)
conservatism_factor ~ c1*Age
dogmatism_factor ~ c2*Age
religiosity_factor ~ c3*Age
# indirect associations (indirect predictor ~ direct predictor)
SpeededIP ~ a1*Age
StrategicIP ~ a2*Age
Percep ~ a3*Age
Caution ~ a4*Age
Discount ~ a5*Age
# indirect associations (outcome ~ indirect predictor)
conservatism_factor ~ b11*SpeededIP + b12*StrategicIP + b13*Percep + b14*Caution + b15*Discount
dogmatism_factor ~ b21*SpeededIP + b22*StrategicIP + b23*Percep + b24*Caution + b25*Discount
religiosity_factor ~ b31*SpeededIP + b32*StrategicIP + b33*Percep + b34*Caution + b35*Discount
# residual covariances
# mediator residual covariance
SpeededIP ~~ StrategicIP + Percep + Caution + Discount
StrategicIP ~~ Percep + Caution + Discount
Percep ~~ Caution + Discount
Caution ~~ Discount
# dependent residual covariance
conservatism_factor ~~ dogmatism_factor + religiosity_factor
dogmatism_factor ~~ religiosity_factor
# effect decomposition
# x = age, m1 = SpeededIP, m2 = StrategicIP, m3 = Percep, m4 = Caution, m5 = Discount
# y1 = conservatism_factor;
ind_x_m1_y1 := a1*b11
ind_x_m2_y1 := a2*b12
ind_x_m3_y1 := a3*b13
ind_x_m4_y1 := a4*b14
ind_x_m5_y1 := a4*b15
ind_x_y1 := ind_x_m1_y1 + ind_x_m2_y1 + ind_x_m3_y1 + ind_x_m4_y1 + ind_x_m5_y1
tot_x_y1 := ind_x_y1 + c1
# y2 = dogmatism_factor
ind_x_m1_y2 := a1*b21
ind_x_m2_y2 := a2*b22
ind_x_m3_y2 := a3*b23
ind_x_m4_y2 := a4*b24
ind_x_m5_y2 := a4*b25
ind_x_y2 := ind_x_m1_y2 + ind_x_m2_y2 + ind_x_m3_y2 + ind_x_m4_y2 + ind_x_m5_y2
tot_x_y2 := ind_x_y2 + c2
# y3 = religiosity_factor
ind_x_m1_y3 := a1*b31
ind_x_m2_y3 := a2*b32
ind_x_m3_y3 := a3*b33
ind_x_m4_y3 := a4*b34
ind_x_m5_y3 := a4*b35
ind_x_y3 := ind_x_m1_y3 + ind_x_m2_y3 + ind_x_m3_y3 + ind_x_m4_y3 + ind_x_m5_y3
tot_x_y3 := ind_x_y3 + c3
'Let’s fit and run this path model in lavaan.
fit1 <- sem(model = model1, data = df,
std.ov = TRUE, mimic = "Mplus",
estimator = "ML", missing = "fiml")
## Warning in lav_data_full(data = data, group = group, cluster = cluster, : lavaan WARNING: 4 cases were deleted due to missing values in
## exogenous variable(s), while fixed.x = TRUE.parest1 <- parameterestimates(fit1)
parest1 %>%
slice(55:75) %>%
select(label, est, se, pvalue, ci.lower, ci.upper) %>%
tt(.)| label | est | se | pvalue | ci.lower | ci.upper |
|---|---|---|---|---|---|
| ind_x_m1_y1 | -0.01473 | 0.0169 | 3.85e-01 | -0.04794 | 0.01847 |
| ind_x_m2_y1 | 0.02956 | 0.0137 | 3.15e-02 | 0.00262 | 0.05650 |
| ind_x_m3_y1 | 0.00765 | 0.0119 | 5.21e-01 | -0.01572 | 0.03103 |
| ind_x_m4_y1 | 0.04944 | 0.0178 | 5.36e-03 | 0.01464 | 0.08423 |
| ind_x_m5_y1 | 0.03731 | 0.0166 | 2.46e-02 | 0.00478 | 0.06985 |
| ind_x_y1 | 0.10923 | 0.0369 | 3.10e-03 | 0.03685 | 0.18161 |
| tot_x_y1 | 0.27269 | 0.0585 | 3.16e-06 | 0.15801 | 0.38738 |
| ind_x_m1_y2 | -0.04499 | 0.0197 | 2.21e-02 | -0.08351 | -0.00647 |
| ind_x_m2_y2 | 0.01060 | 0.0100 | 2.91e-01 | -0.00907 | 0.03027 |
| ind_x_m3_y2 | -0.00642 | 0.0126 | 6.12e-01 | -0.03122 | 0.01837 |
| ind_x_m4_y2 | 0.01187 | 0.0161 | 4.61e-01 | -0.01970 | 0.04343 |
| ind_x_m5_y2 | 0.01522 | 0.0162 | 3.48e-01 | -0.01654 | 0.04698 |
| ind_x_y2 | -0.01373 | 0.0347 | 6.92e-01 | -0.08168 | 0.05423 |
| tot_x_y2 | -0.07333 | 0.0595 | 2.18e-01 | -0.19000 | 0.04333 |
| ind_x_m1_y3 | -0.00921 | 0.0177 | 6.03e-01 | -0.04395 | 0.02552 |
| ind_x_m2_y3 | 0.02789 | 0.0136 | 4.04e-02 | 0.00122 | 0.05456 |
| ind_x_m3_y3 | 0.01286 | 0.0128 | 3.16e-01 | -0.01228 | 0.03800 |
| ind_x_m4_y3 | 0.02405 | 0.0165 | 1.44e-01 | -0.00819 | 0.05630 |
| ind_x_m5_y3 | 0.00584 | 0.0158 | 7.12e-01 | -0.02515 | 0.03683 |
| ind_x_y3 | 0.06143 | 0.0350 | 7.94e-02 | -0.00720 | 0.13006 |
| tot_x_y3 | 0.13190 | 0.0593 | 2.62e-02 | 0.01562 | 0.24817 |
Let’s remind ourselves of what is what:
x = age, m1 = SpeededIP, m2 = StrategicIP, m3 = Percep, m4 = Caution, m5 = Discount, y1 = conservatism, y2 = dogmatism, y3 = religiosity.
The model indicates
Significant, but very small positive indirect associations between Age and Conservatism via the cognitive factors StrategicIP, Caution, and Discount.
No significant indirect associations between Age, Dogmatism and Religiosity via cognitive factors.
We can also check a bit more information using summary.
summary(fit1)
## lavaan 0.6.17 ended normally after 18 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 52
##
## Used Total
## Number of observations 330 334
## Number of missing patterns 2
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## conservatism_factor ~
## Age (c1) 0.163 0.059 2.786 0.005
## dogmatism_factor ~
## Age (c2) -0.060 0.063 -0.953 0.341
## religiosity_factor ~
## Age (c3) 0.070 0.062 1.138 0.255
## SpeededIP ~
## Age (a1) 0.287 0.053 5.436 0.000
## StrategicIP ~
## Age (a2) -0.151 0.054 -2.779 0.005
## Percep ~
## Age (a3) 0.216 0.054 4.024 0.000
## Caution ~
## Age (a4) 0.274 0.053 5.169 0.000
## Discount ~
## Age (a5) -0.072 0.055 -1.305 0.192
## conservatism_factor ~
## SpeeddIP (b11) -0.051 0.058 -0.881 0.378
## StrtgcIP (b12) -0.195 0.058 -3.394 0.001
## Percep (b13) 0.035 0.054 0.650 0.516
## Caution (b14) 0.181 0.055 3.305 0.001
## Discount (b15) 0.136 0.055 2.496 0.013
## dogmatism_factor ~
## SpeeddIP (b21) -0.157 0.062 -2.524 0.012
## StrtgcIP (b22) -0.070 0.061 -1.142 0.253
## Percep (b23) -0.030 0.058 -0.512 0.609
## Caution (b24) 0.043 0.058 0.744 0.457
## Discount (b25) 0.056 0.058 0.955 0.339
## religiosity_factor ~
## SpeeddIP (b31) -0.032 0.062 -0.522 0.602
## StrtgcIP (b32) -0.184 0.061 -3.035 0.002
## Percep (b33) 0.059 0.057 1.035 0.300
## Caution (b34) 0.088 0.058 1.525 0.127
## Discount (b35) 0.021 0.058 0.370 0.711
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## .SpeededIP ~~
## .StrategicIP 0.340 0.055 6.163 0.000
## .Percep -0.177 0.052 -3.393 0.001
## .Caution 0.109 0.051 2.146 0.032
## .Discount -0.167 0.053 -3.129 0.002
## .StrategicIP ~~
## .Percep -0.089 0.053 -1.668 0.095
## .Caution 0.051 0.052 0.985 0.325
## .Discount -0.292 0.056 -5.180 0.000
## .Percep ~~
## .Caution 0.080 0.052 1.542 0.123
## .Discount 0.035 0.053 0.652 0.514
## .Caution ~~
## .Discount -0.032 0.053 -0.608 0.543
## .conservatism_factor ~~
## .dogmatism_fctr 0.235 0.052 4.540 0.000
## .religisty_fctr 0.270 0.052 5.198 0.000
## .dogmatism_factor ~~
## .religisty_fctr 0.074 0.053 1.391 0.164
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|)
## .consrvtsm_fctr -0.006 0.052 -0.111 0.912
## .dogmatism_fctr -0.000 0.055 -0.001 0.999
## .religisty_fctr -0.002 0.054 -0.038 0.970
## .SpeededIP 0.000 0.053 0.000 1.000
## .StrategicIP -0.000 0.054 -0.000 1.000
## .Percep 0.000 0.054 0.000 1.000
## .Caution 0.000 0.053 0.000 1.000
## .Discount -0.000 0.055 -0.000 1.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .consrvtsm_fctr 0.833 0.067 12.530 0.000
## .dogmatism_fctr 0.947 0.076 12.530 0.000
## .religisty_fctr 0.928 0.074 12.530 0.000
## .SpeededIP 0.915 0.071 12.845 0.000
## .StrategicIP 0.974 0.076 12.845 0.000
## .Percep 0.950 0.074 12.845 0.000
## .Caution 0.922 0.072 12.845 0.000
## .Discount 0.992 0.077 12.845 0.000
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|)
## ind_x_m1_y1 -0.015 0.017 -0.870 0.385
## ind_x_m2_y1 0.030 0.014 2.150 0.032
## ind_x_m3_y1 0.008 0.012 0.642 0.521
## ind_x_m4_y1 0.049 0.018 2.785 0.005
## ind_x_m5_y1 0.037 0.017 2.248 0.025
## ind_x_y1 0.109 0.037 2.958 0.003
## tot_x_y1 0.273 0.059 4.660 0.000
## ind_x_m1_y2 -0.045 0.020 -2.289 0.022
## ind_x_m2_y2 0.011 0.010 1.057 0.291
## ind_x_m3_y2 -0.006 0.013 -0.508 0.612
## ind_x_m4_y2 0.012 0.016 0.737 0.461
## ind_x_m5_y2 0.015 0.016 0.939 0.348
## ind_x_y2 -0.014 0.035 -0.396 0.692
## tot_x_y2 -0.073 0.060 -1.232 0.218
## ind_x_m1_y3 -0.009 0.018 -0.520 0.603
## ind_x_m2_y3 0.028 0.014 2.050 0.040
## ind_x_m3_y3 0.013 0.013 1.003 0.316
## ind_x_m4_y3 0.024 0.016 1.462 0.144
## ind_x_m5_y3 0.006 0.016 0.369 0.712
## ind_x_y3 0.061 0.035 1.754 0.079
## tot_x_y3 0.132 0.059 2.223 0.026Also note: Our model perfectly fits the data, but this is because the model is just-identified (i.e., there are 0 degrees of freedom): Every pair of variable is connected - our model is just a different way to describe the full data. This does not mean that the model syntax is incorrect - this is often the case in models with just observed variables and indirect associations.
There is much more in this output that can be of interest. But as we want to focus on the Bayesian analyses here, let’s move on for now.
blavaan: The Bayesian counterpart of lavaan
Using very similar syntax as lavaan, blavaan allows us to fit Bayesian structural equation models (or path models as here).
Let’s check the default priors of blavaan.
dpriors()
## nu alpha lambda beta theta
## "normal(0,32)" "normal(0,10)" "normal(0,10)" "normal(0,10)" "gamma(1,.5)[sd]"
## psi rho ibpsi tau
## "gamma(1,.5)[sd]" "beta(1,1)" "wishart(3,iden)" "normal(0,1.5)"Contrary to brms, in many cases we can use these default priors. We can modify the prior distributions, and in particular, we want to specify the prior of the regressions \(\beta\) as \(N(0,1)\), as we will later use this prior in brms. It should not matter to much, especilly for the other parameters, but this is also good to try out on how to specify priors.
mydp <- dpriors(beta="normal(0,1)")
mydp
## nu alpha lambda beta theta
## "normal(0,32)" "normal(0,10)" "normal(0,10)" "normal(0,1)" "gamma(1,.5)[sd]"
## psi rho ibpsi tau
## "gamma(1,.5)[sd]" "beta(1,1)" "wishart(3,iden)" "normal(0,1.5)"Let’s try to fit the same model as previously.
fit2 <- bsem(model1, dp = mydp, data = df)
## Warning in lav_data_full(data = data, group = group, cluster = cluster, : lavaan WARNING: 4 cases were deleted due to missing values in
## exogenous variable(s), while fixed.x = TRUE.
## Computing post-estimation metrics (including lvs if requested)...
## Warning in tmpmat[lower.tri(tmpmat)] <- lavpartable$est[covpars]: number of items to
## replace is not a multiple of replacement lengthsummary(fit2)
## blavaan 0.5.3 ended normally after 1000 iterations
##
## Estimator BAYES
## Optimization method MCMC
## Number of model parameters 52
##
## Used Total
## Number of observations 330 334
## Number of missing patterns 2
##
## Statistic MargLogLik PPP
## Value -3662.530 0.518
##
## Parameter Estimates:
##
##
## Regressions:
## Estimate Post.SD pi.lower pi.upper Rhat Prior
## conservatism_factor ~
## Age (c1) 0.162 0.059 0.048 0.279 1.000 normal(0,1)
## dogmatism_factor ~
## Age (c2) -0.061 0.065 -0.189 0.065 0.999 normal(0,1)
## religiosity_factor ~
## Age (c3) 0.069 0.063 -0.054 0.191 1.000 normal(0,1)
## SpeededIP ~
## Age (a1) 0.284 0.051 0.183 0.389 0.999 normal(0,1)
## StrategicIP ~
## Age (a2) -0.151 0.055 -0.261 -0.045 1.000 normal(0,1)
## Percep ~
## Age (a3) 0.211 0.052 0.108 0.315 1.000 normal(0,1)
## Caution ~
## Age (a4) 0.272 0.054 0.168 0.375 1.000 normal(0,1)
## Discount ~
## Age (a5) -0.071 0.056 -0.177 0.039 0.999 normal(0,1)
## conservatism_factor ~
## SpeeddIP (b11) -0.052 0.059 -0.167 0.065 1.000 normal(0,1)
## StrtgcIP (b12) -0.197 0.059 -0.313 -0.083 1.001 normal(0,1)
## Percep (b13) 0.036 0.056 -0.073 0.149 1.000 normal(0,1)
## Caution (b14) 0.182 0.056 0.075 0.289 1.000 normal(0,1)
## Discount (b15) 0.135 0.054 0.030 0.246 1.000 normal(0,1)
## dogmatism_factor ~
## SpeeddIP (b21) -0.158 0.066 -0.289 -0.028 1.000 normal(0,1)
## StrtgcIP (b22) -0.072 0.063 -0.193 0.050 1.000 normal(0,1)
## Percep (b23) -0.030 0.060 -0.148 0.089 1.000 normal(0,1)
## Caution (b24) 0.046 0.057 -0.067 0.153 1.000 normal(0,1)
## Discount (b25) 0.054 0.058 -0.057 0.167 1.000 normal(0,1)
## religiosity_factor ~
## SpeeddIP (b31) -0.031 0.062 -0.151 0.094 1.000 normal(0,1)
## StrtgcIP (b32) -0.184 0.063 -0.308 -0.061 0.999 normal(0,1)
## Percep (b33) 0.059 0.060 -0.061 0.176 1.000 normal(0,1)
## Caution (b34) 0.087 0.059 -0.025 0.204 0.999 normal(0,1)
## Discount (b35) 0.021 0.058 -0.095 0.131 1.000 normal(0,1)
##
## Covariances:
## Estimate Post.SD pi.lower pi.upper Rhat Prior
## .SpeededIP ~~
## .StrategicIP 0.341 0.056 0.238 0.456 1.000 lkj_corr(1)
## .Percep -0.176 0.051 -0.279 -0.081 0.999 lkj_corr(1)
## .Caution 0.110 0.052 0.011 0.215 0.999 lkj_corr(1)
## .Discount -0.167 0.054 -0.276 -0.063 1.000 lkj_corr(1)
## .StrategicIP ~~
## .Percep -0.087 0.054 -0.196 0.017 1.000 lkj_corr(1)
## .Caution 0.050 0.054 -0.056 0.156 1.000 lkj_corr(1)
## .Discount -0.296 0.057 -0.416 -0.187 1.000 lkj_corr(1)
## .Percep ~~
## .Caution 0.079 0.053 -0.023 0.182 0.999 lkj_corr(1)
## .Discount 0.035 0.055 -0.074 0.143 0.999 lkj_corr(1)
## .Caution ~~
## .Discount -0.033 0.053 -0.135 0.069 0.999 lkj_corr(1)
## .conservatism_factor ~~
## .dogmatism_fctr 0.242 0.056 0.132 0.349 1.001 lkj_corr(1)
## .religisty_fctr 0.275 0.054 0.177 0.385 1.002 lkj_corr(1)
## .dogmatism_factor ~~
## .religisty_fctr 0.075 0.056 -0.035 0.187 1.001 lkj_corr(1)
##
## Intercepts:
## Estimate Post.SD pi.lower pi.upper Rhat Prior
## .consrvtsm_fctr -0.003 0.054 -0.107 0.103 0.999 normal(0,10)
## .dogmatism_fctr 0.004 0.057 -0.105 0.116 0.999 normal(0,10)
## .religisty_fctr -0.011 0.056 -0.122 0.094 1.000 normal(0,10)
## .SpeededIP 0.003 0.053 -0.099 0.109 1.000 normal(0,10)
## .StrategicIP 0.013 0.055 -0.096 0.119 1.000 normal(0,10)
## .Percep -0.023 0.052 -0.124 0.083 1.000 normal(0,10)
## .Caution 0.001 0.053 -0.103 0.103 1.000 normal(0,10)
## .Discount -0.000 0.055 -0.111 0.109 1.000 normal(0,10)
##
## Variances:
## Estimate Post.SD pi.lower pi.upper Rhat Prior
## .consrvtsm_fctr 0.867 0.071 0.737 1.020 1.000 gamma(1,.5)[sd]
## .dogmatism_fctr 0.987 0.080 0.841 1.155 1.001 gamma(1,.5)[sd]
## .religisty_fctr 0.943 0.074 0.807 1.098 1.000 gamma(1,.5)[sd]
## .SpeededIP 0.926 0.073 0.798 1.080 0.999 gamma(1,.5)[sd]
## .StrategicIP 0.995 0.080 0.853 1.163 1.000 gamma(1,.5)[sd]
## .Percep 0.941 0.075 0.802 1.097 0.999 gamma(1,.5)[sd]
## .Caution 0.953 0.074 0.821 1.114 0.999 gamma(1,.5)[sd]
## .Discount 1.026 0.082 0.875 1.198 1.000 gamma(1,.5)[sd]
##
## Defined Parameters:
## Estimate Post.SD pi.lower pi.upper Rhat Prior
## ind_x_m1_y1 -0.015 0.017 -0.049 0.019
## ind_x_m2_y1 0.030 0.014 0.001 0.058
## ind_x_m3_y1 0.008 0.012 -0.017 0.032
## ind_x_m4_y1 0.049 0.018 0.014 0.085
## ind_x_m5_y1 0.037 0.017 0.004 0.070
## ind_x_y1 0.109 0.038 0.034 0.183
## tot_x_y1 0.271 0.059 0.155 0.387
## ind_x_m1_y2 -0.045 0.021 -0.085 -0.004
## ind_x_m2_y2 0.011 0.011 -0.011 0.032
## ind_x_m3_y2 -0.006 0.013 -0.032 0.019
## ind_x_m4_y2 0.012 0.016 -0.019 0.044
## ind_x_m5_y2 0.015 0.017 -0.018 0.047
## ind_x_y2 -0.013 0.036 -0.084 0.058
## tot_x_y2 -0.074 0.062 -0.195 0.047
## ind_x_m1_y3 -0.009 0.018 -0.043 0.026
## ind_x_m2_y3 0.028 0.014 -0.001 0.056
## ind_x_m3_y3 0.012 0.013 -0.014 0.038
## ind_x_m4_y3 0.024 0.017 -0.010 0.058
## ind_x_m5_y3 0.006 0.016 -0.026 0.037
## ind_x_y3 0.061 0.036 -0.009 0.131
## tot_x_y3 0.130 0.060 0.012 0.248The results are very similar.
Path models in brms
We can do something very similar also using brms! brms is a very flexible and powerful packages that is very flexible with modeling approaches. I think blavaan is the more straight-forward choice as soon as you want to specifiy latent variable models, but for models with just observed variables brms may even be the better pick, even in the case of path models. Let’s try it out and compare the results.
Selecting priors
I already scaled all variables of interest. This will make it a bit simpler to select priors. I will use weakly informative priors for slope and intercept parameters again (see https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations). I do this merely to keep things simple and straightforward here - usually a few more thoughts should be placed in specifying your prior. The selected prior is \(N(0,1)\). Let’s plot to prior to see what kind of an effect it may have on our parameters.
colors <- ipsum_pal()(1)
prior(normal(0, 1)) %>%
parse_dist() %>%
ggplot(aes(xdist = .dist_obj, y = prior)) +
stat_halfeye(.width = 0.683, p_limits = c(.001, .999),
slab_fill = colors[1], slab_alpha = 0.8) +
scale_y_discrete(NULL, breaks = NULL, expand = expansion(add = 0.1)) +
labs(title = "Prior N(0,1)", x = "Parameter Space", y = "") +
coord_cartesian(xlim = c(-4, 4)) +
scale_x_continuous(breaks = seq(-4, 4, by = 1)) +
theme_ipsum() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank()) 
Setting up the model
In brms we have to specify our model a bit different than previously. However, in principle, we are doing nothing else than fitting a multivariate linear regression once again. The way we set up this model should thus remind both of what we specified in the previous blog post, and of the lavaan syntax from above.
model2 <-
# direct and indirect paths to ideological factors
bf(conservatism_factor ~ Age + SpeededIP + StrategicIP + Percep + Caution + Discount) +
bf(dogmatism_factor ~ Age + SpeededIP + StrategicIP + Percep + Caution + Discount) +
bf(religiosity_factor ~ Age + SpeededIP + StrategicIP + Percep + Caution + Discount) +
# paths from age to indirect predictors (cognitive factors)
bf(SpeededIP ~ Age) +
bf(StrategicIP ~ Age) +
bf(Percep ~ Age) +
bf(Caution ~ Age) +
bf(Discount ~ Age)Let’s check how our priors are set by default in brms.
get_prior(model2 + set_rescor(FALSE), data = df)
## Warning: Rows containing NAs were excluded from the model.
## prior class coef group resp dpar nlpar lb ub
## (flat) b
## (flat) Intercept
## (flat) b Caution
## (flat) b Age Caution
## student_t(3, 0, 2.5) Intercept Caution
## student_t(3, 0, 2.5) sigma Caution 0
## (flat) b conservatismfactor
## (flat) b Age conservatismfactor
## (flat) b Caution conservatismfactor
## (flat) b Discount conservatismfactor
## (flat) b Percep conservatismfactor
## (flat) b SpeededIP conservatismfactor
## (flat) b StrategicIP conservatismfactor
## student_t(3, 0, 2.5) Intercept conservatismfactor
## student_t(3, 0, 2.5) sigma conservatismfactor 0
## (flat) b Discount
## (flat) b Age Discount
## student_t(3, 0.1, 2.5) Intercept Discount
## student_t(3, 0, 2.5) sigma Discount 0
## (flat) b dogmatismfactor
## (flat) b Age dogmatismfactor
## (flat) b Caution dogmatismfactor
## (flat) b Discount dogmatismfactor
## (flat) b Percep dogmatismfactor
## (flat) b SpeededIP dogmatismfactor
## (flat) b StrategicIP dogmatismfactor
## student_t(3, 0, 2.5) Intercept dogmatismfactor
## student_t(3, 0, 2.5) sigma dogmatismfactor 0
## (flat) b Percep
## (flat) b Age Percep
## student_t(3, -0.1, 2.5) Intercept Percep
## student_t(3, 0, 2.5) sigma Percep 0
## (flat) b religiosityfactor
## (flat) b Age religiosityfactor
## (flat) b Caution religiosityfactor
## (flat) b Discount religiosityfactor
## (flat) b Percep religiosityfactor
## (flat) b SpeededIP religiosityfactor
## (flat) b StrategicIP religiosityfactor
## student_t(3, -0.6, 2.5) Intercept religiosityfactor
## student_t(3, 0, 2.5) sigma religiosityfactor 0
## (flat) b SpeededIP
## (flat) b Age SpeededIP
## student_t(3, 0, 2.5) Intercept SpeededIP
## student_t(3, 0, 2.5) sigma SpeededIP 0
## (flat) b StrategicIP
## (flat) b Age StrategicIP
## student_t(3, 0, 2.5) Intercept StrategicIP
## student_t(3, 0, 2.5) sigma StrategicIP 0
## source
## default
## default
## default
## (vectorized)
## default
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## default
## default
## (vectorized)
## default
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## default
## default
## (vectorized)
## default
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## default
## default
## (vectorized)
## default
## default
## default
## (vectorized)
## default
## defaultBy default brms sets non-informative flat priors for our slope parameters. We thus have to change this, and can do so directly when fitting the model.
Fitting and summarising the first model
m3 <- brm(model2 + set_rescor(FALSE),
data = df,
prior = c(prior(normal(0, 1), class = "Intercept"),
prior(normal(0, 1), class = "b")),
seed = bayes_seed,
file = "02_fits/fit_m3")We get a warning message:
Warning: Rows containing NAs were excluded from the model.We are not happy with this, but for now we will deal with it later.
Check the results. As usual, we can also use summary(m1.1) or print(m1.1) to do so.
summary(m3, prob = 0.95) # prob = 0.95 is the default
## Family: MV(gaussian, gaussian, gaussian, gaussian, gaussian, gaussian, gaussian, gaussian)
## Links: mu = identity; sigma = identity
## mu = identity; sigma = identity
## mu = identity; sigma = identity
## mu = identity; sigma = identity
## mu = identity; sigma = identity
## mu = identity; sigma = identity
## mu = identity; sigma = identity
## mu = identity; sigma = identity
## Formula: conservatism_factor ~ Age + SpeededIP + StrategicIP + Percep + Caution + Discount
## dogmatism_factor ~ Age + SpeededIP + StrategicIP + Percep + Caution + Discount
## religiosity_factor ~ Age + SpeededIP + StrategicIP + Percep + Caution + Discount
## SpeededIP ~ Age
## StrategicIP ~ Age
## Percep ~ Age
## Caution ~ Age
## Discount ~ Age
## Data: df (Number of observations: 314)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## conservatismfactor_Intercept -0.00 0.05 -0.11 0.11 1.00 6337
## dogmatismfactor_Intercept 0.00 0.06 -0.11 0.12 1.00 6525
## religiosityfactor_Intercept -0.01 0.05 -0.11 0.09 1.00 6099
## SpeededIP_Intercept 0.01 0.05 -0.09 0.12 1.00 6149
## StrategicIP_Intercept 0.01 0.06 -0.10 0.13 1.00 5718
## Percep_Intercept -0.03 0.05 -0.13 0.08 1.00 6003
## Caution_Intercept 0.03 0.05 -0.08 0.14 1.00 5574
## Discount_Intercept 0.01 0.06 -0.10 0.12 1.00 5602
## conservatismfactor_Age 0.16 0.06 0.05 0.29 1.00 4368
## conservatismfactor_SpeededIP -0.05 0.06 -0.17 0.06 1.00 4502
## conservatismfactor_StrategicIP -0.20 0.06 -0.32 -0.08 1.00 5018
## conservatismfactor_Percep 0.04 0.06 -0.08 0.15 1.00 5278
## conservatismfactor_Caution 0.18 0.06 0.07 0.29 1.00 5297
## conservatismfactor_Discount 0.14 0.06 0.02 0.25 1.00 5710
## dogmatismfactor_Age -0.06 0.06 -0.18 0.06 1.00 3949
## dogmatismfactor_SpeededIP -0.16 0.06 -0.28 -0.03 1.00 4213
## dogmatismfactor_StrategicIP -0.07 0.06 -0.19 0.05 1.00 4326
## dogmatismfactor_Percep -0.03 0.06 -0.15 0.09 1.00 4723
## dogmatismfactor_Caution 0.04 0.06 -0.07 0.16 1.00 5841
## dogmatismfactor_Discount 0.05 0.06 -0.07 0.18 1.00 5344
## religiosityfactor_Age 0.07 0.06 -0.05 0.19 1.00 4160
## religiosityfactor_SpeededIP -0.03 0.06 -0.15 0.09 1.00 3946
## religiosityfactor_StrategicIP -0.18 0.06 -0.31 -0.06 1.00 4377
## religiosityfactor_Percep 0.06 0.06 -0.06 0.18 1.00 4796
## religiosityfactor_Caution 0.08 0.06 -0.03 0.20 1.00 5670
## religiosityfactor_Discount 0.02 0.06 -0.10 0.14 1.00 4846
## SpeededIP_Age 0.28 0.05 0.18 0.39 1.00 5559
## StrategicIP_Age -0.14 0.06 -0.26 -0.03 1.00 7492
## Percep_Age 0.21 0.05 0.11 0.32 1.00 5614
## Caution_Age 0.30 0.06 0.18 0.41 1.00 6587
## Discount_Age -0.08 0.06 -0.19 0.03 1.00 6507
## Tail_ESS
## conservatismfactor_Intercept 2623
## dogmatismfactor_Intercept 2985
## religiosityfactor_Intercept 2855
## SpeededIP_Intercept 3067
## StrategicIP_Intercept 2994
## Percep_Intercept 3038
## Caution_Intercept 3158
## Discount_Intercept 3195
## conservatismfactor_Age 3269
## conservatismfactor_SpeededIP 3487
## conservatismfactor_StrategicIP 3386
## conservatismfactor_Percep 3219
## conservatismfactor_Caution 3159
## conservatismfactor_Discount 3473
## dogmatismfactor_Age 3397
## dogmatismfactor_SpeededIP 3290
## dogmatismfactor_StrategicIP 3359
## dogmatismfactor_Percep 3091
## dogmatismfactor_Caution 3100
## dogmatismfactor_Discount 3030
## religiosityfactor_Age 2947
## religiosityfactor_SpeededIP 3058
## religiosityfactor_StrategicIP 3375
## religiosityfactor_Percep 3125
## religiosityfactor_Caution 2985
## religiosityfactor_Discount 3287
## SpeededIP_Age 3144
## StrategicIP_Age 2951
## Percep_Age 3037
## Caution_Age 2954
## Discount_Age 3069
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma_conservatismfactor 0.93 0.04 0.86 1.00 1.00 4946 3050
## sigma_dogmatismfactor 0.99 0.04 0.91 1.07 1.00 5632 2831
## sigma_religiosityfactor 0.97 0.04 0.90 1.05 1.00 6108 2799
## sigma_SpeededIP 0.97 0.04 0.89 1.04 1.00 7028 3235
## sigma_StrategicIP 0.99 0.04 0.92 1.08 1.00 6009 3152
## sigma_Percep 0.96 0.04 0.89 1.04 1.00 5479 3026
## sigma_Caution 0.97 0.04 0.89 1.05 1.00 5749 3318
## sigma_Discount 1.00 0.04 0.93 1.09 1.00 6932 2907
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).This time, we need some trickery to compute the indirect effects. Let’s remind ourselves how we did this previously - what I write out below, is what we specified in the lavaan syntax.
The indirect effect in general is calculated by the multiplication of the path from the direct predictor to the indirect predictor (a)
- In this case, this is the path from Age to SpeededIP, StrategicIP, Percep, Caution, Discount respectively
with the path from the indirect predictor to the outcome (b)
- In this case, this is the path from SpeededIP, StrategicIP, Percep, Caution, Discount to conservatism_factor, dogmatism_factor, or religiosity_factor
And the direct paths (c) are the paths from the direct predictor to the outcome
- In this case, this is the path from Age to conservatism_factor, dogmatism_factor, or religiosity_factor
The overall indirect effect for one outcome is nothing but the sum of all indirect effects to that outcome.
The total effects are nothing but the sum of the indirect and direct effects.
In brms we get the posterior of all of those effects. We can use the posterior draws to calculate these effects - though some manual work is needed. Let’s first extract the draws.
m3_draws <- as_draws_df(m3)
head(m3_draws)
## # A draws_df: 6 iterations, 1 chains, and 41 variables
## b_conservatismfactor_Intercept b_dogmatismfactor_Intercept
## 1 0.034 0.0476
## 2 -0.073 0.0472
## 3 0.067 -0.0393
## 4 -0.057 0.0448
## 5 0.011 -0.0044
## 6 -0.015 0.0076
## b_religiosityfactor_Intercept b_SpeededIP_Intercept b_StrategicIP_Intercept
## 1 0.089 0.0604 0.0062
## 2 -0.102 0.0793 -0.0450
## 3 0.061 -0.0360 0.0651
## 4 -0.091 0.0584 -0.0283
## 5 0.042 0.0086 0.0081
## 6 -0.064 0.0171 0.0299
## b_Percep_Intercept b_Caution_Intercept b_Discount_Intercept
## 1 -0.041 0.0104 -0.032
## 2 -0.090 0.0472 -0.047
## 3 0.026 0.0121 0.065
## 4 -0.091 0.0820 -0.047
## 5 -0.029 0.0161 -0.045
## 6 -0.024 -0.0012 0.057
## # ... with 33 more variables
## # ... hidden reserved variables {'.chain', '.iteration', '.draw'}Now let’s try to calculate the indirect effects.
m3_draws <- m3_draws %>%
mutate(
# Age -> Conservatism
ind_Age_Spe_Con = b_SpeededIP_Age * b_conservatismfactor_SpeededIP,
ind_Age_Str_Con = b_StrategicIP_Age * b_conservatismfactor_StrategicIP,
ind_Age_Per_Con = b_Percep_Age * b_conservatismfactor_Percep,
ind_Age_Cau_Con = b_Caution_Age * b_conservatismfactor_Caution,
ind_Age_Dis_Con = b_Discount_Age * b_conservatismfactor_Discount,
ind_Age_Con = ind_Age_Spe_Con + ind_Age_Str_Con + ind_Age_Per_Con + ind_Age_Cau_Con + ind_Age_Dis_Con,
tot_Age_Con = ind_Age_Con + b_conservatismfactor_Age,
# Age -> Dogmatism
ind_Age_Spe_Dog = b_SpeededIP_Age * b_dogmatismfactor_SpeededIP,
ind_Age_Str_Dog = b_StrategicIP_Age * b_dogmatismfactor_StrategicIP,
ind_Age_Per_Dog = b_Percep_Age * b_dogmatismfactor_Percep,
ind_Age_Cau_Dog = b_Caution_Age * b_dogmatismfactor_Caution,
ind_Age_Dis_Dog = b_Discount_Age * b_dogmatismfactor_Discount,
ind_Age_Dog = ind_Age_Spe_Dog + ind_Age_Str_Dog + ind_Age_Per_Dog + ind_Age_Cau_Dog + ind_Age_Dis_Dog,
tot_Age_Dog = ind_Age_Dog + b_dogmatismfactor_Age,
# Age -> Religiosity
ind_Age_Spe_Rel = b_SpeededIP_Age * b_religiosityfactor_SpeededIP,
ind_Age_Str_Rel = b_StrategicIP_Age * b_religiosityfactor_StrategicIP,
ind_Age_Per_Rel = b_Percep_Age * b_religiosityfactor_Percep,
ind_Age_Cau_Rel = b_Caution_Age * b_religiosityfactor_Caution,
ind_Age_Dis_Rel = b_Discount_Age * b_religiosityfactor_Discount,
ind_Age_Rel = ind_Age_Spe_Rel + ind_Age_Str_Rel + ind_Age_Per_Rel + ind_Age_Cau_Rel + ind_Age_Dis_Rel,
tot_Age_Rel = ind_Age_Rel + b_religiosityfactor_Age
)Use the describe_posterior function to get some summary statistics of what we calculated.
m3_draws %>%
select(45:65) %>%
describe_posterior(.,
centrality = "median",
dispersion = FALSE,
ci = 0.95,
ci_method = "hdi",
rope_range = c(-0.02, 0.02),
rope_ci = 1)
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Summary of Posterior Distribution
##
## Parameter | Median | 95% CI | pd | ROPE | % in ROPE
## ---------------------------------------------------------------------------------
## ind_Age_Spe_Con | -0.01 | [-0.05, 0.02] | 80.58% | [-0.02, 0.02] | 61.68%
## ind_Age_Str_Con | 0.03 | [ 0.00, 0.06] | 99.48% | [-0.02, 0.02] | 30.18%
## ind_Age_Per_Con | 7.16e-03 | [-0.02, 0.03] | 73.52% | [-0.02, 0.02] | 83.62%
## ind_Age_Cau_Con | 0.05 | [ 0.02, 0.09] | 99.90% | [-0.02, 0.02] | 3.15%
## ind_Age_Dis_Con | -9.59e-03 | [-0.03, 0.00] | 91.27% | [-0.02, 0.02] | 83.33%
## ind_Age_Con | 0.06 | [ 0.00, 0.13] | 97.20% | [-0.02, 0.02] | 8.88%
## tot_Age_Con | 0.23 | [ 0.12, 0.34] | 100% | [-0.02, 0.02] | 0%
## ind_Age_Spe_Dog | -0.04 | [-0.09, -0.01] | 99.25% | [-0.02, 0.02] | 10.08%
## ind_Age_Str_Dog | 9.17e-03 | [-0.01, 0.03] | 87.12% | [-0.02, 0.02] | 83.33%
## ind_Age_Per_Dog | -5.80e-03 | [-0.04, 0.02] | 69.17% | [-0.02, 0.02] | 82.83%
## ind_Age_Cau_Dog | 0.01 | [-0.02, 0.05] | 77.25% | [-0.02, 0.02] | 63.70%
## ind_Age_Dis_Dog | -3.16e-03 | [-0.02, 0.01] | 76.20% | [-0.02, 0.02] | 96.75%
## ind_Age_Dog | -0.03 | [-0.10, 0.03] | 83.00% | [-0.02, 0.02] | 30.88%
## tot_Age_Dog | -0.09 | [-0.20, 0.01] | 95.38% | [-0.02, 0.02] | 7.45%
## ind_Age_Spe_Rel | -8.77e-03 | [-0.04, 0.03] | 70.28% | [-0.02, 0.02] | 69.05%
## ind_Age_Str_Rel | 0.02 | [ 0.00, 0.05] | 99.40% | [-0.02, 0.02] | 35.43%
## ind_Age_Per_Rel | 0.01 | [-0.01, 0.04] | 84.92% | [-0.02, 0.02] | 72.17%
## ind_Age_Cau_Rel | 0.02 | [-0.01, 0.06] | 93.10% | [-0.02, 0.02] | 39.05%
## ind_Age_Dis_Rel | -9.98e-04 | [-0.02, 0.01] | 62.18% | [-0.02, 0.02] | 99.08%
## ind_Age_Rel | 0.05 | [-0.01, 0.12] | 94.35% | [-0.02, 0.02] | 13.50%
## tot_Age_Rel | 0.12 | [ 0.02, 0.23] | 98.88% | [-0.02, 0.02] | 2.50%Alrighty, let’s compare to what we obtained previously in lavaan
parest1 %>%
slice(55:75) %>%
select(label, est, pvalue, ci.lower, ci.upper) %>%
mutate(across(where(is.numeric), ~round(.x, 2))) %>%
tt(.)| label | est | pvalue | ci.lower | ci.upper |
|---|---|---|---|---|
| ind_x_m1_y1 | -0.01 | 0.38 | -0.05 | 0.02 |
| ind_x_m2_y1 | 0.03 | 0.03 | 0.00 | 0.06 |
| ind_x_m3_y1 | 0.01 | 0.52 | -0.02 | 0.03 |
| ind_x_m4_y1 | 0.05 | 0.01 | 0.01 | 0.08 |
| ind_x_m5_y1 | 0.04 | 0.02 | 0.00 | 0.07 |
| ind_x_y1 | 0.11 | 0.00 | 0.04 | 0.18 |
| tot_x_y1 | 0.27 | 0.00 | 0.16 | 0.39 |
| ind_x_m1_y2 | -0.04 | 0.02 | -0.08 | -0.01 |
| ind_x_m2_y2 | 0.01 | 0.29 | -0.01 | 0.03 |
| ind_x_m3_y2 | -0.01 | 0.61 | -0.03 | 0.02 |
| ind_x_m4_y2 | 0.01 | 0.46 | -0.02 | 0.04 |
| ind_x_m5_y2 | 0.02 | 0.35 | -0.02 | 0.05 |
| ind_x_y2 | -0.01 | 0.69 | -0.08 | 0.05 |
| tot_x_y2 | -0.07 | 0.22 | -0.19 | 0.04 |
| ind_x_m1_y3 | -0.01 | 0.60 | -0.04 | 0.03 |
| ind_x_m2_y3 | 0.03 | 0.04 | 0.00 | 0.05 |
| ind_x_m3_y3 | 0.01 | 0.32 | -0.01 | 0.04 |
| ind_x_m4_y3 | 0.02 | 0.14 | -0.01 | 0.06 |
| ind_x_m5_y3 | 0.01 | 0.71 | -0.03 | 0.04 |
| ind_x_y3 | 0.06 | 0.08 | -0.01 | 0.13 |
| tot_x_y3 | 0.13 | 0.03 | 0.02 | 0.25 |
The results that we got are very very similar and lead to the same interpretation.
However, before we interpret the estimates at all, it is best practice to check if our model converged properly. Let’s follow the steps of the WAMB checklist this time!
Resources
I am not a statistician, so please take all of what I wrote and coded here with a lot of caution. However, there are many great resources out there if you want to delve deeper into the topic - I am just getting started with this, but here are a few potentially helpful links.
Doing Bayesian Data Analysis in brms and the tidyverse: A brms and tidyverse version by Solomon Kurz of the classic book by Kruschke Doing Bayesian data analysis: A tutorial with R, JAGS, and Stan.
Practical Bayes Part I and Practical Bayes Part II: Great resources to dealing and avoiding common problems in Bayesian models by Micheal Clark
Bayesian Basics: General overview of how to run Bayesian models in R
bayestestR: Vignettes on basic Bayesian models and their interpretation
Stan Wiki: Prior Choice: Go to reference for selecting priors
Course Handouts for Bayesian Data Analysis Class by Mark Lai
And there are of course many many more resources, this is just to get started with some.