Monotonic effects in non-gaussian models

Interfaces brms
1

Please also provide the following information in addition to your question:

  • Operating System: Win 10
  • brms Version: 2.6.0

Hi, I have a question how to interprete / deal with monotonic effects in models, e.g. with count-outcome.

There’s a paper describing how to model these effects with brms, using a linear model. For the linear model, the parameter b of the monotonic effect indicates direction and size (range between lowest and highest category of the ordinal predictor), while the simplex parameters indicate the “normalized” distance between each category.

How would this be interpreted, e.g., in a count-model? Especially if I would like to report incident rate ratios, I would usually exponentiate the “point estimate” (e.g. posterior median). Can I also exponentiate the parameter b of the monotonic effect, and does it then tell me the range of the “ratio”?

Example:

library(brms)
income_options <- c("below_20", "20_to_40", "40_to_100", "greater_100")
income <- factor(sample(income_options, 100, TRUE), levels = income_options, ordered = TRUE)
mean_ls <- c(30, 60, 70, 75)
ls <- mean_ls[income] + rnorm(100, sd = 7)
dat <- data.frame(income, ls, count = rpois(100, 2))

m1 <- brm(ls ~ mo(income), data = dat)
m2 <- brm(count ~ mo(income), data = dat, family = poisson)

summary(m1)
#> Population-Level Effects: 
#>           Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
#> Intercept    29.97      1.79    26.55    33.40       2686 1.00
#> moincome     46.18      2.47    41.27    50.87       2591 1.00
#> 
#> Simplex Parameters: 
#>              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
#> moincome1[1]     0.63      0.04     0.55     0.71       3690 1.00
#> moincome1[2]     0.20      0.05     0.11     0.29       3897 1.00
#> moincome1[3]     0.17      0.04     0.08     0.25       2750 1.00

summary(m2)
#> Population-Level Effects: 
#>           Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
#> Intercept     0.69      0.14     0.40     0.94       1516 1.00
#> moincome      0.12      0.20    -0.26     0.53       1488 1.00
#> 
#> Simplex Parameters: 
#>              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
#> moincome1[1]     0.35      0.23     0.01     0.84       2806 1.00
#> moincome1[2]     0.31      0.22     0.01     0.81       2707 1.00
#> moincome1[3]     0.34      0.23     0.01     0.85       2870 1.00

For the second model, the posterior mean for moincome is 0.12, which is on the log-scale. exp(0.12) ~ 1.13, so the IRR of moincome would be 1.13. Can I conclude that the ratio for the monotonic effect has a range up to 1.13 times higher for the highest category?

This is how the marginal effects plot looks like, which automatically transforms the predictions into the scale of the outcome:

marginal_effects(m2, "income")

Rplot

So, how would I interprete the (exponentiated) coefficient of moincome for the poisson-model? Or does it only make sense to interprete these coefficients by looking at marginal effects plots?

2

Exponentiating makes everything multiplicative that was additive before. This rule of course applies to monotonic effects as well. This has to be kept in mind.

If b is the range of the monotonic effect on the log-scale, then

\exp(\eta + b) - \exp(\eta)

is the range on the expoentiated scale where \eta is the rest of the predictor term (i.e. other predictors) you condition on.

When transforming summary statistics, please keep in mind that this does not work for all statistics. For instance, the mean is not equivariant under non-linear transformations.
It is best to transform the posterior samples and then summarize only after having applied all transformations.

3

Thanks, Paul! This sounds like it is easier to interprete the monotonic effect via marginal effects plots, instead showing a classical table with regression coefficients.

When transforming summary statistics, please keep in mind that this does not work for all statistics.

Is this related to Jensen’s inequality?!?

So you would suggest using

m2 %>% 
  posterior_samples() %>% 
  exp() %>% 
  purrr::map_dbl(mean)

and not

m2 %>% 
  posterior_samples() %>% 
  purrr::map_dbl(mean) %>% 
  exp()

to get transformed “estimates”?

4

I like plotting things more anyway ;-)

The complications because of the transformations are not worse for monotonic effects than for everything else. So if you consider plots easier for monotonic effects that you may consider them easier for all kinds of effects.

Yeah I am sort of suggesting that. Whatever quantity you are interested in should be compute per posterior sample and summarized afterwards.

5

I like plotting things more anyway ;-)

Me too, especially for everything that includes interactions or in general models of non-gaussian families. ;-)