Groups keyboard shortcuts have been updated
Dismiss
See shortcuts

HPD for models with indirect and total effects

18 views
Skip to first unread message

Pedro Ribeiro

unread,
Jun 25, 2024, 8:16:49 PM6/25/24
to blavaan
Hi everyone,

I'm estimating a model in blavaan (version 0.5-5) with indirect and total effects and I have been unable to retrieve the highest posterior density (HPD) interval of the parameters. I'm sending below a token example taken from a previous post on this group:

set.seed(1234)
X <- rnorm(100)
M <- 0.5*X + rnorm(100)
Y <- 0.7*M + rnorm(100)
Data <- data.frame(X = X, Y = Y, M = M)
model <- ' # direct effect
                 Y ~ c*X
               # mediator
                 M ~ a*X
                 Y ~ b*M
               # indirect effect (a*b)
                 ab := a*b
               # total effect
                 total := c + (a*b)'

fit <- bsem(model, data = Data,
            n.chains = 2, burnin = 1000, sample = 1000)
summary(fit)

blavInspect(fit, what = "hpd")

I'm getting the error:

Error in dimnames(x) <- dn : length of 'dimnames' [2] not equal to array extent

The command blavInspect(fit, what = "hpd") works when the indirect and total effects are not included. I was wondering whether this is a bug but please let me know if I'm missing something here! Thanks once again everyone for this amazing package!

Cheers,
Pedro

Mauricio Garnier-Villarreal

unread,
Jun 26, 2024, 6:46:43 AM6/26/24
to blavaan
Pedro

I think this is a bug related to having defined parameters and the new parameter labels. As the HPD is calculated from the MCMC draws, but the defined parameters are calculated outside of it.

For now, you can get the hpd for the all the parameters, except the defined ones from blavaan like this

blavInspect(fit, what = "hpd", add.labels = FALSE)

Or you can calculate them from the posterior draws and get the HPD like this

parTable(fit)[,c("lhs","op","rhs","label","pxnames")]

mc_out <- data.frame(as.matrix(blavInspect(fit, "mcmc", add.labels = FALSE)))
dim(mc_out)
head(mc_out)

mc_out$ab <- mc_out$bet_sign.3. * mc_out$bet_sign.1.
mc_out$total <- (mc_out$bet_sign.3. * mc_out$bet_sign.1.) + mc_out$bet_sign.2.

head(mc_out)

coda:::HPDinterval(as.mcmc(mc_out) )


Hope this helps

Pedro Ribeiro

unread,
Jun 27, 2024, 3:19:45 AM6/27/24
to blavaan
Hi Mauricio,

Thanks so much, this is really helpful!

Cheers,
Pedro
Reply all
Reply to author
Forward
0 new messages