Comparing models with different functional forms of the same y variable using loo_compare

llewmills · February 19, 2024, 10:03am

I have modelled the effect of a cannabis agonist drug on number of days of cannabis use during a 12-week clinical trial. It is a mixed effects model for repeated measures regression, with group (time invariant categorical predictor: placebo vs active) as the primary predictor and time as the repeated measure (levels: baseline, Week 4, Week 8, Week 12). The outcome is number of days of illicit cannabis in the previous 28 days (measured four times as stated).

I can treat the outcome as either continuous, for which I would use a standard gaussian regression or as a bounded count in which case i would use an aggregated binomial. In brms the models looked like this

Gaussian

brm(formula = cu ~ group*week + (1|id),
 	  data = cu,
  	  family = gaussian(),
          prior = c(prior(normal(14, 1.5), class = Intercept),
 			  prior(normal(0, 11), class = b),
 			  prior(cauchy(1,2), class = sd)),
	 iter = 1e4,
 	 warmup = 1e3,
 	 chains = 3,
 	 cores = 3,
 	 control = list(adapt_delta = .99),
	 save_pars = save_pars(all=T),
 	 seed = 1234) -> modCont

Aggregated binomial

brm(formula = cu | trials(set)  ~ group*week + (1|id),
	data = cu,
	binomial(link = logit),
	prior = c(prior(normal(0, 1.5), class = Intercept),
			  prior(normal(0, 1), class = b),
		      prior(cauchy(0,2), class = sd)),
    iter = 1e4,
	warmup = 1e3,
	chains = 3,
	cores = 3,
	control = list(adapt_delta = .99),
	save_pars = save_pars(all=T),
	seed = 1234) -> modAggBinom

Now I used loo with moment matching on each model and hence was able to compare the two models with loo_compare and got the following output and warning message

            elpd_diff se_diff
modCont        0.0       0.0 
modAggBinom -813.4     142.6 
Warning message:
Not all models have the same y variable. ('yhash' attributes do not match)

Now to me this makes sense; you shouldn’t be able to compare models where the outcome variables are in a different form.

BUT the reason I did this is that I have a memory from either a lecture on youtube or a paper where Aki Vehtari suggested that with loo-cv it is possible to compare models where the outcome has different functional form. Did I dream that? Is there a way to compare the performance of my aggregated binomial vs gaussian models? If so what is the procedure? Also if so, where could I find some guidance on how to achieve such a comparison?

avehtari · February 20, 2024, 12:50pm

It is possible to compare continuous and discrete if you are careful, se CV-FAQ: Can cross-validation be used to compare different observation models / response distributions / likelihoods?

It would be useful to see posterior predictive check plots, too

llewmills · February 21, 2024, 4:10am

Thank you @avehtari. So if I understand what it says in the link correctly

(1) What I did in the example I gave is not correct, i.e. directly comparing the continuous vs discrete models via loo_compare?

(2) In order to compare the two I would need to discretise the continuous outcome? If this is the case (a) in what sense are they being compared given that iby discretising the continuous outcomes it will no longer be continuous (b) do you know of an online tutorial where this is done so I can follow along and learn how to do it?

avehtari · February 21, 2024, 8:00am

I have in my to-do list to make a case study illustrating the discretization. The easiest case is when the target is integer counts and the target is not scaled for the continuous model, as then in the discretization the width of each bin is 1, and the approximate probability of bin is the density at the target value times the bin width 1. In that case the density value is also (approximate) probability value and the models can be compared directly. Based on that, it seems like your models should be comparable if the data cu is the same. It is possible that the hash is different if brms stores the data in one model as doubles and in another one as integers. As the elpd_diff is very big, I’m suspicious of the results, but if you post also a plot of the distribution of y and some posterior predictive plots, I might be able to say more. The explanation and the code you posted is not clear on what do you mean by aggregated.

llewmills · February 21, 2024, 9:32am

