Example 8.17: Logistic regression via MCMC
[This article was first published on SAS and R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
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.
To leave a comment for the author, please follow the link and comment on their blog: SAS and R.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.