Specifying Priors for Ordinal (adjacent category) Model

Interfaces brms
,
1

I am attempting to model some ordinal response data with a hierarchical model. I this study, 75 sound tokens were grouped into 5 broader categories with repeated presentation across several conditions. The conditions were nested within 2 groups of participants. The outcome was an ordinal response (1 of 9 pictures).

My brms model is specified as follows:

# Model formula
Mod_fmla <- bf(Outcome ~ group +
                           condition +
                           new_category +
                            (1|participant) +
                           (1|new_category/token) +
                           (1|condition/group))

I am specifying my priors as follows:

priors_acat <- c(
        set_prior("normal(0,1)",  
                  class = "b",
                  coef = "condition80"),
        
        set_prior("normal(0,1)",  
                  class = "b",
                  coef = "conditionCP"),
        
        set_prior("normal(0,1)",  
                  class = "b",
                  coef = "conditionNFC"),
        
        set_prior("normal(0,1)",  
                  class = "b",
                  coef = "groupHI"),
        
        set_prior("normal(0,1)",  
                  class = "b",
                  coef = "new_categoryunpleasant_low"),
        
        set_prior("normal(0,1)",  
                  class = "b",
                  coef = "new_categoryneutral"),
        
        set_prior("normal(0,1)",  
                  class = "b",
                  coef = "new_categorypleasant_low"),
        
        set_prior("normal(0,1)",  
                  class = "b",
                  coef = "new_categorypleasant_high"),
        
        set_prior("normal(0.1,1)",  
                  class = "Intercept"),
        
        set_prior("normal(0,1)",  
                  class = "Intercept",
                  coef = "1"),
        
        set_prior("normal(0,1)",  
                  class = "Intercept",
                  coef = "2"),
        
        set_prior("normal(0.4,1)",  
                  class = "Intercept",
                  coef = "3"),
        
        set_prior("normal(0.4,2)",  
                  class = "Intercept",
                  coef = "4"),
        
        set_prior("normal(5,1)",  
                  class = "Intercept",
                  coef = "5"),
        
        set_prior("normal(6,1)",  
                  class = "Intercept",
                  coef = "6"),
        
        set_prior("normal(1,4)",  
                  class = "Intercept",
                  coef = "7"),
        
        set_prior("normal(0,1)",  
                  class = "Intercept",
                  coef = "8"),
        
        set_prior("cauchy(0,1)",  
                  class = "sd",
                  group = "participant",
                  coef = "Intercept"),
        
        set_prior("cauchy(0,1)",  
                  class = "sd",
                  group = "condition:group",
                  coef = "Intercept"),
        
        set_prior("cauchy(0,1)",  
                  class = "sd",
                  group = "new_category:token",
                  coef = "Intercept")
        
)

In order to check priors, Iā€™m trying a prior predictive check on a ā€œprior-onlyā€ model:

# Prior only
ACAT_prior_Mod <-
        brm(Mod_fmla,
            data = Data,
            family = "acat",
            prior = priors_acat,
            inits = 0,
            iter = 1000,
            warmup = 500,
            chains = 2,
            cores = ncores,
            sample_prior = "only",
            control = list(adapt_delta = 0.8,
                           max_treedepth = 12)
        )

Iā€™m attempting some visual checks of the ā€œprior-onlyā€ model with the following:

ppc_dens_overlay(y = Emo_Data$Valence, 
                 yrep = posterior_predict(Valence_BRMS_acat_prior_Mod, 
                                          nsamples = 25) %>% 
                         na.omit()) +
        scale_x_continuous(limits = c(0,9), breaks = seq(0,9,1)) +
        theme_bw() + 
        theme(panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              axis.line = element_line(colour = "black")
              )

But, the shape of the prior distributions is way off:

prior_check
prior_check724Ɨ620 80.2 KB

Is there a way to specify the priors for category-specific intercepts to better capture the shape of the responses?

2

Just to clarify, in the plot, is the x-axis actually continuous? Or is this 10 different groups?

Would something like a grouped plot work better (look for the ā€˜ppc_stat_groupedā€™ example here)?

1 Like
3

I think the x-axis represents the 9 ordinal responses and the 8 cut-points/thresholds between them.

1 Like
4

I guess step 1 is figuring out the units on these things exactly then.

What is y here: y = Emo_Data$Valence?

