Bayesian RI-CLPM

[Aleija R.-R.]

Hi all,

I am trying to run an RI-CLPM in blavaan and have a handful of questions. I'll start with the "simplest", regarding the fit indices. Based on the 2020 paper, I know that gammahat and BCFI are the best to use, but I don't know if I should be concerned about the considerably "outside the cutoff" values for BRMSEA, adj gammahat, and BTLI. For model evaluation, how should I incorporate those other indices with the good-to-acceptable BCFI and gammahat values? For context, the sample size is ~8,000 with all variables having between 3 and 6 waves of data. Below is also the structure of my current model, which includes a mix of latent and observed variables in the cross-lagged paths, though other models I run with only observed variables have similar issues. Thanks so much.

Posterior summary statistics and highest posterior density (HPD)
 90% credible intervals for devm-based fit indices:

               EAP Median   MAP    SD lower upper
BRMSEA       0.097  0.097 0.097 0.000 0.096 0.097
BGammaHat    0.975  0.975 0.975 0.000 0.975 0.976
adjBGammaHat 0.754  0.754 0.754 0.002 0.751 0.757
BMc          0.859  0.859 0.859 0.001 0.857 0.860
BCFI         0.941  0.941 0.941 0.000 0.940 0.942
BTLI         0.683  0.683 0.684 0.002 0.680 0.687
BNFI         0.940  0.940 0.940 0.000 0.940 0.941

ri_clpm_model = '

########################################
# Mediator Measurement Model
########################################

# measurement model
MED_W2 =~ c1*indicator1_W2 + c2*indicator2_W2 + c3*indicator3_W2
MED_W4 =~ c1*indicator1_W4 + c2*indicator2_W4 + c3*indicator3_W4
MED_W6 =~ c1*indicator1_W6 + c2*indicator2_W6 + c3*indicator3_W6

# equal intercepts across time
indicator3_W2 ~  i1*1
indicator2_W2 ~  1
indicator1_W2 ~  i2*1

indicator3_W4 ~  i1*1
indicator2_W4 ~  1
indicator1_W4 ~  i3*1

indicator3_W6 ~  i1*1
indicator2_W6 ~  1
indicator1_W6 ~  i3*1


# latent means: fix baseline to 0, estimate later waves
MED_W2 ~ 0*1
MED_W4 ~ m4*1
MED_W6 ~ m6*1

# Correlated within-task residuals across time
indicator1_W6 ~~ indicator1_W4 + indicator1_W2
indicator1_W4 ~~ indicator1_W2

indicator3_W6 ~~ indicator3_W4 + indicator3_W2
indicator3_W4 ~~ indicator3_W2

indicator2_W6 ~~ indicator2_W4 + indicator2_W2
indicator2_W4 ~~ indicator2_W2

########################################
# Create between components (RIs)
########################################

RI_predictor =~ 1*predictor_W1 + 1*predictor_W2 + 1*predictor_W3 +
                1*predictor_W4 + 1*predictor_W5 + 1*predictor_W6

RI_outcome =~ 1*outcome_W1 + 1*outcome_W2 + 1*outcome_W3 +
              1*outcome_W4 + 1*outcome_W5 + 1*outcome_W6

RI_med2 =~ 1*mediator2_W2 + 1*mediator2_W4 + 1*mediator2_W6

RI_med1 =~ 1*MED_W2 + 1*MED_W4 + 1*MED_W6

########################################
# Set residual variances of all lvs to 0
########################################

MED_W2 ~~ 0*MED_W2
MED_W4 ~~ 0*MED_W4
MED_W6 ~~ 0*MED_W6

########################################
# Set residual variances of all ovs to 0
########################################

# Predictor observed variables (waves 1-6)
predictor_W1 ~~ 0*predictor_W1
predictor_W2 ~~ 0*predictor_W2
predictor_W3 ~~ 0*predictor_W3
predictor_W4 ~~ 0*predictor_W4
predictor_W5 ~~ 0*predictor_W5
predictor_W6 ~~ 0*predictor_W6