I can’t edit my original post but here is the data. It is real data (already plublished here) from a 12 week clinical trial testing efficacy of a cannabis agonist drug nabiximols on cannabis use among people with cannabis dependence.

id is participant id. 128-level factor.
cu is the outcome, days of illicit cannabis use in the previous 28 days. It was measured four times, at weeks 0 (baseline), 4, 8, and 12. Not all participants stayed in the trial the whole 12 weeks so some have missing data. In the two models I described in the original post this was treated as either a continuous variable (in the model modCont) or as a bounded count (in the model modAggBinom).
set is the possible number of days of cannabis use. All values are therefore 28. This variable is used in the binomial regression syntax for the modAggBinom model described in the original post, i.e. cu | trials(set)
group is the group participants were randomised to, a binary categorical (placebo vs nabiximols).
week is a four-level factor indicating which stage of the trial the measure of cannabis use cu was taken: 0, 4, 8, or 12.

If you wanted to use this data for your case study I would be honoured (your to-do list permitting). This is real data from a published clinical trial and is simple.

Here is the data

id <- factor(c(1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 4, 5, 5, 6, 6, 6, 6, 7, 7, 8, 8, 8, 8, 9, 9, 9, 10, 10, 10, 10, 11, 11, 11, 12, 12, 13, 14, 15, 16, 16, 17, 18, 18, 18, 18, 19, 20, 20, 20, 20, 21, 21, 21, 21, 22, 22, 23, 23, 23, 24, 24, 24, 24, 25, 25, 25, 25, 26, 27, 27, 28, 28, 28, 28, 29, 30, 30, 30, 30, 31, 31, 32, 32, 32, 32, 33, 33, 33, 34, 34, 34, 35, 35, 36, 36, 37, 37, 37, 37, 38, 39, 39, 39, 39, 40, 40, 40, 41, 42, 42, 42, 42, 43, 43, 43, 43, 44, 44, 45, 45, 46, 46, 46, 46, 47, 47, 47, 47, 48, 48, 49, 49, 49, 50, 50, 50, 50, 51, 51, 51, 52, 52, 52, 52, 53, 53, 53, 53, 54, 54, 55, 55, 55, 55, 56, 57, 57, 57, 57, 58, 58, 58, 58, 59, 59, 59, 59, 60, 60, 60, 60, 61, 61, 61, 62, 63, 63, 64, 64, 64, 65, 65, 65, 65, 66, 66, 66, 66, 67, 67, 67, 67, 68, 68, 68, 69, 69, 69, 69, 70, 70, 70, 70, 71, 71, 71, 71, 72, 73, 73, 73, 73, 74, 74, 74, 75, 76, 76, 76, 76, 77, 77, 77, 77, 78, 78, 78, 79, 79, 79, 79, 80, 80, 80, 80, 81, 81, 81, 81, 82, 82, 83, 83, 84, 84, 84, 85, 85, 85, 86, 86, 86, 86, 87, 87, 87, 87, 88, 88, 88, 88, 89, 89, 89, 89, 90, 90, 90, 90, 91, 91, 91, 91, 92, 92, 92, 92, 93, 93, 93, 93, 94, 94, 94, 94, 95, 95, 95, 95, 96, 96, 96, 96, 97, 97, 97, 98, 98, 98, 98, 99, 99, 99, 99, 100, 101, 101, 101, 102, 102, 102, 102, 103, 103, 103, 103, 104, 104, 105, 105, 105, 105, 106, 106, 106, 106, 107, 107, 107, 107, 108, 108, 108, 108, 109, 109, 109, 109, 110, 110, 111, 111, 112, 112, 112, 112, 113, 113, 113, 113, 114, 115, 115, 115, 115, 116, 116, 116, 116, 117, 117, 117, 117, 118, 118, 119, 119, 119, 119, 120, 120, 120, 120, 121, 121, 121, 122, 123, 123, 123, 123, 124, 124, 124, 125, 125, 125, 125, 126, 126, 126, 126, 127, 127, 128))

