Null model in blavaan?

[Mpop]

Latest version of JASP can do Bayesian SEM, which I think it does atop of blavaan, which can be found in its Process module. It outputs BIC for one thing. Problem is that BIC is a relative measure and there is no option to do any null model to compare it to. And so I want to do a null model in blavaan to compare my model to (comparing BICs). My model details the interrelations between 4 variables each with 11 data points. For the corresponding null model I have the following. It works up to a point but it isn't managing to pull out needed information to compute the BIC, I think because it isn't using the correct syntax to pull this info out. Can anyone please help? 

Plus any comment(s) on the wisdom, or lack thereof, of using this as a null model welcomed. 

Thank you. 

-------------------------------

# Install and load the required package
if (!require(blavaan)) install.packages("blavaan", dependencies = TRUE)
library(blavaan)

# Example dataset with 4 variables and 11 observations
set.seed(123)
data <- data.frame(
  var1 = sample(1:100, 11, replace = TRUE),
  var2 = sample(1:100, 11, replace = TRUE),
  var3 = sample(1:100, 11, replace = TRUE),
  var4 = sample(1:100, 11, replace = TRUE)
)

# Define a simple null model (no relationships)
null_model <- '
  var1 ~~ var1
  var2 ~~ var2
  var3 ~~ var3
  var4 ~~ var4
'

# Fit the model using Bayesian estimation with blavaan
fit_null <- blavaan::blavaan(null_model, data = data)

# Print model fit summary
summary_fit_null <- summary(fit_null)
print(summary_fit_null)

# Extract the necessary values for BIC and PPP
# BIC and PPP are part of the fit object summary
marginal_log_likelihood <- summary_fit_null$fit["MargLogLik"]
PPP_value <- summary_fit_null$fit["PPP"]

# Number of observations
n <- nrow(data)

# Number of parameters: 4 variances (one for each variable)
k <- 4  # For this simple null model with variances

# Compute Bayesian BIC (from marginal log-likelihood)
BIC_bayesian <- log(n) * k - 2 * marginal_log_likelihood
cat("Bayesian BIC for the null model:", BIC_bayesian, "\n")

# Output the Posterior Predictive P-value (PPP)
cat("Posterior Predictive P-value (PPP) for the null model:", PPP_value, "\n")

--------------------------------

[Mauricio Garnier-Villarreal]

hi

These indices are automatically calculated  by blavaan with

library(lavaan)
fitMeasures(fit_null)

Also, I dont recommend to use BIC for model comparison here, as they dont account for the variability of the posterior distributions properly. So, in this tutorial you can see how I recommend to do model comparison

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

[Mpop]

Thank you. 

Have amended the code, which seems to be working now, what do you think of it now (see below)? 

I know it reports some indices more than once but that is ok for me. 

---------------------------

# Install and load the required package
if (!require(blavaan)) install.packages("blavaan", dependencies = TRUE)
library(blavaan)

# Example dataset with 4 variables and 11 observations
set.seed(123)
data <- data.frame(
  var1 = sample(1:100, 11, replace = TRUE),
  var2 = sample(1:100, 11, replace = TRUE),
  var3 = sample(1:100, 11, replace = TRUE),
  var4 = sample(1:100, 11, replace = TRUE)
)

# Define a simple null model (no relationships)
null_model <- '
  var1 ~~ var1
  var2 ~~ var2
  var3 ~~ var3
  var4 ~~ var4
'
# Fit the model using Bayesian estimation with blavaan

fit_null <- blavaan(null_model, data = data)

# Print model fit summary
summary_fit_null <- summary(fit_null)
print(summary_fit_null)

# Report additional fit indices (WAIC, BIC, etc.)
fit_measures <- fitmeasures(fit_null, fit.measures = c("waic", "bic", "margloglik"))
print(fit_measures)

# Use lavaan (which blavaan builds upon) to get a collection of indices reported in a group
library(lavaan)
fitMeasures(fit_null)

----------------------------

I'm most interested in BIC because it is my understanding that knowing that for the alternative hypothesis and null models I can then get a BF10 value (using =exp((BIC_differential)/-2)). Which I can then multiply with another BF10 factor I have from another experiment. 

Moreover, because when I run the alternative hypothesis model in JASP I get the following error messages with the returned results: 

Warning: WAIC estimate unreliable - at least one effective parameter estimate (p_waic) larger than 0.4. We recommend using LOO instead.

Warning: LOO estimate unreliable - at least one observation with shape parameter (k) of the generalized Pareto distribution higher than 0.5. 

But I've been told by EJ (one of developers of JASP) that in that case the BIC should still be reliable. 

Of course, any comment(s) on anything here very gratefully received and appreciated. 

Thank you again.

[Mauricio Garnier-Villarreal]

Hi

There are 2 sides in the discussion about Bayesian hypothesis

- The Bayes Factor, which seeks to reduce everything to a BF test. JASP is designed from this perspectiev

- The posterior distribution side, which seeks to properly defined posterior distributions and ask questions to it

Here is where EJ and I disagree, I clearly side on the posterior distribution side. You can see the critics of the BF in these references. The BF needs strong priors to be reliable, and on large models with many parameters , like SEM, its hard to find and justify priors for the whole model. I see BF as a specific tool, instead of general bayesian comparison

Schad, D. J., Nicenboim, B., Bürkner, P.-C., Betancourt, M., & Vasishth, S. (2022). Workflow techniques for the robust use of bayes factors. Psychological Methods. https://doi.org/10.1037/met0000472

