Discretizing a hurdle model in order to compare it to a gaussian model: Part 2

Modeling
1

In this post a few days ago I attempted to get some advice on the procedure for discretizing a hurdle model and a gaussian model of the same outcome data so as to compare the two (recommended by @avehtari in his reply here). Professor Vehtari recommended using the process outlined in Section 4 of his Nabiximols case study.

I attempted to conduct the procedure on the data included in the case study, however this was not such a great candidate for a hurdle model as the outcome variable was an integer count variable, which, as described in the Section 4 of the case study, would be better modelled as a betabinomial. One of the great properties of the binomial family, when used to model non-negative integer data is that, because the width of each interval on the continuous scale is always 1, the probability value from the binomial model is the same as the density value from the gaussian model, and thus you can compare the betabinomial model to the gaussian directly without any need for extra discretisation steps. However this is not true of a hurdle model on continuous data (more below).

I have some new data from a clinical trial I am analysing and this time the outcome variable is an ideal candidate for a hurdle model. It is continuous, highly right-skewed, non-negative and includes many zeros. I include the data here

d.RData (5.7 KB)

The outcome (named outcome) is measured a maximum of six times through the trial at weeks 0 (baseline), 3, 6, 9, 12, and at a follow-up interview 12 weeks after the medication has been stopped (some participants drop out early and hence have less than six measures. The measurement occasion information is contained in a factorial variable called weekFac. There are two types of medication (contained in variable group): placebo and 'exp. Participant ID is contained in the id` variable. I want to examine group differences at each measurement occasion

Code to prepare environment, load the dataset, and remove missing values.

# load packages
library(tidyr)
library(matrixStats)
library(dplyr)
library(tibble)
library(modelr)
library(loo)
library(brms)
library(rstan)
library(posterior)
options(posterior.num_args = list(digits = 2))
library(pillar)
options(pillar.negative = FALSE)
library(ggplot2)
library(bayesplot)
theme_set(bayesplot::theme_default(base_family = "sans", base_size = 14))
library(tidybayes)
library(patchwork)
library(tinytable)
options(
  tinytable_format_num_fmt = "significant_cell",
  tinytable_format_digits = 2,
  tinytable_tt_digits = 2
)

# load data
load(file = "d.RData")

# remove nas
d <- d |>
        drop_na(outcome)

The outcome is on a range from 0-13 and looks like this

Rplot
Rplot722×518 5.07 KB

Now I run the two models, hurdle lognormal and gaussian, and run each through loo.

# gaussian model
fit_normal <- brm(formula = outcome ~ group*weekFac + (1 | id),
                  data = d,
                  family = gaussian(),
                  save_pars = save_pars(all = TRUE),
                  seed = 1234,
                  refresh = 0)

fit_normal <- add_criterion(fit_normal,
                            criterion = "loo",
                            save_psis = TRUE,
                            moment_match = TRUE)


loo_normal <- loo(fit_normal)

# hurdle model
fit_hurdle <- brm(formula = bf(outcome ~ group*weekFac + (1 | id),
                               hu ~ group*weekFac,
                               decomp = "QR"),
                  data = d,
                  family = hurdle_lognormal(),
                  save_pars = save_pars(all = TRUE),
                  seed = 1234,
                  refresh = 0)

fit_hurdle <- add_criterion(fit_hurdle,
                            criterion = "loo",
                            save_psis = TRUE,
                            moment_match = TRUE)

loo_hurdle <- loo(fit_hurdle)

Now I am going to ignore all the model refinement and diagnostics for the purposes of this post. What I need advice on is the process of discretisation.

And when I compare them

loo_compare(loo_normal, loo_hurdle)

# output
#      model elpd_diff se_diff p_worse diag_diff       diag_elpd
# fit_hurdle       0.0     0.0      NA           10 k_psis > 0.7
# fit_normal    -426.3   149.8    1.00           12 k_psis > 0.7

The hurdle is superior. HOWEVER, from what I gather, as-is, this loo comparison is not valid because the hurdle lognormal contains both a discrete logit model for ‘0 vs everything else’ and a continuous portion modelling the non-zero values. So Professor Vehtari advised doing “the model comparison in discrete space by completely discretizing both models” and suggested using the Nabiximols case study as a guide.

Now, I have gone through Section 4 of the case study, which I think is the relevant section. And show my attempts at the discretization below.

The closest I can get is to use the log_lik() and expand_grid() functions.

llNormal <- log_lik(object = fit_normal, 
                    newdata = d[i, ] |> 
                                select(!outcome) |> 
                                  expand_grid(outcome = seq(from = 0, 
                                                            to = 13, 
                                                            by = 1)))

llHurdle <- log_lik(object = fit_hurdle, 
                    newdata = d[i, ] |> 
                                select(!outcome) |> 
                                  expand_grid(outcome = seq(from = 0, 
                                                            to = 13, 
                                                            by = 1)))

It turns out you can pass the output of the log_lik function directly into the loo() function

loo_hurdle_ll <- loo(llHurdle)
loo_normal_ll <- loo(llNormal)

And compare the two via loo_compare()

loo_compare(loo_normal_ll, loo_hurdle_ll)

# output
#  model elpd_diff se_diff p_worse diag_diff       diag_elpd
# model2       0.0     0.0      NA            1 k_psis > 0.7
# model1   -1071.9   303.1    1.00   N < 100 11 k_psis > 0.7

Now it looks like once again the now-discretized (?) hurdle model is superior to the gaussian. But I am fairly sure I am doing something wrong. In the case study Professor Vehtari uses the colLogSumExps() function from the matrixStats package to integrate across intervals, but I can’t figure out how to then use that output to compare models.

So I am flying blind basically, with little idea how to proceed. I have an output, which is better than an error message I suppose, but I have no idea if the output can be trusted. I need advice from people who know (much) more than me about whether what I did is right and, if not, how to do it properly.

2

It took some time to carefully think this through

  • If there clearly are many zeros, we know hurdle model will have an advantage even before doing elpd_loo comparison. If the non-zeros are truly continuous valued they all should be unique (probability of observing exactly the same value twice drawn from a continuous distribution is 0) and having more than 0 already states that a continuous model is wrong model.
  • Comparing continuous, discrete, and continuous-discrete models is easier when the data are truly discrete
  • If the data are continuous or continuous-discrete, the results can be sensitive to specific discretization
  • It’s easiest to start with discretizing to intervals having (mostly width 1). Two options:
    a) Both models share the same discretization intervals [0.5,1.5), [1.5, 2.5), etc. (it doesn’t matter which interval end point is open) And for all continuous valued observations above 0.5 we can use the same logic as in nabiximols case study, that is density times interval width 1 is probability. For smooth densities we don’t need to care that continuous values are somewhere in the interval and those continuous values can be used as they are. For normal model we can also discretize using interval [-0.5, 0.5). For the hurdle model for 0’s we get probabilities directly, and then values in interval (0, 0.5) get probabilities by multiplying density by 2 (as the interval has width 0.5).
    b) Discretize with intervals (-Infty,0],(0,1],(1,2], etc. For positive range intervals the same works as in a. For hurdle model probability of 0 comes from the model directly. For normal model, use CDF of normal as the probability of 0.
  • The option a) is more sensitive to making the interval shorter and option b) is making a bigger mismatch between the fitted normal model and how it is evaluated.
  • We haven’t yet discussed what is the data collection process for 0’s. If the 0’s are due to some detection threshold or other rounding down due to accuracy of the measurement, it would make sense to always model it as left censored data. Hurdle model is one way without explicitly stating the detection threshold. For normal model the above mentioned cdf approach would match left-censoring model, or even better would could use log-normal and then 0’s would be presented in the data as left censord observation with known or unknown censoring point.
  • I guess I’m providing too many options for you, but then you did not yet provide enough information about where the data comes from. We can have also a zoom call to discuss the data collection process.
