You can’t compare densities and probabilities directly. Thus you can’t compare model given continuous and discrete observation models, unless you compute probabilities in intervals from the continuous model (also known as discretising continuous model).
The answer is complete in theory, but doesn’t tell how to do the discretization in practice. The first part of this notebook shows how easy that discretization can be in a special case of counts.
Data comes from a study [Nabiximols for the Treatment of Cannabis Dependence: A Randomized Clinical Trial] by Lintzeris et al. (2019), and was posted in that Discourse thread by Mills (one of the co-authors). The data includes 128 participants (id) in two groups (group = Placebo or Nabiximols).
Participants received 12-week treatment involving weekly clinical reviews, structured counseling, and flexible medication doses—up to 32 sprays daily (tetrahydrocannabinol, 86.4 mg, and cannabidiol, 80 mg), dispensed weekly.
Primary outcome was self-reported number of days using illicit cannabis during the 12-week period.
The number of cannabis used days (cu) was for previous 28 days (set) asked after 0, 4, 8, and 12 weeks (week as a factor).
Mills provided two brms(Bürkner 2017) models with specified priors. The first one is a normal regression model with varying intercept for each participant (id).
fit_normal <-brm(formula = cu ~ group*week + (1| id),data = cu_df,family =gaussian(),prior =c(prior(normal(14, 1.5), class = Intercept),prior(normal(0, 11), class = b),prior(cauchy(1,2), class = sd)),save_pars =save_pars(all=TRUE),seed =1234,refresh =0)fit_normal <-add_criterion(fit_normal, criterion="loo", save_psis=TRUE)
Warning: Found 6 observations with a pareto_k > 0.7 in model 'fit_normal'. It
is recommended to set 'moment_match = TRUE' in order to perform moment matching
for problematic observations.
The second provided models is binomial model with the number of trials being 28 for each outcome (cu)
fit_binomial <-brm(formula = cu |trials(set) ~ group*week + (1| id),data = cu_df,binomial(link = logit),prior =c(prior(normal(0, 1.5), class = Intercept),prior(normal(0, 1), class = b),prior(cauchy(0,2), class = sd)),save_pars =save_pars(all=TRUE),seed =1234,refresh =0)fit_binomial <-add_criterion(fit_binomial, criterion="loo", save_psis=TRUE)
Warning: Found 110 observations with a pareto_k > 0.7 in model 'fit_binomial'.
It is recommended to set 'moment_match = TRUE' in order to perform moment
matching for problematic observations.
At this point I ignore the fact that there many high Pareto-k values indicating unreliable PSIS-LOO estimates (Vehtari et al. 2022), and first discuss comparison of predictive probabilities and predictive probability densities.
I pick the id=2 and compute the posterior predictive probability of each possible outcome (0,1,2,\ldots,28). At this point we don’t need to consider LOO-CV issues, and can focus on comparison of posterior predictive probabilities. We can use log_lik() function to compute the log predictive probabilities. In the following plot the area of each bar is equal to the probability of corresponding count outcome and these probabilities sum to 1. They grey bar illustrates the posterior predictive probability of cu=12.
In case of normal model, the observation model is continuous, and we compute posterior predictive densities. I use log_lik() to compute log predictive densities. I compute the predictive density of cu=12 and predictive density in fine grid for other values (although the observation model is continuous, in practice with computers, we can evaluate the density only in finite number of points). The following plot shows the predictive density and grey line at cu=12 with end of line corresponding to density at cu=12.
We can’t compare probabilities and densities directly, but we can discretize the density to get probabilities. As the outcomes are integers (0,1,2,\ldots,28), we can compute probabilities for intervals ((-0.5, 0.5), (0.5,1.5),
(1.5,2.5),\ldots,(27.5,28.5)). The following plot shows the vertical lines for the edges of these intervals.
Warning: Removed 190 rows containing missing values or values outside the scale range
(`geom_line()`).
We can then integrate the density over the interval. For example, integrate predictive density from 11.5 to 12.5 to get a probability that 11.5 < cu < 12.5 (it doesn’t matter here if the interval ends are open or closed). To simplify the computation, we approximate the density as piecewise constant function shown in the next plot.
Warning: Removed 190 rows containing missing values or values outside the scale range
(`geom_line()`).
Now the probability of each interval is approximated by the height times the width of a bar. The height is the density in the middle of the interval and width of the bar is 1, and thus the probability value is the same as the density value! In this case this is simple as the counts are integers and the distance between counts is 1. The following plot shows the probability of cu=12.
Warning: Removed 190 rows containing missing values or values outside the scale range
(`geom_line()`).
Thus, in this case, the LOO-ELPD comparison is valid.
The normal model predictions are not constrained between 0 and 28 (or -0.5 and 28.5), and thus the probabilities for cu\in (0,1,2,\ldots,28) do not necessarily sum to 1.
sum(p3)
[1] 0.9549362
We could switch to truncated normal, but for this data, it would not work well, and as we have better options, I skip illustration of that.
Before continuing with better models, I illustrate the case where the integration interval is not 1, and the comparison would require more work. It is common with continuous targets to normalize the target to have zero mean and unit standard deviation.
cu_scaled_df <- cu_df |>mutate(cu_scaled = (cu -mean(cu)) /sd(cu))
We fit the normal model with scaled target.
fit_normal_scaled <-brm(formula = cu_scaled ~ group*week + (1| id),data = cu_scaled_df,family =gaussian(),prior =c(prior(normal(0, 1.5), class = Intercept),prior(normal(0, 11), class = b),prior(cauchy(1,2), class = sd)),save_pars =save_pars(all=TRUE),seed =1234,refresh =0)fit_normal_scaled <-add_criterion(fit_normal_scaled, criterion="loo")
Warning: Found 5 observations with a pareto_k > 0.7 in model
'fit_normal_scaled'. It is recommended to set 'moment_match = TRUE' in order to
perform moment matching for problematic observations.
Now the densities for the scaled target are much higher, and the scaled model seems to be much better.
loo_compare(fit_normal, fit_normal_scaled)
Warning: Not all models have the same y variable. ('yhash' attributes do not
match)
However to make a fair comparison, we need to take into account the scaling we used. To get the probabilities for integer counts, the densities need to be multiplied by the discretization binwidth 0.093. Correspondingly in log scale we add \log(0.093) to each log-density, and the total adjustment is
nrow(cu_df)*log(sd(cu_df$cu))
[1] 914.2574
Adding this to eldp_loo of the scaled outcome model, makes the two models to have almost the same elpd_loo (the small difference comes from not scaling the priors).
5 Model checking
Back to the original models. The normal model was better than the binomial model in the LOO comparison, but that doesn’t mean it is a good model. Before trying to solve the high Pareto-k issues in LOO computation, we can use posterior predictive checking.
We see that binomial model has too many PIT values near 0 and 1, which indicates the posterior predictive intervals are too narrow, and the model is not appropriately handling the overdispersion compared to binomial model, even with included varying intercept term (1 | id).
How about the normal model?
pp_check(fit_normal, type="pit_ecdf")
Not bad. There are fewer PIT values near 0 and 1 than expected, but this can be because of double use of data in posterior predictive checking. We can get a more accurate PIT-ECDF plot by using LOO predictive distributions (even if we ignore the few high Pareto-k warnings).
At the moment fixed LOO-PIT for discrete outcomes is not yet in bayesplot, and I’m defining the necessary functions here.
This looks quite good, which is a bit surprising considering the prediction are not constrained between 0 and 28.
6 Model extension
Binomial model doesn’t have overdispersion term, but the normal model does. How about using beta-binomial model which is an overdispersed version of the binomial model.
fit_betabinomial <-brm(formula = cu |trials(set) ~ group*week + (1| id),data = cu_df,beta_binomial(),prior =c(prior(normal(0, 1.5), class = Intercept),prior(normal(0, 1), class = b),prior(cauchy(0,2), class = sd)),save_pars =save_pars(all=TRUE),seed =1234,refresh =0)fit_betabinomial <-add_criterion(fit_betabinomial, criterion="loo", save_psis=TRUE)
Warning: Found 18 observations with a pareto_k > 0.7 in model
'fit_betabinomial'. It is recommended to set 'moment_match = TRUE' in order to
perform moment matching for problematic observations.
Warning: Not all models have the same y variable. ('yhash' attributes do not
match)
The beta-binomial model beats the normal model big time. The difference is so big that it is unlikely that fixing Pareto-k warnings matter.
7 Improved LOO computation
Since in case of high Pareto-k warnings, the estimate is usually optimistic, it is sufficient if we use moment matching (Paananen et al. 2021) for beta-binomial model.
Truncated normal might get closer to beta-binomial than normal, but as I said it had difficulties here, and since beta-binomial is proper discrete observation model, I continue with that.
8 Model checking
PIT-ECDF calibration plot looks quite good, but there is a slight S-shape that might be explained by double use of data in posterior predictive checking.
pp_check(fit_betabinomial, type="pit_ecdf")
LOO-PIT-ECDF calibration looks very nice and also better than for the normal model.
The data included many counts of 0 and 28, and we can further check whether we might need to include zero-inflation or 28-inflation model component. As illustrated in Roaches case study PIT’s may sometimes be weak to detect zero-inflation, but reliability diagram (Dimitriadis, Gneiting, and Jordan 2021) can focus there in more detail.
Calibration check with reliability diagram for 0 vs others
th<-0rd=reliabilitydiag(EMOS =pmin(E_loo(0+(posterior_predict(fit_betabinomial)>th),loo(fit_betabinomial)$psis_object)$value,1),y =as.numeric(cu_df$cu>th))autoplot(rd)+labs(x="Predicted probability of non-zero",y="Conditional event probabilities")+ bayesplot::theme_default(base_family ="sans", base_size=14)
Calibration check with reliability diagram for 28 vs others
th<-27rd=reliabilitydiag(EMOS =pmin(E_loo(0+(posterior_predict(fit_betabinomial)>th),loo(fit_betabinomial)$psis_object)$value,1),y =as.numeric(cu_df$cu>th))autoplot(rd)+labs(x="Predicted probability of 28",y="Conditional event probabilities")+ bayesplot::theme_default(base_family ="sans", base_size=16)
Although the red line did not completely stay within blue envelope, these look good.
9 Prior sensitivity analysis
The priors provided by Mills seems a bit narrow, so I do prior-likelihood sensitivity analysis using powerscaling approach (Kallioinen et al. 2023).
We have prior-data conflict for all global parameters. It is possible that the priors had been chosen based on substantial prior information and then the posterior should be influenced by the prior, too, and the conflict is not necessarily a problem.
10 Alternative priors
However, I decided to try slightly wider priors, especially as the data seem to be quite informative
fit_betabinomial2 <-brm(formula = cu |trials(set) ~ group*week + (1| id),data = cu_df,beta_binomial(),prior =c(prior(normal(0, 3), class = Intercept),prior(normal(0, 3), class = b),prior(normal(0, 3), class = sd)),save_pars =save_pars(all=TRUE),seed =1234,refresh =0)
All the above models included the baseline week=0 in the interaction term with group, which does not make sense as the group should not affect the baseline. I modify the data and models by moving the baseline week=0cus to be pre-treatment covariate.
Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
Warning: Found 5 observations with a pareto_k > 0.7 in model
'fit_betabinomial2b'. We recommend to set 'reloo = TRUE' in order to calculate
the ELPD without the assumption that these observations are negligible. This
will refit the model 5 times to compute the ELPDs for the problematic
observations directly.
To compare the new model against the previous ones, we exclude the pointwise elpds for week 0.
The quick fix to exclude week 0, does not change the yhash and we get a warning we can ignore. The new model is clearly better than the previous best model.
By comparing phi parameter posteriors, we see that the overdispersion of the beta.binomial in the new model is smaller (smaller phi means bigger overdispersion).
Calibration check with reliability diagrams for 0 vs others and 28 vs others look better than for the previous model.
th<-0rd=reliabilitydiag(EMOS =pmin(E_loo(0+(posterior_predict(fit_betabinomial2b)>th),loo(fit_betabinomial2b)$psis_object)$value,1),y =as.numeric(cu_df_b$cu>th))autoplot(rd)+labs(x="Predicted probability of non-zero",y="Conditional event probabilities")+ bayesplot::theme_default(base_family ="sans", base_size=16)
th<-27rd=reliabilitydiag(EMOS =pmin(E_loo(0+(posterior_predict(fit_betabinomial2b)>th),loo(fit_betabinomial2b)$psis_object)$value,1),y =as.numeric(cu_df_b$cu>th))autoplot(rd)+labs(x="Predicted probability of 28",y="Conditional event probabilities")+ bayesplot::theme_default(base_family ="sans", base_size=16)
13 Treatment effect
As the treatment is included in interaction term, looking at the univariate posterior marginal is not sufficient for checking whether the treatment helps.
We can examine the effect of the treatment by looking at the posterior predictive distribution and the expectation of the posterior predictive distribution for a new person (new id). We use tidybayes(Kay 2023) and ggdist(Kay 2024).
We start by examining the posterior predictive distribution. We predict cu for a new id=129 given cu_baseline=28 (median value). The distribution is wide as it included the aleatoric uncertainty from the unknown id specific intercept and beta-binomial distribution.
The difference between groups is the biggest in week 12. The next plot shows the comparison of predicted cu between the groups in different weeks.
cu_df_b |>data_grid(group, week, cu_baseline=28, id=129, set=28) |>add_predicted_draws(fit_betabinomial2b, allow_new_levels=TRUE) |>compare_levels(.prediction, by=group) |>ggplot(aes(x = .prediction)) +facet_grid(. ~ week, switch ="y", axes="all_x", labeller = label_both) +stat_dotsinterval(quantiles=100)+coord_cartesian(expand =FALSE, clip="off") +theme(strip.background =element_blank(), strip.placement ="outside") +labs(x="Difference in cu given placebo vs Nabiximols", y="") +geom_vline(xintercept=0, linetype="dashed") +theme(axis.line.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank(),axis.line.x=element_blank())
For a new individual the posterior predictive distribution indicates 90% possibility of smaller cu with Nabiximols than with placebo.
The following plots shows the expectation of the posterior predictive distribution for cu for a new id=129 given cu_baseline=28 (expectation as in average over many new individuals). The distribution is much narrower as it includes only the epistemic uncertainty, but not the aleatoric uncertainty (from the unknown id specific intercept and beta-binomial distribution).
Finally we can check how much there is difference in the conclusions, if we had used the continuos normal model or underdispersed binomial model.
The normal model underestimates the magnitude of the change and is overconfident having much narrower posterior. In this case, the conclusion that Nabiximols is beneficial is the same, but too narrow posterior might in some other case lead to an overoptimistic conclusion.
The binomial model is overconfident having much narrower posterior. In this case, the conclusion that Nabiximols is beneficial is the same, but too narrow posterior might in some other case lead to an overoptimistic conclusion.
Above, posterior predictive checking, LOO predictive checking, and LOO-ELPD comparisons were important for model checking and comparison, and only after I trusted that there is no significant discrepancy between the model predictions and the data, I did look at the treatment effect. I used LOO model comparison to assess the the more elaborate models did provide better predictions, but I did not use it above for assessing the predictive performance of using the treatment effect in the model.
Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
Warning: Found 3 observations with a pareto_k > 0.7 in model
'fit_betabinomial3b'. We recommend to set 'reloo = TRUE' in order to calculate
the ELPD without the assumption that these observations are negligible. This
will refit the model 3 times to compute the ELPDs for the problematic
observations directly.
And compare the models with and without group variable:
There is no difference in predictive performance. This is probably due to
1. the actual effect not being very far from 0,
2. the aleatoric uncertainty (modelled by the beta-binomial)
for a new 4-week period is big, that is the predictive
distribution is very wide and due to the constrained range
has also thick tails (actually U shape), which makes the log
score not to be sensitive to tail values as usually
I don’t usually use Bayes factors for many reasons, but for completness showing here that bridgesampling(Gronau, Singmann, and Wagenmakers 2020) estimated BF doesn’t either see difference with or without treatment models.
Bridge sampling estimate of the log marginal likelihood: -606.3994
Estimate obtained in 8 iteration(s) via method "normal".
bridgesampling::bf(br2b, br3b)
Estimated Bayes factor in favor of br2b over br3b: 0.41407
As the predictive distribution is wide with thick tails, we can also focus on comparing absolute error of using means of the predictive distributions as point estimates. This approach has the benefit that we can improve accuracy of the finite MCMC sample by dropping the sampling from the varying intercept population distribution and from the beta-binomial distribution (using posterior_epred(, re_formula=NA). E_loo is used to go from posterior predictive draws to the mean of leave-one-out predictive distribution.
Probability that the leave-one-out predictive absolute error is smaller with model 2b (with group variable) than with model 3b (without group variable) using normal approximation
Bürkner, Paul-Christian. 2017. “brms: An R Package for Bayesian Multilevel Models Using Stan.”Journal of Statistical Software 80 (1): 1–28. https://doi.org/10.18637/jss.v080.i01.
Dimitriadis, Timo, Tilmann Gneiting, and Alexander I. Jordan. 2021. “Stable Reliability Diagrams for Probabilistic Classifiers.”Proceedings of the National Academy of Sciences 118 (8): e2016191118.
Gronau, Quentin F., Henrik Singmann, and Eric-Jan Wagenmakers. 2020. “bridgesampling: An R Package for Estimating Normalizing Constants.”Journal of Statistical Software 92 (10): 1–29. https://doi.org/10.18637/jss.v092.i10.
Kallioinen, Noa, Topi Paananen, Paul-Christian Bürkner, and Aki Vehtari. 2023. “Detecting and Diagnosing Prior and Likelihood Sensitivity with Power-Scaling.”Statistics and Computing 34 (1): 57.
Lintzeris, Nicholas, Anjali Bhardwaj, Llewellyn Mills, Adrian Dunlop, Jan Copeland, Iain McGregor, Raimondo Bruno, et al. 2019. “Nabiximols for the Treatment of Cannabis Dependence: A Randomized Clinical Trial.”JAMA Internal Medicine 179 (9): 1242–53.
Paananen, Topi, Juho Piironen, Paul-Christian Bürkner, and Aki Vehtari. 2021. “Implicitly Adaptive Importance Sampling.”Statistics and Computing 31 (16).
Säilynoja, Teemu, Paul-Christian Bürkner, and Aki Vehtari. 2022. “Graphical Test for Discrete Uniformity and Its Applications in Goodness-of-Fit Evaluation and Multiple Sample Comparison.”Statistics and Computing 32: 1573–375.
Vehtari, Aki, Jonah Gabry, Mans Magnusson, Yuling Yao, Paul-Christian Bürkner, Topi Paananen, and Andrew Gelman. 2024. “Loo: Efficient Leave-One-Out Cross-Validation and WAIC for Bayesian Models.”https://mc-stan.org/loo/.
Vehtari, Aki, Andrew Gelman, and Jonah Gabry. 2017. “Practical Bayesian Model Evaluation Using Leave-One-Out Cross-Validation and WAIC.”Statistics and Computing 27: 1413–32.
Vehtari, Aki, Daniel Simpson, Andrew Gelman, Yuling Yao, and Jonah Gabry. 2022. “Pareto Smoothed Importance Sampling.”arXiv Preprint arXiv:1507.02646. https://arxiv.org/abs/1507.02646v6.
