[Feature Request] tidy and glance methods?

[emstruong]

Hi Ed,

The tidy() method for blavaan doesn't seem to return confidence intervals. The glance() method also seems to be broken.

I was thinking that it might be a good idea to extend the method, such that it is a cross of the tidy method for lavaan and also the tidy method for brms.

Would you be open to a PR and are there any specific requests you have for the methods?

library(lavaan)
#> This is lavaan 0.6-20
#> lavaan is FREE software! Please report any bugs.
library(broom)
library(blavaan)
#> Loading required package: Rcpp
#> This is blavaan 0.5-9.1372
#> On multicore systems, we suggest use of future::plan("multicore") or
#>   future::plan("multisession") for faster post-MCMC computations.


HS.model <- ' visual  =~ x1 + x2 + x3
              textual =~ x4 + x5 + x6
              speed   =~ x7 + x8 + x9 '
     
fit <- blavaan::bcfa(HS.model, data=HolzingerSwineford1939, seed=959)
#> 
#> SAMPLING FOR MODEL 'stanmarg' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.000275 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.75 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 1500 [  0%]  (Warmup)
#> Chain 1: Iteration:  150 / 1500 [ 10%]  (Warmup)
#> Chain 1: Iteration:  300 / 1500 [ 20%]  (Warmup)
#> Chain 1: Iteration:  450 / 1500 [ 30%]  (Warmup)
#> Chain 1: Iteration:  501 / 1500 [ 33%]  (Sampling)
#> Chain 1: Iteration:  650 / 1500 [ 43%]  (Sampling)
#> Chain 1: Iteration:  800 / 1500 [ 53%]  (Sampling)
#> Chain 1: Iteration:  950 / 1500 [ 63%]  (Sampling)
#> Chain 1: Iteration: 1100 / 1500 [ 73%]  (Sampling)
#> Chain 1: Iteration: 1250 / 1500 [ 83%]  (Sampling)
#> Chain 1: Iteration: 1400 / 1500 [ 93%]  (Sampling)
#> Chain 1: Iteration: 1500 / 1500 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 1.676 seconds (Warm-up)
#> Chain 1:                3.302 seconds (Sampling)
#> Chain 1:                4.978 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'stanmarg' NOW (CHAIN 2).
#> Chain 2: 
#> Chain 2: Gradient evaluation took 0.000186 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.86 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2: 
#> Chain 2: 
#> Chain 2: Iteration:    1 / 1500 [  0%]  (Warmup)
#> Chain 2: Iteration:  150 / 1500 [ 10%]  (Warmup)
#> Chain 2: Iteration:  300 / 1500 [ 20%]  (Warmup)
#> Chain 2: Iteration:  450 / 1500 [ 30%]  (Warmup)
#> Chain 2: Iteration:  501 / 1500 [ 33%]  (Sampling)
#> Chain 2: Iteration:  650 / 1500 [ 43%]  (Sampling)
#> Chain 2: Iteration:  800 / 1500 [ 53%]  (Sampling)
#> Chain 2: Iteration:  950 / 1500 [ 63%]  (Sampling)
#> Chain 2: Iteration: 1100 / 1500 [ 73%]  (Sampling)
#> Chain 2: Iteration: 1250 / 1500 [ 83%]  (Sampling)
#> Chain 2: Iteration: 1400 / 1500 [ 93%]  (Sampling)
#> Chain 2: Iteration: 1500 / 1500 [100%]  (Sampling)
#> Chain 2: 
#> Chain 2:  Elapsed Time: 1.712 seconds (Warm-up)
#> Chain 2:                3.041 seconds (Sampling)
#> Chain 2:                4.753 seconds (Total)
#> Chain 2: 
#> 
#> SAMPLING FOR MODEL 'stanmarg' NOW (CHAIN 3).
#> Chain 3: 
#> Chain 3: Gradient evaluation took 0.000194 seconds
#> Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.94 seconds.
#> Chain 3: Adjust your expectations accordingly!
#> Chain 3: 
#> Chain 3: 
#> Chain 3: Iteration:    1 / 1500 [  0%]  (Warmup)
#> Chain 3: Iteration:  150 / 1500 [ 10%]  (Warmup)
#> Chain 3: Iteration:  300 / 1500 [ 20%]  (Warmup)
#> Chain 3: Iteration:  450 / 1500 [ 30%]  (Warmup)
#> Chain 3: Iteration:  501 / 1500 [ 33%]  (Sampling)
#> Chain 3: Iteration:  650 / 1500 [ 43%]  (Sampling)
#> Chain 3: Iteration:  800 / 1500 [ 53%]  (Sampling)
#> Chain 3: Iteration:  950 / 1500 [ 63%]  (Sampling)
#> Chain 3: Iteration: 1100 / 1500 [ 73%]  (Sampling)
#> Chain 3: Iteration: 1250 / 1500 [ 83%]  (Sampling)
#> Chain 3: Iteration: 1400 / 1500 [ 93%]  (Sampling)
#> Chain 3: Iteration: 1500 / 1500 [100%]  (Sampling)
#> Chain 3: 
#> Chain 3:  Elapsed Time: 1.971 seconds (Warm-up)
#> Chain 3:                3.447 seconds (Sampling)
#> Chain 3:                5.418 seconds (Total)
#> Chain 3: 
#> Computing post-estimation metrics (including lvs if requested)...

