Strange behavior with scale() and poly() in brms model

Hello everyone, I’m trying to figure out how to include a polynomial in one of my models (related to a previous post), and it seems to have some strange behavior when it interacts with scale().

I have a predictor that goes from 1:10, and I would like to model both linear and quadratic effects of this predictor. When using poly(variable, 2) in the model, it seems to work all right-- the model converges, conditional_effects gives out reasonable predictions, etc. However, the estimated coefficients are quite large, since the poly() orthogonalization makes the predictor values much smaller than their original values.

I have tried to scale the predictors using scale(poly(variable,2)) in the model, and this also converges and now gives coefficients on par with my other predictors (e.g. between 0 and 1, ish). However, when I try to use conditional_effects() to get predictions, I get this error:

Error in poly(position, 2) : 
  'degree' must be less than number of unique points

From looking on other forums, it seems like this usually happens when the degree of the polynomial is unreasonably high, but here it is only 2. Also, since brms will happily run the model and converge, I am guessing the problem is with conditional_effects(). Any thoughts?

Note: my model is a logistic one; based on my earlier model it looked like conditional_effects() has some functionality for converting back to the original scale of the variables, so perhaps this is interfering in some way?

  • Operating System: Mac OS Catalina
  • brms Version: 2.15.0

brms does some magic to make functions such a poly and scale work on new data (because these functions always need the old data as well to recover the scaling). I assume something failes when combining the two. Can you provide a minimal reprex on Issues · paul-buerkner/brms · GitHub?

Hello Paul,

Thanks for looking into this! I have posted the min reprex there; I hope it is sufficient (I am not used to writing issue posts on Github). Let me know if you need anything else from me!

Hi Paul,

Is there a fix for this bug on the horizon, or another way to get the transformed statistics (estimate, cred. interval) in probability space instead of logit space? I’ll need to report some stats for a paper soon.

You can use hypothesis() to nonlinearly transform any combination of parameters of interest and this will report the estimate and CIs etc. E.g., for a logistic model with the formula Y ~ 1 + X + Z, if you want the probability for a condition identified by the Intercept + XLevel2 + ZLevel3:

hypothesis(your_mdl, "inv_logit(Intercept + XLevel2 + ZLevel3) > 0")

Thanks @andymilne, looks like fixef() will also do what I want. Sorry to bother!

If you use fixef(), make sure to transform the draws before computing the summary statistic. An advantage of using hypothesis() is that this is all taken care of “under-the-hood”.