Tendeiro, J. N., & Kiers, H. A. L. (2019). A review of issues about null hypothesis Bayesian testing. Psychological Methods, 24(6), 774–795. https://doi.org/10.1037/met0000221

In the link I added before you can see how I recommend to do BSEM model comparison. I prefer the LOO or WAIC methods as they are general comparison methods that account for the full posterior distrubution and is less sensitive to priors. Can read more about it in these references

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

Vehtari, A., Gelman, A., & Gabry, J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing, 27(5), 1413–1432. https://doi.org/10.1007/s11222-016-9696-4

Vehtari, A., Mononen, T., Tolvanen, V., Sivula, T., & Winther, O. (n.d.). Bayesian Leave-One-Out Cross-Validation Approximations for Gaussian Latent Variable Models.

Further, I am not clear how you are using the null model here. The null model is needed to calcluate some approximate fit indices, but it seems odd to use it directly in model comparison. You can see how to calucltae approximate fit indices in this tutorial

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

You need to get the full output of the LOO calcluation to see how many observations have problematic pareto k estimates

[Ed Merkle]

A couple more thoughts here: if you are specifically interested in the Bayes factor, I would use the margloglik metric from fitMeasures() instead of BIC. Some detail about how blavaan does the approximation are at the end of the appendix of this paper:

https://www.jstatsoft.org/article/view/v085i04

Additionally, I find that the WAIC warnings sometimes go away if you use priors that are specific to your model, as opposed to using the software defaults. Some details about setting priors are here:

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

Finally, if you are considering using WAIC, maybe also consider DIC (which is reported by blavaan), whose computations are simpler than WAIC but improved over BIC.

Ed

--

You received this message because you are subscribed to the Google Groups "blavaan" group.

To unsubscribe from this group and stop receiving emails from it, send an email to blavaan+u...@googlegroups.com.

To view this discussion visit https://groups.google.com/d/msgid/blavaan/a2a5d22b-bfb7-4067-8ca0-f62de0d5974fn%40googlegroups.com.

[Mpop]

So I'm running the alternative hypothesis model in JASP (which runs blavaan). And the null model in blavaan using R code (shown above in earlier message). Why I'm using BIC as oppossed to WAIC or LOO is that when I run the alternative hypothesis model in JASP (with default parameters) I get the following error messages, wherein it is my (very amateur) understanding that the BIC value outputted is still reliable in this instance. So, I'm actually kind of forced to use BIC here. Of course, any comment(s) upon this very welcomed and grateful for. 

Warning: WAIC estimate unreliable - at least one effective parameter estimate (p_waic) larger than 0.4. We recommend using LOO instead.

Warning: LOO estimate unreliable - at least one observation with shape parameter (k) of the generalized Pareto distribution higher than 0.5. 

[Mauricio Garnier-Villarreal]

I still am not clear what you are trying to do. For specific JASP issues we cannot help, and might be that my recommendation would be to leave JASP and do all the analysis in R syntax

- Could you clarify why do you want to compare your model with this null model? This is not a common meaningful comparison. But it can be useful for estimation of some approximate fit indices

This warning is not a critical error, but a warning that requires you check further the estimation, or adjust the model. As Ed said, this could go away by adding better priors

Warning: LOO estimate unreliable - at least one observation with shape parameter (k) of the generalized Pareto distribution higher than 0.5. 

We also need to see more details about this, can you run this code and copy here the output, so we can see the the full warning report from LOO?

Where fit is your blavaan model.

library(blavaan)


blav_loo <- function (object1, ...) {
 
  require(loo)
 
  lavopt1 <- blavInspect(object1, "Options")
  if ((lavopt1$test == "none" & lavopt1$target != "stan") ) {
    stop("blavaan ERROR: loo cannot be estimated when test='none'")
  }
  targ1 <- lavopt1$target

  if (targ1 == "stan" && blavInspect(object1, "meanstructure")) {
    ll1 <- loo::extract_log_lik(object1@external$mcmcout)
  }
  else if (blavInspect(object1, "categorical") && lavopt1$test !=
           "none") {
    if ("llnsamp" %in% names(lavopt1)) {
      cat("blavaan NOTE: These criteria involve likelihood approximations that may be imprecise.\n",
          "You could try running the model again to see how much the criteria fluctuate.\n",
          "You can also manually set llnsamp for greater accuracy (but also greater runtime).\n\n")
    }
    ll1 <- object1@external$casells
  }
  else {
    lavopt1$estimator <- "ML"
    ll1 <- blavaan:::case_lls(object1@external$mcmcout,
                              blavaan:::make_mcmc(object1@external$mcmcout),
                    object1)
  }
  nchain1 <- blavInspect(object1, "n.chains")
  niter1 <- nrow(ll1)/nchain1
  cid1 <- rep(1:nchain1, each = niter1)
  ref1 <- relative_eff(exp(ll1), chain_id = cid1)

  loo1 <- loo(ll1, r_eff = ref1, ...)
  waic1 <- waic(ll1)
 
  return(list(loo=loo1, waic=waic1))

}


blav_loo(fit)

[Mpop]

RE: " Could you clarify why do you want to compare your model with this null model?"

My mission is to get to a BF10 value for this BSEM model which can be got from its BIC value as compared to the BIC value of a null model. That is the point of the null model. 

[Mauricio Garnier-Villarreal]

This is not how I would recommend to use these models.

But, if you want to focus on it, yes, the BIC would still be applicable/comparable