tidy(fit)
#> # A tibble: 24 × 5
#>    term          op    estimate std.lv std.all
#>    <chr>         <chr>    <dbl>  <dbl>   <dbl>
#>  1 visual =~ x1  =~       1      0.888   0.760
#>  2 visual =~ x2  =~       0.575  0.511   0.431
#>  3 visual =~ x3  =~       0.759  0.674   0.591
#>  4 textual =~ x4 =~       1      0.991   0.849
#>  5 textual =~ x5 =~       1.12   1.11    0.855
#>  6 textual =~ x6 =~       0.932  0.924   0.838
#>  7 speed =~ x7   =~       1      0.595   0.546
#>  8 speed =~ x8   =~       1.24   0.738   0.718
#>  9 speed =~ x9   =~       1.19   0.711   0.692
#> 10 x1 ~~ x1      ~~       0.578  0.578   0.423
#> # ℹ 14 more rows
tidy(fit, conf.int = TRUE)
#> # A tibble: 24 × 5
#>    term          op    estimate std.lv std.all
#>    <chr>         <chr>    <dbl>  <dbl>   <dbl>
#>  1 visual =~ x1  =~       1      0.888   0.760
#>  2 visual =~ x2  =~       0.575  0.511   0.431
#>  3 visual =~ x3  =~       0.759  0.674   0.591
#>  4 textual =~ x4 =~       1      0.991   0.849
#>  5 textual =~ x5 =~       1.12   1.11    0.855
#>  6 textual =~ x6 =~       0.932  0.924   0.838
#>  7 speed =~ x7   =~       1      0.595   0.546
#>  8 speed =~ x8   =~       1.24   0.738   0.718
#>  9 speed =~ x9   =~       1.19   0.711   0.692
#> 10 x1 ~~ x1      ~~       0.578  0.578   0.423
#> # ℹ 14 more rows
glance(fit)
#> Error in `rename()`:
#> ! Can't rename columns that don't exist.
#> ✖ Column `rmsea.ci.upper` doesn't exist.

^{Created on 2025-12-14 with reprex v2.1.1}

[ecmerkle]

Thanks for the suggestion, and it sounds like a good improvement!

I imagine the problem is that some column names are different, as compared to lavaan. For example, see colnames(summary(fit)). I am not sure that those columns should be renamed because it could break something else, but maybe duplicate columns could be created.

[ecmerkle]

Also, for something like rmsea: there is a Bayesian rmsea available via blavFitIndices(), but it would add a good deal of computation time for many models. I think we should try to avoid adding more computation time.

[emstruong]

Hi @ecmerkle ,

If you meant renaming what summary() itself returns, I agree with not renaming anything to prevent breaking.

But what about something like this (very quickly written for conceptual purposes).

tidy.blavaan <- function(BLAVAAN_MODEL, conf.int = TRUE, conf.level = 0.95) {

dataframe <- blavInspect(BLAVAAN_MODEL)

dataframe |>
dplyr::mutate(...) |>
dplyr::rename(...) |>
as_tibble()

}

In other words, let the vast majority of the computational aspects be handled through blavaan's internal functions. However, merely use the tidy() function to extract the relevant information and format it in a data structure that better matches what most tidy() functions in the broader broom() ecosystem, namely: brms and lavaan.

What we can do for fit indices via glance() is to add a user-facing argument for what fit indices to get and make the defaults something both informative and computationally cheap.

In case you are not already aware, perhaps these broom-specific glossaries could be helpful: https://www.tidymodels.org/learn/develop/broom/#glossaries

[ecmerkle]

