I have fitted a categorical logistic regression in rstanarm. I’m now looking to obtain posterior estimates for odds ratios and probabilities.
Assuming the Stan model is in a variable called fit, getting point estimates is easy:
odds <- exp(fit$coefficients)
probs <- odds / (1 + odds)
The same can’t be said for credible intervals. Typically, one could exponentiate the interval for the log odds; however, credible intervals aren’t usually invariant under monotonic transformation. If they were, one could do:
odds_cri <- exp(posterior_interval(fit, prob = 0.95))
prob_cri <- odds_cri / (1 + odds_cri)
Posterior SDs are another story altogether.
The model could be fitted in Stan itself via rstan, but more work from the modeller will be needed to set things up for valid inference.
Is there an easy way to do this in rstanarm? I’ve not been able to find anything in the documentation (vignettes, etc.).
I’m not sure what you mean here. If you have a credible interval (a, b) for a random variable X, and you have a strictly monotonic upward function f, then the then you know the credible interval of the same coverage for Y = f(X) is indeed (f(a), f(b)). The densities change, but the boundaries adapt to handle it.
This is perhaps easiest to see from the MCMC estimators. Let’s say I have M draws
x_1, \ldots, x_M \sim p_X(x).
We estimate the central posterior interval at \alpha \in (0, 1) coverage by sorting the elements of x_{1:M} into z_1, \ldots, z_M in ascending order and using \left(z_{\lfloor M \cdot \alpha / 2 \rfloor}, z_{M - \lfloor M \cdot \alpha / 2 \rfloor}\right) as our estimate of the bounds of the central \alpha posterior interval of X.
Now take a strictly monotone increasing f:\mathbb{R} \rightarrow \mathbb{R} (i.e., a < b implies f(a) < f(b)). Then my transformed draws are
f(x_1), \ldots, f(x_M) \sim p_Y(y)
and it turns out the sort of f(x_1), \ldots, f(x_M) is just f(z_1), \ldots, f(z_M) because f is monotonic. Hence my estimated bounds for the central \alpha interval of Y = f(X) is \left(f(z_{\lfloor M \cdot \alpha / 2 \rfloor}), f(z_{M - \lfloor M \cdot \alpha / 2 \rfloor})\right).
I’m sure it’s also possible to do with calculus by applying change of variables and then integrating between the new bounds. I’m just better at sampling than calculus :-).
Ah OK.
This is reassuring and makes things a lot easier. I had read a Cross Valiated post here: bayesian - Are HPD intervals invariant to reparameterization? - Cross Validated which piqued my attention.
I can now get all this information from the objects outputted via rstanarm::stan_glm(). However, is there a simple way to get the posterior standard deviations for the odds and probabilities?
I could write a Stan model from scratch and declare odds and probabilities within a generated quantities block, but that would likely be a lot of extra work. Hopefully rstanarm can output what I’m needing.