# How can I use credibility intervals in Bayesian logistic regression?

I'm using Gibbs sampling to learn the distributions of coefficients for a multinomial logistic regression model. At the end, I end up using the mean values of distributions of coefficients, and the resulting logistic regression is used as a classifier.

I'm trying to find out advantages of having probability distributions for coefficients and the response variable, but I can't really see the way to leverage credibility intervals. What can I do with these distributions that I end up with, other than using their mean values?

ok, I think I've failed to express my question clearly. Here is the update

Let's assume that I have y = b1*x1 + b2*x2 + b3*x3 in my hands. for all b1,b2,b3 I have normal distributions, with mean and quantiles. These distributions come from the mcmc results. I can use the means of b1,b2,b3 and that'd give me a classifier. I can also use %2.5 quantile values and 97.5 quantile values for b1,b2,b3 and so on.

What would be the probabilistic interpretation of these other equations? Can I produce a smoother classifier this way, rather than using only means? The use of credibility intervals is quite clear to me when they are used for a single variable, but in this case, I'm talking about N variables (5 actually) each with their own credibility intervals. I'm having trouble getting the semantics of this setup. I have not seen any papers etc that discusses this, and any pointers would be appreciated.

I wouldn't use the means at all for the classifier. You don't need to apply "corrections" or to "smooth out" a Bayesian solution, it is the optimal one for the prior information and data that you have actually used. But the means can be useful for giving you a feel for which combinations of regressor variables are likely to lead to classifying towards a particular category.

However this can be a horrendously complicated beast for multinomial regression, as you have a matrix of betas to interpret (one column for each category, except for the reference, which can be thought of as having all betas "estimated" as zero with zero standard error).

Given that this seems to be an attempt at a intuitive way to understand what your classifier is doing, let me propose another. I will delete this section if this is not what you were intending.

You have your MCMC samples of the beta matrix: call this $\beta_{ij}^{(b)}$ where $i=1,\dots,R$ denotes the multinomial category, $j=1,\dots,p$ denotes the regressor variable (the $X$), and $b=1,\dots,B$ denotes the $bth$ MCMC sampled value. If the categories have different $X$ variables, then simply set those excluded variables' betas to zero in the matrix: $\beta_{ik}^{(b)}=0$ for all $b$ if variable $k$ was not part of the model fit to the ith category, and $\beta_{Rj}^{(b)}=0$ for all $j$ and $b$.

The first thing you need is a set of covariates to use $X_{mj}\;\;\;\;m=1,\dots,M$, where $m$ is the "observation number" and $M$ is the number of predictions you are going to make. The data used to fit your model should do for this purpose, so $M=\text{sample size}$. You now calculate the linear predictor for each category for each prediction for each MCMC sample: $$y_{im}^{(b)}=\sum_{j=1}^{p}X_{mj}\beta_{ij}^{(b)}$$

(this may be quicker to code up as a matrix/array operation). Note that $y_{Rm}^{(b)}=0$ for all $b$ and $m$. Then convert this into a probability for the $mth$ observation belonging to the $ith$ category/class, call this quantity $Classify(m,i)$.

$$Classify(m,i)=\frac{1}{B}\sum_{b=1}^{B}\frac{\exp\left(y_{im}^{(b)}\right)}{\sum_{l=1}^{R}\exp\left(y_{lm}^{(b)}\right)}$$

Now you plot the value of $Classify(i,m)$ against $X_{mj}$, so you will have a total of $R\times p$ plots. Looking at these should give you a feel for what the classifier is doing in relation to the regressor variables.

Note that when it comes to actually classifying a new variable, you only need $Classify(i,m)$ in order to do this - all other quantities from the MCMC are irrelevant for the purpose of classification. What you do need though is a loss matrix which describes the loss incurred from classifying into category $i_{est}$ when the true category is actually $i_{true}$, this will be a $R\times R$ matrix, usually zero on the diagonal and positive everywhere else. This can be very important if correctly identifying "rare" classes is crucial compared to correctly identifying "common" classes.

• Let me try to rephrase it like this: I have a linear predictor for each category, where the betas are results of point estimations (when I use the mean values from their distributions or use any other maximization technique) How do I use probability distributions for betas rather than their point estimates in linear predictors? Or is there a way to use them at all? – mahonya May 30 '11 at 9:36
• @user3280 - the procedure I outlined does exactly that - the summation over the Gibbs samples (the b's) is an approximation to integrating over that parameter space. This is using the probability distribution of the betas to get predictions. But remember that the ultimate goal is prediction here, so you shouldn't really care about the betas per se, only about the predictions. the betas are a secondary concern. – probabilityislogic May 30 '11 at 19:30

Don't use the mean of the sampled coefficients for making predictions, instead compute the predictions for logistic regression models with all of the sampled coefficient vectors and take the mean of those predictions (or better still treat the predictions for all sampled coefficient vectors as the posterior distribution of the probability of class membership - the spread of that distribution is a useful indicator of how confident the classifier is about the probabilistic classification).

The distribution of the sampled coeffcient vector gives an impression of how well the training data constrain the value of that parameter, so if the distribution is broad, we can't be confident of the "true" value of that coefficient (as explained by Manoel Galdino). However, the key advantage of having a distribution of plausible coefficient vectors is that it provides you with a rational way to get a distribution of plausible values for the probability of class membership, which is what we really want.

Often using a Bayesian approach, we are not really interested in the coefficients of the model, but in the function implemented by the model.

• I've updated the question, also could you please explain your last sentence a little bit? Thanks! – mahonya May 29 '11 at 11:34

I don't know if I am understanding correctly your question. But I guess you may use the posterior density to assess the uncertainty around point estimates like the mean. You may plot a histogram, calculate standard deviations. This is easy to do, if you have the MCMC output. Just take the values sampled (after a burnin period) and compute the means, standard deviations and plot histograms or densities.

Another advantage of a posterior distribution is that you can assess your uncertainty on tails of the distribution and also if the distribution is simetric around the mode/mean. Confidence intervals in general assume that the distribution is simetric and that outliers are rare...

What you're looking for, and what the other respondants have proposed, is called the posterior predictive distribution. It takes into account the inherent uncertainty of the parameter estimates.

You can either use the samples from the MCMC run, or you can approximate it from the mean and covariance of the posterior distribution of the parameters by use of the probit function. See pages 218-220 of Chris Bishop's book "Pattern Recognition and Machine Learning" for an overview of how this can be done.