Comparing brms Univariate vs. Joint Models using PSIS-LOO. Is ELPD additive?

Modeling
, ,
1

Hi all!

First time poster, long time struggler :)

I am fairly new to Bayesian modelling (though I’ll be starting a PhD in it come October). I’m currently working on a modelling project using brms for my MSc dissertation, where I am analysing several continuous water quality parameters (e.g. phosphate, nitrate, ammonia), each with their own distribution, covariates, and smooth terms.

There are known dependencies between these pollutants, so I’m exploring whether treating them as conditionally independent (as in traditional univariate approaches) might be inadequate. To investigate this, I’ve:

  • Fit univariate models separately for each pollutant, and
  • Fit a joint model using brms with “set_rescor = FALSE”, where each outcome has its own likelihood and predictors (e.g., ammonia is included as a predictor in the nitrate submodel, reflecting known mechanistic links).

Now I’d like to compare the predictive performance of the joint and univariate approaches using psis_loo().

From what I understand, when set_rescor = FALSE, the joint log-likelihood is just the sum of the individual log-likelihoods per response. So my question is:

Is it valid to sum the individual loo() results from the univariate models and compare that total ELPD to the ELPD from the joint model?

I’ve seen suggestions that this might be appropriate when responses are conditionally independent (which I am assuming is the case with my univariate models), but I’m struggling to find any clear evidence that this would be suitable and I want to be confident before including this in my dissertation.

Any thoughts or corrections would be hugely appreciated!!

1 Like
2

Welcome to Stan discourse,

I have answer for this, but curious to first know more about what is the connection to the question Joint Models using PSIS-LOO. Is ELPD additive?? Two different first posters are asking the same question about the same data and models the same day. Is this some course assignment?

3

@avehtari Thanks for the welcome!

I was also very confused when I saw the other post, as it is quite literally just a section of my original post copied over.

My supervisor and I have designed this project ourselves for my thesis so no one else is working on this specific problem.

I had assumed it was probably just an accidental repost, maybe they were wanting to reply? I’m really not sure.

4

ok, here’s my answer

Expected log scores for independent predictions are additive. However, we need to be careful considering what is left out. When using brms mvbind multivariate model, in the joint model one observation for each outcome is dropped, and thus for example with two outcomes, you could also call this leave-one-pair-out. In independent models only one observation is left out the time, but as you can sum independent log scores, the result is the same as if you would drop two paired observations at the same time.

For example, using the brms vignette example Estimating Multivariate Models with brms • brms

library(brms)

data("BTdata", package = "MCMCglmm")
head(BTdata)

bform1 <- 
  bf(mvbind(tarsus, back) ~ sex + hatchdate) +
  set_rescor(FALSE)
fit1 <- brm(bform1, data = BTdata)
fit1 <- add_criterion(fit1, "loo")

fit1a <- brm(tarsus ~ sex + hatchdate, data = BTdata)
fit1b <- brm(back ~ sex + hatchdate, data = BTdata)
fit1a <- add_criterion(fit1a, "loo")
fit1b <- add_criterion(fit1b, "loo")

loo(fit1)
loo(fit1a)
loo(fit1b)

For all these models, LOO sees 828 observations, and for the multivariate model we’re leaving out two observations (one tarsus and one back observation), while in the independent models we’re leaving out one observation (either one tarsus or one back observation).

In the above simplified example, there are no shared parameters or priors in fit1, and thus the sum of elpd_loo’s for fit1a and fit1b is the same as elpd_loo for fit1. This laso means that, instead of making those separate models, you can use mvbind to create a multivariate model where the models are independent and you avoid the need to explicitly summing the elpd_loos.

By having a different formula with shared parameters or priors, you may see better elpd_loo for the joint model.

4 Likes
5

Sorry for the slow reply @avehtari, and thank you so much for your detailed response. Your explanation and the example have cleared this up for me and I very much appreciate it! Feeling a lot more confident with this approach now :)

1 Like