# Outcome observed variables (waves 1-6)
outcome_W1 ~~ 0*outcome_W1
outcome_W2 ~~ 0*outcome_W2
outcome_W3 ~~ 0*outcome_W3
outcome_W4 ~~ 0*outcome_W4
outcome_W5 ~~ 0*outcome_W5
outcome_W6 ~~ 0*outcome_W6

# Mediator2 observed variables (waves 2, 4, 6)
mediator2_W2 ~~ 0*mediator2_W2
mediator2_W4 ~~ 0*mediator2_W4
mediator2_W6 ~~ 0*mediator2_W6

########################################
# Create within components
########################################

# Predictor within-person
wp_predictor_W1 =~ 1*predictor_W1
wp_predictor_W2 =~ 1*predictor_W2
wp_predictor_W3 =~ 1*predictor_W3
wp_predictor_W4 =~ 1*predictor_W4
wp_predictor_W5 =~ 1*predictor_W5
wp_predictor_W6 =~ 1*predictor_W6

# Outcome within-person
wp_outcome_W1 =~ 1*outcome_W1
wp_outcome_W2 =~ 1*outcome_W2
wp_outcome_W3 =~ 1*outcome_W3
wp_outcome_W4 =~ 1*outcome_W4
wp_outcome_W5 =~ 1*outcome_W5
wp_outcome_W6 =~ 1*outcome_W6

# Mediator2 within-person
wp_med2_W2 =~ 1*mediator2_W2
wp_med2_W4 =~ 1*mediator2_W4
wp_med2_W6 =~ 1*mediator2_W6

# Mediator1 within-person
wp_med1_W2 =~ 1*MED_W2
wp_med1_W4 =~ 1*MED_W4
wp_med1_W6 =~ 1*MED_W6

########################################
# Estimate lagged effects between within-person centered variables
########################################

# Predictor → Outcome
wp_outcome_W3 ~ c1p*wp_predictor_W1
wp_outcome_W4 ~ wp_predictor_W2
wp_outcome_W5 ~ c2p*wp_predictor_W3
wp_outcome_W6 ~ wp_predictor_W4

# Outcome → Predictor
wp_predictor_W3 ~ c3p*wp_outcome_W1
wp_predictor_W4 ~ wp_outcome_W2
wp_predictor_W5 ~ c4p*wp_outcome_W3
wp_predictor_W6 ~ wp_outcome_W4

# Predictor/Outcome → Mediators at wave 2
wp_med2_W2 ~ a1*wp_predictor_W1 + a3*wp_outcome_W1
wp_med1_W2 ~ a2*wp_predictor_W1 + a4*wp_outcome_W1

# Wave 2 → Wave 3 (bidirectional)
wp_outcome_W3 ~ b1*wp_med2_W2 + b2*wp_med1_W2
wp_predictor_W3 ~ b3*wp_med2_W2 + b4*wp_med1_W2

# Mediation effects (wave 1→2→3)
ind_dir1 := a1*b1
ind_dir2 := a2*b2
tot_dir1 := c1p + (a1*b1) + (a2*b2)

ind_revdir1 := a3*b3
ind_revdir2 := a4*b4
tot_revdir1 := c3p + (a3*b3) + (a4*b4)

# Wave 3 → Wave 4
wp_med2_W4 ~ a5*wp_predictor_W3 + a7*wp_outcome_W3
wp_med1_W4 ~ a6*wp_predictor_W3 + a8*wp_outcome_W3

# Wave 4 → Wave 5 (bidirectional)
wp_outcome_W5 ~ b5*wp_med2_W4 + b6*wp_med1_W4
wp_predictor_W5 ~ b7*wp_med2_W4 + b8*wp_med1_W4

# Mediation effects (wave 3→4→5)
ind_dir3 := a5*b5
ind_dir4 := a6*b6
tot_dir2 := c2p + (a5*b5) + (a6*b6)

ind_revdir3 := a7*b7
ind_revdir4 := a8*b8
tot_revdir2 := c4p + (a7*b7) + (a8*b8)

# Wave 5 → Wave 6
wp_med2_W6 ~ wp_predictor_W5 + wp_outcome_W5
wp_med1_W6 ~ wp_predictor_W5 + wp_outcome_W5