group <- factor(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0),
				levels = 0:1,
				labels = c("placebo", "nabiximols"))

week <- factor(c(0, 4, 8, 12, 0, 4, 8, 0, 4, 8, 0, 4, 8, 12, 0, 4, 0, 4, 8, 12, 0, 4, 0, 4, 8, 12, 0, 4, 8, 0, 4, 8, 12, 0, 4, 8, 0, 4, 0, 0, 0, 0, 4, 0, 0, 4, 8, 12, 0, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 0, 4, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 0, 4, 0, 4, 8, 12, 0, 0, 4, 8, 12, 0, 4, 0, 4, 8, 12, 0, 4, 8, 0, 4, 12, 0, 8, 0, 4, 0, 4, 8, 12, 0, 0, 4, 8, 12, 0, 4, 8, 0, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 0, 4, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 0, 4, 12, 0, 4, 8, 12, 0, 4, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 0, 4, 8, 12, 0, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 12, 0, 0, 4, 0, 4, 8, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 0, 4, 8, 12, 0, 4, 8, 0, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 0, 4, 0, 4, 8, 0, 4, 8, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 0, 4, 8, 12, 0, 4, 8, 12, 0, 0, 4, 8, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 0, 4, 0, 4, 8, 12, 0, 4, 8, 12, 0, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 12, 0, 0, 4, 8, 12, 0, 4, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 0))

cu <- c(13, 12, 12, 12, 28, 0, NA, 16, 9, 2, 28, 28, 28, 28, 28, NA, 28, 28, 17, 28, 28, NA, 16, 0, 0, NA, 28, 28, 28, 28, 17, 0, NA, 28, 27, 28, 28, 26, 24, 28, 28, 28, 25, 28, 26, 28, 18, 16, 28, 28, 7, 0, 2, 28, 2, 4, 1, 28, 28, 16, 28, 28, 24, 26, 15, 28, 25, 17, 1, 8, 28, 24, 27, 28, 28, 28, 28, 28, 27, 28, 28, 28, 28, 20, 28, 28, 28, 28, 12, 28, NA, 17, 15, 14, 28, 0, 28, 28, 28, 0, 0, 0, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 21, 24, 28, 27, 28, 28, 26, NA, 28, NA, 20, 2, 3, 7, 28, 1, 19, 8, 21, 7, 28, 28, 20, 28, 28, 28, 24, 20, 17, 11, 25, 25, 28, 26, 28, 24, 17, 16, 27, 14, 28, 28, 28, 28, 28, 28, 14, 13, 4, 24, 28, 28, 28, 21, 28, 21, 26, 28, 28, 0, 0, 28, 23, 20, 28, 20, 16, 28, 28, 28, 10, 1, 1, 2, 28, 28, 28, 28, 18, 22, 9, 15, 28, 9, 1, 20, 18, 20, 24, 28, 28, 28, 19, 28, 28, 28, 28, 28, 28, 28, 28, 28, 4, 14, 20, 28, 28, 0, 0, 0, 28, 20, 9, 24, 28, 28, 28, 28, 28, 21, 28, 28, 14, 24, 28, 23, 0, 0, 0, 28, NA, 28, NA, 28, 15, NA, 12, 25, NA, 28, 2, 0, 0, 28, 10, 0, 0, 28, 0, 0, 0, 23, 0, 0, 0, 28, 0, 0, 0, 28, 0, 0, 0, 28, 2, 1, 0, 21, 14, 7, 8, 28, 28, 28, 0, 28, 28, 20, 18, 24, 0, 0, 0, 28, 15, NA, 28, 1, 1, 2, 28, 1, 0, 0, 28, 28, 14, 21, 25, 19, 16, 13, 28, 28, 28, 28, 28, 28, 28, 27, 19, 21, 18, 1, 0, 0, 28, 28, 28, 28, 28, 24, 27, 28, 18, 0, 3, 8, 28, 28, 28, 9, 20, 25, 20, 12, 19, 0, 0, 0, 27, 28, 0, 0, 0, 20, 17, 16, 14, 28, 7, 0, 1, 28, 24, 28, 25, 23, 20, 28, 14, 16, 7, 28, 28, 26, 28, 28, 26, 28, 28, 28, 24, 20, 28, 28, 28, 28, 28, 8, 6, 4, 28, 20, 28)

