Comparing loo of discretised versions of hurdle lognormal and normal models

Modeling
1

I am trying to evaluate the relative performance of a hurdle lognormal model against a standard linear/gaussian model, conducted on the same data, a right-skewed integer 0-100 scale with a lot of 0’s and a normal-ish distribution for the non-zero options.

In a reply to this post @avehtari suggests that to compare the performance of a hurdle lognormal model to a gaussian model it is necessary to “do the model comparison in discrete space by completely discretizing both models” . He suggests that “for integer valued targets the discretization is super easy” and suggested working through (I assume Section 4 of) the Nabiximols case study.

Unsurprisingly I have not found it as easy as @avehtari, not least because I don’t know what doing the model ‘in discrete space’ means. I need some assistance. I will use the Nabiximols data as an example for illustration purposes.

Here are the data, comparing days of cannabis use in the previous 28-day period in two groups - the placebo and nabiximols groups - at 0, 4, 8, and 12 weeks during a 12-week clinical trial. Hurdle model is probably not the best model for these data but I need to know how to do the discretisation procedure on a hurdle lognormal for my 0-100 outcome variable.

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_df <- data.frame(id, group, week, cu, set)

# remove nas
cu_df <- cu_df |>
           drop_na(cu)

And here are the models and comparing them with loo_compare()

# gaussian model
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,
                            moment_match = TRUE)

loo_hurdle <- loo(fit_normal)

# hurdle model
fit_hurdle <- brm(formula = bf(cu ~ group*week + (1 | id),
                               hu ~ group*week,
                               decomp = "QR"),
                  data = cu_df,
                  family = hurdle_lognormal(),
                  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_hurdle <- add_criterion(fit_hurdle, 
                            criterion = "loo", 
                            save_psis = TRUE,
                            moment_match = TRUE)

loo_hurdle <- loo(fit_hurdle)

loo_compare(loo_normal, loo_hurdle)

# output)
      model elpd_diff se_diff p_worse diag_diff      diag_elpd
 fit_normal       0.0     0.0      NA                         
 fit_hurdle    -114.2    24.4    1.00           1 k_psis > 0.7

The normal model appears superior, with |elpd_diff|/se_diff > 2

Now let’s compare discretised version of each model.

i <- 2
llHurdle <- log_lik(object = fit_hurdle, 
                    newdata = cu_df[i, ] |> 
                 select(!cu) |> 
                   expand_grid(cu = seq(from = 0, 
                                        to = 28, 
                                        by = 1)))

loo_hurdle_ll <- loo(llHurdle)

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

loo_normal_ll <- loo(llNormal)

Now when I try to compare the log-likelihoods

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                         
 model1     -11.3    10.6    0.86   N < 100 2 k_psis > 0.7

The normal is still better but no longer notably so. I just need some advice concerning whether I have done this right. For example I don’t know if running the log_lik() function through loo_compare() via loo() is the way to go.

Guidance much appreciated.

2

In your earlier post you did not say what is your target variable and since you were including hurdle lognormal, I assumed you have continuous target in interval [0, infty), which is the difficult case.

If the target is integer count of days, then you are in the super easy case where you can compare the models as discussed in that case study. That case study shows that for integer counts the discretization intervals have all width 1, and then density times 1 is the probability. Thus you can use loo_compare directly.

I don’t understand how did you end up having 100 as the highest value, if the target is days of cannabis use in the previous 28-day period

And why not use the models that are in the Nabiximols case study for cu?

Maybe I’m still missing some key information

3

Thank you for replying @avehtari. The outcome I am really interested in is a 0-100 integer variable with a lot of zeros and a normal-looking distribution across the non-zero values of the outcome. It seemed like the sort of distribution that might suit a hurdle lognormal model. I used the Nabiximols data as an example - because it is also an integer outcome variable - to try and understand how the discretisation procedure works for comparing hurdle lognormal to gaussian via loo_compare(), with my intention to apply what I learned to the 0-100 integer variable.

My apologies for not being clear.

So are you saying that I can compare a hurdle lognormal directly to a gaussian if the outcome variable in both cases is a non-negative integer with an upper limit, in the same way that you can compare a beta-binomial to a gaussian, based on the reasoning outlined in the nabiximols case study?

Now that I see your reply I am wondering if a betabinomial might work better even than a hurdle lognormal.

4

If you have a lot of zeros, maybe zero-inflated beta-binomial?

5

I see, so hurdle lognormal is more for distributions where there appears to be a different process for the zeros vs non-zero, and for the non-zeros. And for continuous unbounded data. Thank you again. I will investigate the zi betabinomial.

6

I think I will post a better example on a different post, one with a continuous, non-negative outcome with many zeroes, where there is a need for discretisation in the non-logistic portion of the hurdle as well as the gaussian.