3

Thanks so much Aki. I can share that the outcome is average grams per day of a substance. And pretty continuous. If I table the outcome it looks like this (not full output)

table(d$outcome)

                  0 0.00119047619047619 0.00214285714285714 0.00238095238095238 
                 81                   1                   1                   1 
             0.0025 0.00297619047619048 0.00428571428571429 0.00476190476190476 
                  1                   1                   1                   7 
0.00530120481927711 0.00595238095238095 0.00952380952380952  0.0104761904761905 
                  1                   1                   3                   1 
 0.0107142857142857  0.0119047619047619              0.0125  0.0128205128205128 
                  1                   2                   2                   1 
 0.0133333333333333  0.0142857142857143  0.0157894736842105  0.0166666666666667 
                  1                   1                   1                   2 
 0.0177777777777778  0.0178571428571429  0.0189285714285714   0.019047619047619

So there are some values where there are more than one instance, but generally the modal number of instances is 1.

Given that there is a hard lower limit of 0 but not theoretical upper limits I guess the intervals would be right censored? I don’t really know though.

The thing that I still don’t know is the process of discretisation: the steps to get to the point where I can compare the hurdle to the gaussian.

For example, I don’t know if we discretise the outcome before it goes into the model i.e. if we bin any values between 1.5 and 2.4999999999 into a categorical “1.5to2.49” value, between 2.5 and 3.4999999999 into a categorical “2.5to3.49” value and then regress that new outcome into a on the predictors. But I can’t see how that would make sense as the approach because then the models would stop being gaussian and hurdle.

So I assume the discretisation happens via some kind of sampling from the posterior, something like this from the Nabiximols case study

ll3 <- log_lik(object = fit_normal, 
               newdata = d[i, ] |> 
                           select(!outcome) |> 
                             expand_grid(outcome = seq(-0.5, 13.5, by = 1)))

# now to get total probability for each 1-unit-wide interval 
p3 <- exp( colLogSumExps(lx = ll3) - log(S) )

But this doesn’t look right either as expand_grid(outcome = seq(-0.5, 13.5, by = 1) does not produce intervals, just chains of estimates at each interval boundary. I imagine each interval would have to go up to infinitessimally less than lower limit of the next interval? So I do not know how I would specify that in expand grid? Is the process to sample chains of point estimates for each of a series of many very small increments within each interval and then integrate across them? And I am also unaware of how to integrate over each interval once I obtain all these point estimates. exp( colLogSumExps(lx = ll3) - log(S) ) in the code above looks as if it gives probability at a single points estimate rather than integrating across an array of point estimates as I would expect for an interval.

Basically, as you can see, I am unsure of next steps to get to the point where we can compare the two discretised models. But I am very curious to find out how. There are several outcomes other than this one from the current experiment that would suit a hurdle, so it would be great to figure out the practical steps to achieving the comparison.