5

The outcome variable (y) is a self-reported rating of valence, entered as a categorical factor (integers 1:9) in response to auditory and visual tokens.

My understanding is that model treats the outcome as a continuous latent variable with discrete cut-points at the thresholds. The structure of the family = "adjacent category" model seems similar to the family = "cumulative" model with thresholds set to ā€œequidistantā€.

If Iā€™m looking at it correctly, I think the plot is showing a valence rating of 5 is most common/probable and that the extremes are least likely.

In reading on possible options, it seems a dirichlet prior might be appropriate, but I get an error because it has to be specified as a vector and I canā€™t see how to set a prior as a vector with set_prior

6

Hi,

use pp_check(Valence_BRMS_acat_prior_Mod, type="hist")

7

Thank you. The histogram plots look like this for 10 posterior samples:

PP_Hist
PP_Hist724Ɨ620 19.4 KB

Still not a great match between Y and Yrep

8

Eh, did you use the call I gave you above, i.e., pp_check()? I was thinking that you should sample from your priors only and then plot a number of draws in one diagram, so

pp_check(Valence_BRMS_acat_prior_Mod, type = "hist", nsamples = 100)
9

Yes, thatā€™s the call I used. But it returns a separate facet for each sample. It defaults to 10 samples, if I pass nsamples = 100, I get 100 facets

10

And you ran brm() with sample_priors="only"?

11

Yes, my brm call:

Valence_BRMS_acat_prior_Mod <-
        brm(Emo_Mod_fmla,
            data = Emo_Data,
            family = "acat",
            prior = priors_acat,
            inits = 0,
            iter = 1000,
            warmup = 500,
            chains = 2,
            cores = ncores,
            sample_prior = "only",
            control = list(adapt_delta = 0.8,
                           max_treedepth = 12)
        )
12

Are you perhaps running bayesplot's pp_check? Try with brms::pp_check(). If that doesnā€™t work then I have no clue since it works here for me. Perhaps @bbbales2 knows if thereā€™s been some changes that I donā€™t know about :)

13

Iā€™m able to plot as a single figure with bayesplotā€™s type = "bars",

bayesplot_bars
bayesplot_bars724Ɨ620 10.8 KB

But, the priors still look off.

Is there a way to specify the threshold priors so they are better aligned with the observed data?

14

Ah, crap, Iā€™m really sorry @JLC this is all my fault. It should be bars of courseā€¦ :)

So, increase nsamples now. By changing the priors you would like to see a more uniform behavior from the priors. To me it looks strange that some categories are predicted to have very low counts (6,7,8) and some (0 and 9) to have very high counts (but larger intervals too).

Check also if itā€™s ok to start with a ā€˜0ā€™ as your data does. I know that sometimes a count must start from ā€˜1ā€™.

2 Likes
15

No problem! Any and all guidance is incredibly helpful!

Iā€™m now mostly confused as to how the priors can/should be specified.

I have tried setting global priors for the betas, intercepts, and sdā€™s, but still see high predicted probabilities at the extremes (1 and 9) and low predictions for everything in between.

priors_global <- c(
        set_prior("student_t(3, 0, 2.5)",
                  class = "b"),
        
        set_prior("normal(0,2)",
                  class = "Intercept"),
        
        set_prior("student_t(3, 0, 2.5)",
                  class = "sd")
)

pp_check_global
pp_check_global724Ɨ620 10.4 KB

How can I shape the priors so that I suggest lower predicted probability for the tails and more for the middle?

16

Ah, yes, thereā€™s really not one answer to this. I would recommend you to set some priors (like youā€™ve done now) and then start seeing what happens when you, e.g., start increase the deviation on the intercept, then on the \beta params, then on the SD.

17

BTW, make sure to have a fairly large nsamples set since youā€™d like to see the uncertainty properly.

1 Like
18

Thanks again!

Iā€™m working with 100 samples, for now.

I think the main challenge is that I donā€™t have a solid grasp of the mathematical relationship between the intercepts, so itā€™s difficult to construct anything that gives me any sort of uniform predictions over the response choices. Everything I try ends up with a horseshoe shape.

19

Continue analyzing the priors and see what happens. A mathematical relationship is usually too complex for me to parse, hence I prefer to plot.

20

Is there any way to remove NAā€™s from the predictions?