Extracting predictions for each outcome in multivariate ordinal model

General
1

I am trying to extract posterior predictions from a brms multivariate ordinal regression object but I can’t figure out how to extract the predictions for each outcome separately

Here is the model

bform_ord_mvr_ppq <- bf(mvbind(psych, phys, qol) ~ ats_baseline*yearsFromStart + (1|p|encID))

fit_ordinal_mvr <- brm(bform_ord_mvr_ppq,
                       family = cumulative("probit"),
                       data = mvrOrd,
                       seed = 1,
                       refresh = 0,
                       chains = 4,
                       warmup = 1e3,
                       iter = 3e3)

Now when I try to extract the posterior predictions for each of five yearsFromStart values at each of five ats_baselinelevels.

library(modelr)
library(tidybayes)

tibble(encID=688) %>%
    data_grid(yearsFromStart = seq(0,1,0.25),
              ats_baseline=c(0,7,14,21,28),
              encID) %>%
      add_epred_draws(object = final_model,
                      allow_new_levels = TRUE) %>%
        arrange(yearsFromStart, ats_baseline, .draw, .category) -> epredOrdDraws_5

The object looks like this (I left out the .chainand .iterationcolumns so it would fit on the page; they are all NA anyway)

  head(epredOrdDraws_5 %>% 
         select(yearsFromStart,
                ats_baseline,
                encID,
                .row,
                .draw,
                .category,
                .epred), 
       n = 20)


# A tibble: 20 × 7
# Groups:   yearsFromStart, ats_baseline, encID, .row, .category [7]
   yearsFromStart ats_baseline encID  .row .draw .category  .epred
            <dbl>        <dbl> <dbl> <int> <int> <fct>       <dbl>
 1              0            0   688     1     1 1         0.00101
 2              0            0   688     1     1 1         0.00114
 3              0            0   688     1     1 1         0.00137
 4              0            0   688     1     1 2         0.00261
 5              0            0   688     1     1 2         0.00119

Now when I do the gaussian version of this multivariate model, using the same outcome variables, the .categorycolumn has only three values: psych, phys, and qol and I could subset the object into three datasets to graph predictions for each of the three outcomes. But with this model the .categorycolumn has eleven values, one for each of the 11 breakpoints in the cumulative link model.

So my question is, does anyone know a way to subset the predictions into the three outcomes, similar to the gaussian model?

2

You can add resp = "psych" in add_epred_draws.

2 Likes