How to improve priors for a modified Beta regression (custom ordered beta regression)?

Modeling
, ,
1

I want to improve my modeling by using weakly informative priors, but I struggle with logit conversions (and probably have many other knowledge gaps that are less obvious to me). I have fitted the model, but I still have the sense that I’m doing something wrong with my priors and experimenting with priors has given me substantially different bulk-ESS.

I am using ‘brms’, trying to fit a custom family (ordered beta regression, as specified here).

I don’t think I can share my data, but I will describe my variables: y (varies from 0 to 1), category (factor two levels), judge (factor 455 levels), target (factor 48 levels). The targets are nested into factor (but all have a unique identifier). My formula looks like this:

bf(y ~ 0 + Intercept + category + (0 + category | judge) + (1 | gr(target, by= category) ) )

The 0 + Intercept is a quirk of how the family is specified. I believe the specifics of the custom family & stan code are relevant:

ord_beta_reg <- custom_family(
  "ord_beta_reg", 
  dpars = c("mu","phi", "cutzero", "cutone"), 
  links = c("identity", "log", NA, NA), 
  lb = c(NA, 0, NA, NA), 
  type = "real")


real ord_beta_reg_lpdf(real y, real mu, real phi, real cutzero, real cutone) {   
  vector[2] thresh;   
  thresh[1] = cutzero;   
  thresh[2] = cutzero + exp(cutone); 
  
  if(y==0) {      
    return log1m_inv_logit(mu - thresh[1]);  
  } else if(y==1) {     
    return log_inv_logit(mu - thresh[2]);  
  } else {      
    return log_diff_exp(log_inv_logit(mu   - thresh[1]), log_inv_logit(mu - thresh[2])) +                
      beta_lpdf(y|exp(log_inv_logit(mu) + log(phi)),exp(log1m_inv_logit(mu) + log(phi)));    
  }  
      }

My latest thinking regarding weakly informative priors for the scale of my data is as follows:

# max effect = max(error) - min(error) / max(b) - min(b) -> in qlogis most reasonable values ~10/10 = 1
# divide max effect by 2.58
# because 99% distribution is within about 2.58 standard units from the mean
# 1/2.58 = .38
# inspired by http://svmiller.com/blog/2021/02/thinking-about-your-priors-bayesian-analysis/
For b -> "normal(0, 0.38)"

# Intercept qlogis(mean(y)), 2.5*sd(y)
For b Intercept -> "normal(-0.89, 0.58)"

# skeptical about random effect correlation being close to 1, -1
For random effect cor -> "lkj(4)"

# custom induced Dirichlet priors for cutpoints
cutone -> "induced_dirichlet([1,1,1]', 0, 2,cutone,cutzero)"
cutzero -> "induced_dirichlet([1,1,1]', 0, 1,cutone,cutzero)"

# precision parameter for Beta, taken from Kubinec (2022)
phi -> "exponential(.1)"

# standard deviation for all group effects, gamma(2, 1/sd(y)), taken from Chung Et Al. (2013)
sd -> "gamma(2, 4.34)"

I’d appreciate comments on whether the rationale for constraining my priors is sensible, as well as links to pertinent resources. Thank you!

Kubinec (2022) → Ordered Beta Regression: A Parsimonious, Well-Fitting Model for Continuous Data with Lower and Upper Bounds | Political Analysis | Cambridge Core
Chung Et Al. (2013) → A Nondegenerate Penalized Likelihood Estimator for Variance Parameters in Multilevel Models | SpringerLink

2

Generally speaking, you only want to constrain priors if 1) you have issues of identifiability in the model or 2) you have a theoretical reason to use strong priors. Otherwise, weakly informative priors are superior for the stability of inference.

I would encourage you to use the R package ordbetareg so that you can use standard forms of syntax and additional features of the distribution.