I think I misunderstood the original post. I thought you wanted to modify blavaan so that the lavaan tidy methods also work with blavaan objects. Now I see that you want to implement new methods that are specific to blavaan.

I think this would be useful, so long as it isn't a ton of extra code to maintain and there aren't lots of extra dependencies. (With infinite time, I would try to reduce the number of blavaan dependencies. But I am guessing that the packages you would use are already indirect blavaan dependencies via the bayesplot package.)

[emstruong]

No worries :)

Yes, I'm pretty confident that I can keep the dependencies under control and stay largely within the realm of base R.

I only used dplyr in my example as I think that's the simplest way to explain my overall design.

Although, it would be nice if we could increment the minimum supported version of R from 3.5 to 4.1. This would let both me and new contributors use the base R pipe, which makes for an easy reading and writing experience.

[ecmerkle]

I am ok with bumping the minimum to 4.1, which I think is already needed due to other dependencies.

[emstruong]

HI @ecmerkle , happy belated new year.

Do you know if I have to install the package (devtools::install()) everytime I want to test if my changes work or if there's some workaround? Because I noticed that I can't actually use the loaded version of blavaan (devtools::load_all()) to fit CFA models, unless I install the codebase... I assume it's because of the need to compile the models...

In fact, is it possible to use cmdstanr instead of rstan when installing the package? I feel like that would make development a lot faster...

[ecmerkle]

Happy 2026! I typically use pkgload::load_all(., export_all = FALSE) during development. The first call to that will take a long time because it compiles the big Stan model. After that, it is fast unless you are directly making changes to the Stan model (then you may need recompile = TRUE).

There is also a "gotcha" involving debug builds, which will slow down the Stan model that is compiled unless you set an environment variable. See r-lib/pkgbuild#154

I will consider the move to cmdstanr some more, especially because rstan is not really being updated any more. It will place more responsibility on the user to compile the Stan file and also to get cmdstan on their system outside of R. I agree that it could improve development in some ways, though I also worry that it will lead to users having trouble with installing cmdstan and getting it recognized in R.

[emstruong]

Thanks @ecmerkle ! I'll give what you recommended a shot.

Oh, I was thinking of using cmdstanr solely for development purposes or for users 'in-the-know'. I think that RStan should probably be the default for most blavaan users, or else you'll probably get a lot of questions about, "Well, how do I get cmdstanr running in the first place?". I don't think it's that obvious.

I don't know much about RStan, but AFAICT, it seems to use a version of Stan that is many many versions behind the latest release. I think that's surprising. Even setting computational speed aside, that implies that the users of brms and blavaan who are not-in-the-know might suffer from older bugs for a very long time.

[ecmerkle]

Yes, the rstan issue often comes up in Stan Discourse. See https://discourse.mc-stan.org/t/stan-stickers/11385/14

[emstruong]

O.o

[emstruong]

Hi @ecmerkle , I found out that apparently load_all() now has a debug argument, so then I'm pretty sure there's no longer a need to manually set an environmental variable.

So then, some copy-pastable commands are

# First time compiling
pkgload::load_all(".", export_all = FALSE, compile = TRUE, debug = FALSE)

# Already compiled
pkgload::load_all(".", export_all = FALSE, debug = FALSE, compile = FALSE)

---

For full disclosure, I asked Claude Code to write up a draft function to handle everything we discussed here and I'm reviewing it myself now, line-by-line. I did this both because I'm a bit tied-up with other things and also because I was wondering whether it really could succeed.

If you have any concerns about this, please let me know.

[emstruong]

Hi Ed,

Regarding tidy(), I didn't like what Claude wrote, so I made things much simpler and made it depend directly on your summary() function.

Regarding glance(), I think what Claude wrote might be OK. Although, I do have a few questions...

Regarding RHat, I see that some recent research sets a 1.01 threshold. Although, I'm aware that older research sets a 1.05 threshold. Do you think a rhat of $\leq 1.01$ in blavaan should be used to determine convergence?

Also, I see that there seem to be two different fit index functions fitMeasures(), which is a re-branded form of blav_fit_measures(); and blavFitIndices().

For the purposes of this glance() function, do you think we should try to have all of fit indices we can get?

For what it's worth, this is what glance() does with lavaan()

library(lavaan)
#> This is lavaan 0.6-21
#> lavaan is FREE software! Please report any bugs.
library(broom)
# fit model
cfa.fit <- cfa(
  "F =~ x1 + x2 + x3 + x4 + x5",
  data = HolzingerSwineford1939, group = "school"
)

