Specifying Priors for Ordinal (adjacent category) Model

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

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

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)
        )

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 :)

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

But, the priors still look off.

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

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

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")
)

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

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.

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

1 Like

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.

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

Is there any way to remove NA’s from the predictions?

Does this perhaps happen for some type of priors? If yes that’s an indication of prior combinations that are not good.

It did, indeed.

I did find a workaround by calling the bayseplot function directly and omitting NAs

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

Ah, but be careful. The NAs could indicate that your priors are too broad!

See this post where Aki explains is nicely:

Thanks!

I did read that post when searching for possible approaches. I’m just filtering the NAs as a hack to try to see what, if any, changes result from modifying the priors.

It seems no matter what I do, there’s a horseshoe:

priors_acat <- c(
        # set_prior("student_t(3, 0, 2.5)",
        #           class = "b"),
        
        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(-2,0.05)",  
                  class = "Intercept",
                  coef = "1"),
        
        set_prior("normal(0,0.1)",  
                  class = "Intercept",
                  coef = "2"),
        
        set_prior("normal(0,0.1)",  
                  class = "Intercept",
                  coef = "3"),
        
        set_prior("normal(0,0.1)",  
                  class = "Intercept",
                  coef = "4"),
        
        set_prior("normal(0,0.1)",  
                  class = "Intercept",
                  coef = "5"),
        
        set_prior("normal(0,0.1)",  
                  class = "Intercept",
                  coef = "6"),
        
        set_prior("normal(0,0.1)",  
                  class = "Intercept",
                  coef = "7"),
        
        set_prior("normal(2, 0.05)",  
                  class = "Intercept",
                  coef = "8"),
        
        set_prior("cauchy(0,1)",  
                  class = "sd",
                  group = "participant",
                  coef = "Intercept"),
        
        set_prior("cauchy(0,1)",  
                  class = "sd",
                  group = "group:condition",
                  coef = "Intercept")
        
)

Gives this:

And a large change (increasing the location from 0 to 5 for intercepts 4,5,and 6), still gives the same shape on predictions:

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(2,2)",  
                  class = "Intercept",
                  coef = "1"),
        
        set_prior("normal(2,1)",  
                  class = "Intercept",
                  coef = "2"),
        
        set_prior("normal(4,1)",  
                  class = "Intercept",
                  coef = "3"),
        
        set_prior("normal(5,1)",  
                  class = "Intercept",
                  coef = "4"),
        
        set_prior("normal(5,1)",  
                  class = "Intercept",
                  coef = "5"),
        
        set_prior("normal(1,1)",  
                  class = "Intercept",
                  coef = "6"),
        
        set_prior("normal(1,1)",  
                  class = "Intercept",
                  coef = "7"),
        
        set_prior("normal(-2,2)",  
                  class = "Intercept",
                  coef = "8"),
        
        set_prior("cauchy(0,1)",  
                  class = "sd",
                  group = "participant",
                  coef = "Intercept"),
        
        set_prior("cauchy(0,1)",  
                  class = "sd",
                  group = "group:condition",
                  coef = "Intercept")        
)

Hi

looking at your priors I have to say they look a bit strange. Some are negative some positive, some has a mean that differs a lot from others, etc. Additionally you seem to set a prior for each chef when very often one sets (if one doesn’t have prior information!) one type of prior for all, e.g., \beta, \alpha, etc.

Could you put up your model here, with data if possible so one can try this? Also, acat could be an option, but have you looked into cumulative?

1 Like

Thanks again!

I found that when setting a global prior for the intercepts, I ended up with the horseshoe pattern in predictions, so I started trying to specify each intercept directly.

The data are attached here:
Emotion_Data.csv (434.6 KB)

My model is specified as:

# List of needed packages
Pkgs <- c("tidyverse", "magrittr", "brms", 
          "strengejacke", "rstan", "bayesplot", "parallel", 
          "bayestestR", "tidybayes", "modelr")

# Load packages
lapply(Pkgs, require, c = T)

## Set computing options
ncores = detectCores()
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

# Read in data
Emo_Data <- read.csv("Emotion_Data.csv") %>%
        mutate(group = factor(group, 
                              levels = c("NH", "HI")),
               condition = factor(condition, 
                                  levels = c("60", "80", "CP", "NFC")),
               new_category = factor(new_category, 
                                     levels = c("unpleasant_high", "unpleasant_low", 
                                                "neutral", "pleasant_low", "pleasant_high"))) %>%
        filter(!is.na(Valence))

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

# Set priors
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(1, 0.01)",  
                  class = "Intercept",
                  coef = "1"),
        
        set_prior("normal(2,0.01)",  
                  class = "Intercept",
                  coef = "2"),
        
        set_prior("normal(3,0.01)",  
                  class = "Intercept",
                  coef = "3"),
        
        set_prior("normal(4,0.01)",  
                  class = "Intercept",
                  coef = "4"),
        
        set_prior("normal(5,0.01)",  
                  class = "Intercept",
                  coef = "5"),
        
        set_prior("normal(6,0.01)",  
                  class = "Intercept",
                  coef = "6"),
        
        set_prior("normal(7,0.01)",  
                  class = "Intercept",
                  coef = "7"),
        
        set_prior("normal(8,0.01)",  
                  class = "Intercept",
                  coef = "8"),
        
        set_prior("cauchy(0,1)",  
                  class = "sd",
                  group = "participant",
                  coef = "Intercept"),
        
        set_prior("cauchy(0,1)",  
                  class = "sd",
                  group = "group:condition",
                  coef = "Intercept"),
        
        set_prior("cauchy(0,1)",  
                  class = "sd",
                  group = "new_category:token",
                  coef = "Intercept")        
)


Valence_BRMS_acat_Mod <-
        brm(Emo_Mod_fmla,
            data = Emo_Data,
            family = "acat",
            prior = priors_acat,
            inits = 0,
            iter = 2000,
            warmup = 1000,
            chains = 2,
            cores = ncores,
            control = list(adapt_delta = 0.9,
                           max_treedepth = 12)
        )

With the priors for intercepts set fairly tightly to the response category boundaries, I get the following with a prior-only predictive check:

edit to note that the model takes a very long time to run!

As @torkar pointed out, you might be better off specifying the cumulative family as opposed to acat. The way you’ve set up your priors on the thresholds would make sense if they were partitioning a latent variable (i.e. what cumulative does). I’ve no idea what the priors mean in relation to the acat model, where afaik you’re modeling the relative probability between category k and k + 1.

Additionally, your priors for standard deviations are cauchy, and the thick tails might be causing problems. Thinking about it, that could be the reason why you get a lot of extreme observations via prior predictive checks - because your model might be putting a lot of density into some really weird areas of the parameter space. Try setting the priors to normal or student-t with 3+ degrees of freedom and see what happens.