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). The data are used to illustrate various workflow aspects including comparison of models, but the analysis here do not match exactly the analysis in the paper or in the follow-up paper by Lintzeris et al. (2020), and further improvements in the analysis could be made by using additional data not included in the data used here.
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.
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 3 observations with a pareto_k > 0.7 in model 'fit_normal'. We
recommend 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 112 observations with a pareto_k > 0.7 in model 'fit_binomial'.
We recommend 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 can 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.9533272
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))
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 1 observations with a pareto_k > 0.7 in model
'fit_normal_scaled'. We recommend 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
Let’s get back to the first 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.
Histograms of posterior predictive replicates show that the binomial model predicts lower probability for both 0 and 28 than what is the observed proportion.
This effect can be seen also in PIT-ECDF calibration check plot (Säilynoja, Bürkner, and Vehtari 2022) which shows that the predictive distribution is too narrow.
pp_check(fit_binomial, type="pit_ecdf")
Figure 3
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).
The normal model posterior predictive replicates look even more different than the observed data, including predicting outcomes less than 0 and larger than 28, but the higher probability of 0 and 28 makes the difference that the normal model wins in the LOO comparison.
This effect can be seen also in PIT-ECDF calibration check plot which shows that the predictive distribution is too wide, but still better than the binomial model one.
pp_check(fit_normal, type="pit_ecdf")
Figure 5
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 14 observations with a pareto_k > 0.7 in model
'fit_betabinomial'. We recommend to set 'moment_match = TRUE' in order to
perform moment matching for problematic observations.
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
The posterior predictive replicates from beta-binomial look very much like the original data.
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")
Figure 8
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)
Figure 10
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)
Figure 11
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).
There are prior-data conflicts 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 3 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 3 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 which 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)
Figure 13
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)
Figure 14
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.
For a new individual the posterior predictive distribution indicates 90% possibility of smaller cu with Nabiximols than with placebo. I have used dot plots with 100 dots, as dot plots are better than kernel density estimates and histograms for showing spikes (as here at 0) and make it easy to estimate tail probabilities by counting the dots.
For comparison I plot 1) the ggplot2 default KDE with stat_density, 2) ggdist KDE with stat_slabinterval, and 3) ggplot default histogram with stat_histogram for the week 12 difference. The ggplot2 KDE oversmooths a lot, the ggdist KDE and ggplot2 histogram smooth less, but it is not as easy to estimate the tail probabilities as with dots plot.
To get more accuracy we can remove the aleatoric uncertainty and focus on expected effect. 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 (using re_formula=NA excludes the uncertainty in the unknown id specific intercept and using _epred gives just the expectation excluding the uncertainty in beta-binomial distribution).
Finally we check how much there is difference in the conclusions, if we had used the continuos normal model. We include also models that are closer to the model reported by Lintzeris et al. (2019), which used sum of days of cannabis use across the 12-week trial (the follow-up paper by Lintzeris et al. (2020) considered several 4 week intervals). We build normal and beta-binomial models which match the model in the paper, expect that we don’t include site factor as this was not available for us. #| results: hide #| cache: false
fit_normal2b <-brm(formula = cu ~ group*week + cu_baseline + (1| id),data = cu_df_b,gaussian(),prior =c(prior(normal(0, 3), class = Intercept),prior(normal(0, 3), class = b)),save_pars =save_pars(all=TRUE),seed =1234,refresh =0)
Running MCMC with 4 sequential chains...
Chain 1 finished in 0.5 seconds.
Chain 2 finished in 0.5 seconds.
Chain 3 finished in 0.5 seconds.
Chain 4 finished in 0.5 seconds.
All 4 chains finished successfully.
Mean chain execution time: 0.5 seconds.
Total execution time: 2.3 seconds.
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_normal2b'. 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.
The beta-binomial model beats the normal model big time.
loo_compare(fit_normal2b, fit_betabinomial2b)
Warning: Not all models have the same y variable. ('yhash' attributes do not
match)
We’re not able to compare the models predicting total count in 12 weeks and models predicting counts in three 4-week period using elpd_loo as the targets are different.
We plot here the difference either for the last 4 week period or for the all weeks depending on the model.
The normal models underestimate the magnitude of the change and are overconfident having much narrower posteriors. When we compare posteriors for effect in weeks 9-12 given the models which included week as a covariate, the conclusion about benefit of Nabiximols is the same with normal and beta-binomial model despite that the normal model underestimates the effect. When we compare posteriors for effect in weeks 1-12 given the models which did not include week as a covariate, the conclusion on benefit of Nabiximols is clearly different. In this case, we know based on LOO-CV comparisons and posterior predictive checking that the beta-binomial has much better performance and matches the data much better, and we can safely drop the normal models from consideration. The conclusions using the models with week as a covariate or all counts summed together are similar, although focusing in the last four week period gives sharper posterior which is plausible as the effect seems to be increasing in time.
The results here are illustrating that the choice of data model can matter for the conclusions, but also the best models shown here could be even further improved by using the additional data used by Lintzeris et al. (2020).
15 Predictive performance comparison
Above, posterior predictive checking, LOO predictive checking, and LOO-ELPD comparisons were important for model checking and comparison, and only after we trusted that there is no significant discrepancy between the model predictions and the data, we did look at the treatment effect. We used LOO model comparison to assess the the more elaborate models did provide better predictions, but we did not use it above for assessing the predictive performance of using the treatment effect in the model.
Let’s now build the mode without treatment group variable and check whether there is significant difference in predictive performance.
Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
Warning: Found 4 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 4 times to compute the ELPDs for the problematic
observations directly.
We compare the models with and without group variable:
There is a negligible difference in predictive performance for predicting cannabis use in a 4-week period for a new person. 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 in tails.
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
By dropping out the aleatoric part from the LOO-CV predictions we do see clear difference in the predictive performance if the treatment group variable is dropped, but still the randomness in the actual observations makes the probability of the difference to be further away from 1 than we focused in the expected counterfactual predictions above.
We don’t usually use Bayes factors for many reasons, but for completness we used bridgesampling(Gronau, Singmann, and Wagenmakers 2020) to get estimated Bayes factor.
Bridge sampling estimate of the log marginal likelihood: -606.4643
Estimate obtained in 10 iteration(s) via method "normal".
bridgesampling::bf(br2b, br3b)
Estimated Bayes factor in favor of br2b over br3b: 0.50562
Like LOO-CV, Bayes factor does not see difference between models with or without treatment group variable. Bayes factor is a ratio of marginal likelihoods, and marginal likelihoods can be formulated as sequential log score predictive performance quantity, which explains why BF is also weak to rank models in case of high aleatoric uncertainty.
References
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.
Lintzeris, Nicholas, Llewellyn Mills, Adrian Dunlop, Jan Copeland, Iain Mcgregor, Raimondo Bruno, Adrienne Kirby, et al. 2020. “Cannabis Use in Patients 3 Months After Ceasing Nabiximols for the Treatment of Cannabis Dependence: Results from a Placebo-Controlled Randomised Trial.”Drug and Alcohol Dependence 215 (108220).
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.