# summarize model fit with tidiers
glance(cfa.fit)
#> # A tibble: 1 × 17
#>    agfi   AIC   BIC   cfi chisq  npar rmsea rmsea.conf.high  srmr   tli
#>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>           <dbl> <dbl> <dbl>
#> 1 0.971 4473. 4584. 0.766  99.3    30 0.244           0.288 0.115 0.533
#> # ℹ 7 more variables: converged <lgl>, estimator <chr>, ngroups <int>,
#> #   missing_method <chr>, nobs <int>, norig <int>, nexcluded <int>

^{Created on 2026-02-17 with reprex v2.1.1}

And this is what tidy looks like for blavaan right now

pkgload::load_all(".", export_all = TRUE, debug = TRUE, compile = FALSE)
#> ℹ Loading blavaan
#> Loading required package: Rcpp
#> 
#> This is blavaan 0.5-10.1386
#> 
#> On multicore systems, we suggest use of future::plan("multicore") or
#>   future::plan("multisession") for faster post-MCMC computations.
data(HolzingerSwineford1939, package = "lavaan")

HS.model <- 'visual =~ x1 + x2 + x3'
fit <- bcfa(HS.model, data = HolzingerSwineford1939, seed = 123,
            n.chains = 1, sample = 300)
#> 
#> SAMPLING FOR MODEL 'stanmarg' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 9.6e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.96 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:   1 / 800 [  0%]  (Warmup)
#> Chain 1: Iteration:  80 / 800 [ 10%]  (Warmup)
#> Chain 1: Iteration: 160 / 800 [ 20%]  (Warmup)
#> Chain 1: Iteration: 240 / 800 [ 30%]  (Warmup)
#> Chain 1: Iteration: 320 / 800 [ 40%]  (Warmup)
#> Chain 1: Iteration: 400 / 800 [ 50%]  (Warmup)
#> Chain 1: Iteration: 480 / 800 [ 60%]  (Warmup)
#> Chain 1: Iteration: 501 / 800 [ 62%]  (Sampling)
#> Chain 1: Iteration: 580 / 800 [ 72%]  (Sampling)
#> Chain 1: Iteration: 660 / 800 [ 82%]  (Sampling)
#> Chain 1: Iteration: 740 / 800 [ 92%]  (Sampling)
#> Chain 1: Iteration: 800 / 800 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.415 seconds (Warm-up)
#> Chain 1:                0.187 seconds (Sampling)
#> Chain 1:                0.602 seconds (Total)
#> Chain 1: 
#> Computing post-estimation metrics (including lvs if requested)...
tidy(fit)
#>             term op  estimate  std.error conf.low conf.high  rhat      ess
#> 1     visual=~x1 =~ 1.0000000 0.00000000       NA        NA    NA       NA
#> 2     visual=~x2 =~ 0.8211175 0.14140096    0.566     1.097 0.999 218.8016
#> 3     visual=~x3 =~ 1.2327730 0.24051891    0.846     1.801 1.007 107.6644
#> 4         x1~~x1 ~~ 0.8851181 0.11774344    0.663     1.136 1.002 165.3633
#> 5         x2~~x2 ~~ 1.0759911 0.09874856    0.906     1.288 0.997 126.6143
#> 6         x3~~x3 ~~ 0.5868126 0.14518098    0.273     0.866 0.997 112.2273
#> 7 visual~~visual ~~ 0.4772758 0.11255847    0.262     0.715 1.006 164.1344
#>             prior
#> 1                
#> 2    normal(0,10)
#> 3    normal(0,10)
#> 4 gamma(1,.5)[sd]
#> 5 gamma(1,.5)[sd]
#> 6 gamma(1,.5)[sd]
#> 7 gamma(1,.5)[sd]

^{Created on 2026-02-17 with reprex v2.1.1}

[ecmerkle]

It is looking good! Some responses:

About rhat: I think that 1.01 should be fine for target="stan" and cmdstan. I need to look back at target = "jags", which is less efficient and which may use older Rhat computations.

About fit indices: the blavFitIndices() command is useful but typically takes a good deal of additional time to run after (and separately from) model estimation. I am not sure you want a user to estimate a model, then the glance command takes minutes or longer because it is running blavFitIndices() in the background. In that case, the user expects to glance but they instead get a long hard stare!

[emstruong]

OK, well I'll have to pause again, unfortunately. But I'll try to get back to this.

All that's left is revising:

• Claude's glance code,
• Claude's tidy documentation
• Claude's glance documentation

Then I can submit a PR.