set <- rep(28, length(cu))

cu <- data.frame(id, group, week, cu, set)

The brms code I provided in the original post should work on this dataset without need of any alteration.

llewmills · February 21, 2024, 9:36am

I use the term aggregated binomial regression but I am starting to realise few other people use this term. It is simply a binomial regression where the outcome is a bounded count (days’ illicit cannabis use out of 28) rather than a bernoulli [0,1] outcome.

I’m not sure I have the skills to implement your discretization suggestion without some guidance, not without a few days’ or weeks` reading and some practice at the very least. Personally I would find it incredibly useful to be able to do this. Days of substance use data is almost never distributed normally (most often bimodal, with modes at 0 days used/28 and 28 days used/28). It makes analysis quite challenging.

avehtari · February 21, 2024, 10:06am

daysCannUse is not defined, is it the same as cu?

llewmills · February 21, 2024, 10:30am

Yes that’s right. i’m sorry about that. Have edited the post with the data (on my phone while walking the dog). Should work now.

avehtari · February 21, 2024, 10:56am

Based on a quick experiment, you should have used family=beta_binomial(). I’ll get back to you with more details and explanation of comparison to continuous model.

llewmills · February 21, 2024, 11:22am

Amazing thank you!

avehtari · February 26, 2024, 6:40pm

Here’s a link to a case study that illustrates the discretization and I couldn’t resist to add additional workflow for model and prior checking https://users.aalto.fi/~ave/casestudies/Nabiximols/nabiximols.html

Let me know if the explanation is clear enough, and I can add more if needed

llewmills · February 26, 2024, 9:40pm

Absolutely increible @avehtari! So thorough. You never need to resist the urge to go into greater depth. A total honour to have you analyse our data. Thank you so much! I must say I feel very relieved that the results you found with the beta-binomial substantively match those we found in the original study. Best of all, you have given me a roadmap for how to model the data in a study we are currently setting up with essentially the exact same outcome but with a different drug. I feel like I should send you a bottle of scotch to thank you.

avehtari · February 27, 2024, 7:46am

You can thank me in acknowledgments and send me a link to the paper. Feel free to ask questions about the other parts of the workflow, too

llewmills · February 27, 2024, 10:12am

Deal, though I must say I came out ahead in this arrangement. I don’t anticipate a paper coming out till 2026-27; addiction trials take a LONG time. Thank you so much again. I haven’t had a chance to absorb everything yet but I will follow up with questions as they inevitably occur to me.

llewmills · February 28, 2024, 12:37am

I noticed a typo in your notebook. In the ‘Comparison of discrete and continuous models section’ is a sentence that reads ‘The normal model predictions are not constrained between 0 and 28 (or −0.5 and 28.5), and thus the probabilities for cu ∈(0,1,2,…,28) do not necessarily sum to 1’. Now when I copied the text into this reply the numbers appeared but in the html document itself they are invisible. Just thought you might like to know. :)

avehtari · February 28, 2024, 8:41am

This shows fine in my browser

I did notice that it can take up to 10s that they come up. The html is self-contained, and with all the figures quite big, and mathjax can then be slow, too

Topic		Replies	Views
Loo for a subset of the data General loo , brms	4	621	August 4, 2021
Loo to compare ordered vs continuous response models brms loo	2	501	November 14, 2018
Interpreting varying and unstable model comparison results brms loo , interpret-results	2	722	August 9, 2018
Model comparison- big data Modeling loo	5	753	July 24, 2018
Comparing loo for two complex multivariate models with large data brms loo	9	1050	October 1, 2020

Comparing models with different functional forms of the same y variable using loo_compare

Related Topics