Hi,

First of all, I want to thank Paul and everyone involved with the development for the amazing work done with brms. Due to its flexibility and ease of use brms has quickly become my favorite package for pretty much all statistical analyses.

Currently I am trying to do some model comparison between Structural Equation Models (SEMs) and some models that do not include a structural part (i.e. only direct effects) in brms. In other words, I am interested in whether the structural component adds relevant information to the model or not. Using LOOIC for this does not come about very well as the responses in the SEM and non-SEM models are different (there is at least one additional response in the SEMs) and I am wondering if anyone would have any suggestion of how could model comparison in such a situation be performed?

A â€śquick and dirtyâ€ť example of the problems using LOOIC:

library(brms)

x <- rnorm(100, mean=0, sd=1)

z <- 0.5x + rnorm(100, mean=0, sd=1)x + 0.5*z + rnorm(100, mean=0, sd=1)

y <- 0.5data <- data.frame(y=y, x=x, z=z)

zmod <- bf(z~x)

ymod1 <- bf(y~x+z)

ymod2 <- bf(y~z)yfit1 <- brm(ymod1 + zmod + set_rescor(FALSE),

family=gaussian,

data=data)yfit2 <- brm(ymod2 + zmod + set_rescor(FALSE),

family=gaussian,

data=data)yfit3 <- brm(ymod1,

family = gaussian,

data=data)

The models find the coefficients for each of the relationships and, at least in my test run, all coefficients in yfit1 differ from zero and standard deviations of one for both z and y are recognized:

Population-Level Effects:

Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat

y_Intercept 0.16 0.11 -0.04 0.38 4000 1.00

z_Intercept -0.07 0.10 -0.27 0.12 4000 1.00

y_x 0.35 0.15 0.07 0.63 4000 1.00

y_z 0.56 0.11 0.35 0.77 4000 1.00

z_x 0.66 0.12 0.43 0.90 4000 1.00

Family Specific Parameters:

Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat

sigma_y 1.06 0.08 0.93 1.24 4000 1.00

sigma_z 0.99 0.07 0.86 1.14 4000 1.00

Then the model comparison:

loo_yfit <- loo(yfit1,yfit2,yfit3)

loo_yfit

`LOOIC SE`

yfit1 581.05 18.59

yfit2 561.86 22.52

yfit3 298.78 14.70

yfit1 - yfit2 19.19 29.18

yfit1 - yfit3 282.27 12.67

yfit2 - yfit3 263.07 27.10

Comparison between yfit1 and yfit2 seems to be more-or-less OK, but yfit3 is obviously not comparable due to different response variables (or one additional response actually) and brms rightly warns about this:

Warning message:

Model comparisons are likely invalid as the response values of at least two models do not match.

I may be banging my head against the wall here, but is there any way of conducting some kind of formal model comparison in such a situation in brms, or am I just forced to stare at the Cr.I.s of the coefficients and conduct some stepwise selection method, or revert to some other package (e.g. lavaan)? I realize that SEMs are probably not among the core development directions of brms, but just thought someone might have come across the same issue and solved it one way or the other.

Operating System: Windows 7

brms Version: 2.3.0

Cheers,

Aapo K.