########################################
# Estimate covariance between within-person variables at first wave
########################################

wp_predictor_W1 ~~ wp_outcome_W1
wp_med2_W2 ~~ wp_med1_W2 + wp_predictor_W2 + wp_outcome_W2
wp_med1_W2 ~~ wp_predictor_W2 + wp_outcome_W2

########################################
# Estimate covariances between residuals of within-person variables
########################################

wp_predictor_W2 ~~ wp_outcome_W2

wp_predictor_W3 ~~ wp_outcome_W3

wp_predictor_W4 ~~ wp_outcome_W4 + wp_med2_W4 + wp_med1_W4
wp_outcome_W4 ~~ wp_med2_W4 + wp_med1_W4
wp_med2_W4 ~~ wp_med1_W4

wp_predictor_W5 ~~ wp_outcome_W5

wp_predictor_W6 ~~ wp_outcome_W6 + wp_med2_W6 + wp_med1_W6
wp_outcome_W6 ~~ wp_med2_W6 + wp_med1_W6
wp_med2_W6 ~~ wp_med1_W6


########################################
# Estimate variance and covariance of random intercepts
########################################

RI_predictor ~~ RI_outcome + RI_med2 + RI_med1 + RI_predictor
RI_outcome ~~ RI_med2 + RI_med1 + RI_outcome
RI_med2 ~~ RI_med1 + RI_med2
RI_med1 ~~ RI_med1

########################################
# Estimate (residual) variance of within-person centered variables
########################################
wp_predictor_W1 ~~ wp_predictor_W1 # Variances
wp_outcome_W1 ~~ wp_outcome_W1

wp_predictor_W2 ~~ wp_predictor_W2
wp_outcome_W2 ~~ wp_outcome_W2
wp_med1_W2 ~~ wp_med1_W2
wp_med2_W2 ~~ wp_med2_W2

wp_predictor_W3 ~~ wp_predictor_W3 # Residual variances
wp_outcome_W3 ~~ wp_outcome_W3

wp_predictor_W4 ~~ wp_predictor_W4
wp_outcome_W4 ~~ wp_outcome_W4
wp_med1_W4 ~~ wp_med1_W4
wp_med2_W4 ~~ wp_med2_W4

wp_predictor_W5 ~~ wp_predictor_W5
wp_outcome_W5 ~~ wp_outcome_W5

wp_predictor_W6 ~~ wp_predictor_W6
wp_outcome_W6 ~~ wp_outcome_W6
wp_med1_W6 ~~ wp_med1_W6
wp_med2_W6 ~~ wp_med2_W6


########################################
# Within-person autoregressions
########################################

# Mediator2
wp_med2_W6 ~ wp_med2_W4
wp_med2_W4 ~ wp_med2_W2

# Mediator1
wp_med1_W6 ~ wp_med1_W4
wp_med1_W4 ~ wp_med1_W2

# Predictor
wp_predictor_W6 ~ wp_predictor_W5
wp_predictor_W5 ~ wp_predictor_W4
wp_predictor_W4 ~ wp_predictor_W3
wp_predictor_W3 ~ wp_predictor_W2
wp_predictor_W2 ~ wp_predictor_W1

# Outcome
wp_outcome_W6 ~ wp_outcome_W5
wp_outcome_W5 ~ wp_outcome_W4
wp_outcome_W4 ~ wp_outcome_W3
wp_outcome_W3 ~ wp_outcome_W2
wp_outcome_W2 ~ wp_outcome_W1

########################################
# Additional constraints
########################################

# Set correlations between RIs and within-person factors at exogenous wave at 0

RI_predictor ~~ 0*wp_predictor_W1
RI_outcome ~~ 0*wp_outcome_W1
RI_med2 ~~ 0*wp_med2_W2
RI_med1 ~~ 0*wp_med1_W2
'

cat("Model fit began at:", format(Sys.time()), "\n")

model_bayes = blavaan::bcfa(
  model    = ri_clpm_model,
  data     = my_data,
  n.chains = 3,
  burnin   = 1000,
  sample   = 5000,
  meanstructure = T,
  int.ov.free = F,
  bcontrol = list(cores = 3),
  control  = list(
    adapt_delta   = 0.9,
    max_treedepth = 10
  )
)

