| This post was kindly contributed by SAS and R - go there to comment and to read the full post. |
In examples 8.15 and 8.16 we considered Firth logistic regression and exact logistic regression as ways around the problem of separation, often encountered in logistic regression. (Re-cap: Separation happens when all the observations in a category share a result, or when a continuous covariate predicts the outcome too well. It results in a likelihood maximized when a parameter is extremely large, and causes trouble with ordinary maximum likelihood approached.) Another option is to use Bayesian methods. Here we focus on Markov chain Monte Carlo (MCMC) approaches to Bayesian analysis.
SAS
SAS access to MCMC for logistic regression is provided through the bayes statement in proc genmod. There are several default priors available. The normal prior is the most flexible (in the software), allowing different prior means and variances for the regression parameters. The prior is specified through a separate data set. We begin by setting up the data in the events/trials syntax. Then we define a fairly vague prior for the intercept and the effect of the covariate: uncorrelated, and each with a mean of zero and a variance of 1000 (or a precision of 0.001). Finally, we call proc genmod to implement the analysis.
data testmcmc;
x=0; count=0; n=100; output;
x=1; count=5; n=100; output;
run;
data prior;
input _type_ $ Intercept x;
datalines;
Var 1000 1000
Mean 0 0
;
run;
title "Bayes with normal prior";
proc genmod descending data=testmcmc;
model count/n = x / dist=binomial link=logit;
bayes seed=10231995 nbi=1000 nmc=21000
coeffprior=normal(input=prior) diagnostics=all
statistics=summary;
run;
In the forgoing, nbi is the length of the burn-in and nmc is the total number of Monte Carlo iterations. The remaining options define the prior and request certain output. The diagnostics=all option generates many results, including posterior autocorrelations, Gelman-Rubin, Geweke, Raftery-Lewis, and Heidelberger-Welch diagnostics. The summary statistics are presented below; the diagnostics are not especially encouraging.
Posterior Summaries
Standard
Parameter N Mean Deviation
Intercept 21000 -20.3301 10.3277
x 21000 17.2857 10.3368
Posterior Summaries
Percentiles
Parameter 25% 50% 75%
Intercept -27.6173 -18.5558 -11.9025
x 8.8534 15.5267 24.6024
It seems that perhaps this prior is too vague. Perhaps we can make it a little more precise. A log odds ratio of 10 implies an odds ratio > 22,000, so perhaps we can accept a prior variance of 25, with about 95% of the prior weight between -10 and 10.
data prior;
input _type_ $ Intercept x;
datalines;
Var 25 25
Mean 0 0
;
run;
ods graphics on;
title "Bayes with normal prior";
proc genmod descending data=testmcmc;
model count/n = x / dist=binomial link=logit;
bayes seed=10231995 nbi=1000 nmc=21000
coeffprior=normal(input=prior) diagnostics=all
statistics=(summary interval) plot=all;
run;
Posterior Summaries
Standard
Parameter N Mean Deviation
Intercept 21000 -6.5924 1.7958
x 21000 3.5169 1.8431
Posterior Summaries
Percentiles
Parameter 25% 50% 75%
Intercept -7.6347 -6.3150 -5.2802
x 2.1790 3.2684 4.5929
Posterior Intervals
Parameter Alpha Equal-Tail Interval HPD Interval
Intercept 0.050 -10.8101 -3.8560 -10.2652 -3.5788
x 0.050 0.5981 7.7935 0.3997 7.4201
These are more plausible values, and the diagnostics are more favorable.
In the above, we added the keyword interval to generate confidence regions, and used the ods graphics on statement to enable ODS graphics and the plot=all option to generate the graphical output shown above.
R
There are several packages in R that include MCMC approaches. Here we use the MCMCpack package, which include the MCMClogit() function. It appears not to accept the weights option mentioned previously, so we generate data at the observation level to begin. Then we run the MCMC.
events.0=0 # for X = 0
events.1=5 # for X = 1
x = c(rep(0,100), rep(1,100))
y = c(rep(0,100-events.0), rep(1,events.0),
rep(0, 100-events.1), rep(1, events.1))
library(MCMCpack)
logmcmc = MCMClogit(y~as.factor(x), burnin=1000, mcmc=21000, b0=0, B0=.04)
The MCMClogit() accepts a formula object and allows the burn-in and number of Monte Carlo iterations to be specified. The prior mean b0 can be specified as a vector if different for each parameter, or as a scalar, as show. Similarly, the prior precision B0 can be a matrix or a scalar; if scalar, the parameters are uncorrelated in the prior.
> summary(logmcmc)
Iterations = 1001:22000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 21000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
(Intercept) -6.570 1.816 0.01253 0.03139
as.factor(x)1 3.513 1.859 0.01283 0.03363
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
(Intercept) -10.6591 -7.634 -6.299 -5.229 -3.831
as.factor(x)1 0.6399 2.147 3.292 4.599 7.698
plot(logmcmc)
The result of the plot() is shown below. These results and those from SAS are reassuringly similar. Many diagnostics are available through the coda package. The codamenu() function allows simple menu-based access to its tools. 
| This post was kindly contributed by SAS and R - go there to comment and to read the full post. |