[Ed Merkle]

TK,

I can imagine that you will get different answers here depending on who you ask. A couple things that come to mind for me are below, though I acknowledge that they are not simple solutions:

- I think of the work of McNiesh and Wolf showing that "one size fits all" cutoff values don't work so well, and they recommend dynamic cutoffs (for example, see the paper below). I imagine that a similar thing applies to these Bayesian metrics, though the dynamic cutoffs will be more difficult to obtain.

https://doi.org/10.3758/s13428-022-01847-y

- I have recently found it useful to do model checking by incorporating the latent variables and treating the model more like a regression. Some examples are in my NCME workshop slides at the url below.

https://ecmerkle.github.io/blavaan/articles/resources.html

Mauricio sometimes lurks on here, and I'd also be interested to hear his opinion.

Ed

[Aleija R.-R.]

Thanks for the feedback and the references, Ed! A question related to my model, though not explicitly related to the above question, has to do with the test argument mentioned in the original blavaan paper. That paper mentions that for datasets with missing data, one can set test = "none" to avoid the PPP, which is not all that useful when missing data is in play. I did this with a model and it led to issues generating any sort of fit indices using blavFitIndices. This makes sense, but I was wondering if there was a way to eliminate the PPP calculations while still retaining the ability calculate global and incremental fit indices. Thanks for all that you and Mauricio do for the community!

[Mauricio Garnier-Villarreal]

TK

Following on Eds comments. I think this is a good example of understanding what each index is doing, and the problem of general cutoffs. From our paper an the one from Fan and Sivo, I wouldnt recommend to use RMSEA in general, as it is heavily affected by model size, so a larger model will generally be preffered, and a good small model would look "worse" than it actually is. I can imagine that a model like this (with a bunch of constraints) would be smaller and so be penalized by RMSEA. Of course the practical problem is that reviewers might asked about it. 

I would suggest to cite our paper and Fan an Sivo, for why only gammahat and CFI are reported, there is a section on our paper about the interpretation and use of fit indices, and how they were created as effect sizes instead of "pass/fail" tests and how we should go back to use them as effect sizes.  And check the Bayesian Modification indices ( Garnier-Villarreal, M., & Jorgensen, T. D. (2025 ), as a way to show that no model modification would improve the model substantially

I like the general concept of dynamic cutoffs, but in practice is hard to do, because you have to justify the decisions to the "modified conditions" for new cutoffs, so its easy to see it as cherry picking conditions that are favorable. And the package simsem was originally designed to do this dynamic cutoff as a montecarlo simulation

Hope this helps 

Fan, X., & Sivo, S. A. (2007). Sensitivity of fit indices to model misspecification and model types. Multivariate Behavioral Research, 42(3), 509–529. https://doi.org/10.1080/00273170701382864

Garnier-Villarreal, M., & Jorgensen, T. D. (2020). Adapting Fit Indices for Bayesian Structural Equation Modeling: Comparison to Maximum Likelihood. Psychological Methods, 25(1), 46–70. https://doi.org/10.1037/met0000224

Garnier-Villarreal, M., & Jorgensen, T. D. (2025). Evaluating Local Model Misspecification with Modification Indices in Bayesian Structural Equation Modeling. Structural Equation Modeling: A Multidisciplinary Journal, 32(2), 304–318. https://doi.org/10.1080/10705511.2024.2413128

[Aleija R.-R.]

Hi Mauricio,

Thanks for the reply! In referring to my current model as small, does this mean that the "size" of a model is tied more to the number of constraints it has than to the number of parameters and/or variables it has? I've typically thought of my current model as large because it has many parameters and variables, but your message suggests that my model is actually small because it has many constraints. Thanks for your input and for the references you posted above!

[Mauricio Garnier-Villarreal]

TK

For RMSEA, the size of the model is in function of the degrees of freedom. So models with higher df will have lower RMSEA. By having many constraints you will decrease the df of a model. I think you can find the model df from the object of the blavfitindices() function