Journal of Statistical Software
June 2011,Volume 42,Issue 9.https://www.wendangku.net/doc/618588846.html,/
MCMCpack :Markov Chain Monte Carlo in R
Andrew D.Martin
Washington University
in St.Louis Kevin M.Quinn University of California,Berkeley Jong Hee Park University of Chicago
Abstract
We introduce MCMCpack ,an R package that contains functions to perform Bayesian
inference using posterior simulation for a number of statistical models.In addition to
code that can be used to ?t commonly used models,MCMCpack also contains some
useful utility functions,including some additional density functions and pseudo-random
number generators for statistical distributions,a general purpose Metropolis sampling
algorithm,and tools for visualization.
Keywords :Bayesian inference,Markov chain Monte Carlo,R .
1.Introduction
The Bayesian paradigm for statistical inference is appealing to researchers on both theoret-ical (Je?reys 1998;Bernardo and Smith 1994)and practical (Gelman,Carlin,Stern,and Rubin 2003)grounds.The interest in Bayesian methods in the social sciences is growing,and a number of researchers are using these approaches in substantive applications dealing with:deliberative bodies (Clinton,Jackman,and Rivers 2004;Martin and Quinn 2002),eco-nomic performance (Western 1998),income dynamics (Geweke and Keane 2000),legislative redistricting (Gelman and King 1990),mass voting (King,Rosen,and Tanner 1999),party competition (Quinn,Martin,and Whitford 1999;Quinn and Martin 2002),social networks (Ho?,Raftery,and Handcock 2002;Ho?and Ward 2004),and historical changes (Park 2010).Despite this interest,the social scienti?c community has yet to take full advantage of Bayesian techniques when approaching research problems.This is due primarily to two distinct prob-lems.The ?rst problem,which has been largely solved,was the inability to compute the high dimensional integrals necessary to characterize the posterior distribution for most mod-els.To a large extent,this has been solved by the advent of Markov chain Monte Carlo (MCMC)methods (Tanner and Wong 1987;Gelfand and Smith 1990;Besag,Green,Higdon,
2MCMCpack:Markov Chain Monte Carlo in R
and Mengersen1995)and the dramatic increases in computing power over the past twenty years.For a comprehensive treatment of MCMC methods,see Robert and Casella(2004).
MCMC methods are widely considered the most important development in statistical com-puting in recent history.MCMC has allowed statisticians to?t essentially any probability model—including those not even considered a few years ago.Unfortunately,statisticians have been about the only people who have been willing,and able,to write the computer code necessary to use MCMC methods to?t probability models.
Which brings us to the second problem;namely,the lack of?exible,yet easy-to-use software for social scientists unwilling or unable to devote substantial time and energy writing custom software to?t models via MCMC.Since reasonably e?cient MCMC algorithms exist to sample from the posterior distribution for most classes of models,developing software to meet the needs of social scientists is feasible.
MCMCpack(Martin,Quinn,and Park2011)is an R(R Development Core Team2011b) package that contains functions to perform Bayesian inference.It provides a computational environment that puts Bayesian tools(particularly MCMC methods)into the hands of social science researchers so that they(like statisticians)can?t innovative models of their choosing. Just as the advent of general purpose statistical software(like SPSS and SAS)on mainframe and then personal computers led to the widespread adoption of statistical approaches in the social sciences,providing easy-to-use general purpose software to perform Bayesian inference will bring Bayesian methods into mainstream social science.
MCMCpack currently contains the eighteen statistical models:linear regression models(linear regression with Gaussian errors,a singular value decomposition regression,and regression for a censored dependent variable),discrete choice models(logistic regression,multinomial logistic regression,ordinal probit regression,and probit regression),measurement models (a one-dimensional IRT model,a k-dimensional IRT model,a k-dimensional ordinal factor model,a k-dimensional linear factor model,a k-dimensional mixed factor model,and a k-dimensional robust IRT model),a model for count data(a Poisson regression model),models for ecological inference(a hierarchical ecological inference model and a dynamic ecological inference model),and time-series models for change-point problems(a binary change-point model,a probit change-point model,an ordinal probit change-point model,and a Poisson change-point model).Many of these models,especially the measurement models,are otherwise intractable unless one uses MCMC.
The package also contains the density functions and pseudo-random number generators for the Dirichlet,inverse Gamma,inverse Wishart,noncentral Hypergeometric,and Wishart dis-tributions.These functions are particularly useful for visualizing prior distributions.Finally, MCMCpack contains a number of utility functions for creating graphs,reading and writing data to external?les,creating mcmc objects,and manipulating variance-covariance matrices. The coda package is currently used for posterior visualization and summarization(Plummer, Best,Cowles,and Vines2006).MCMCpack is available from the Comprehensive R Archive Network at https://www.wendangku.net/doc/618588846.html,/package=MCMCpack.
The remainder of this paper is organized as follows.In Section2we discuss the package environment and features of MCMCpack.The following sections will show how model?tting functions in MCMCpack are implemented in a Gaussian linear model,and a Poisson change-point model.We conclude with a discussion of possible MCMCpack future developments.
Journal of Statistical Software3
2.The MCMCpack environment
We have chosen to make the R system for statistical computation and graphics the home environment for our software.R has a number of features that make it ideal for our purposes. It is open-source,free software that is distributed under the GNU General Public License.It is an extremely powerful programming environment that is designed for statistical data analysis and data visualization.R is very similar to the S language(Becker,Chambers,and Wilks 1988),and provides an easy-to-use interface to compiled C,C++or Fortran code,as well as excellent facilities for debugging,pro?ling,and documenting native R programs.R is already the general purpose tool of choice for most applied statisticians and is well documented and supported(Venables and Ripley2000,2002;R Development Core Team2011a,c).Evidence of the move toward R among social scientists can be found in the growing number of texts designed for social science graduate students that explicitly advocate R usage(Fox2002;Gill 2002),and the decision of the Inter-University Consortium for Political and Social Research to require students to use R as a means of integrating its most advanced courses in their summer training program(Anderson,Fox,Frankin,and Gill2003).Over the last?ve years R has become the lingua franca of applied statisticians working in the social sciences.
2.1.Design philosophy
In building the MCMCpack package we have attempted to adhere to a design philosophy that emphasizes:(a)widespread availability;(b)model-speci?c,computationally e?cient MCMC algorithms;(c)compiled C++code to maximize computational speed;(d)an easy-to-use, standardized model interface that is very similar to the standard R model?tting functions; and(e)compatibility with existing R packages(such as coda)for convergence assessment, posterior summarization,and data visualization.
From a purely practical perspective,the most important design goal has been the imple-mentation of MCMC algorithms that are model-speci?c.The major advantage of such an approach is that the sampling algorithms,being hand-crafted to particular classes of models, can be made dramatically more e?cient than black box approaches,such as those found in the WinBUGS software,while remaining robust to poorly conditioned or unusual data.In addition to using reasonably computationally e?cient sampling algorithms,the MCMCpack model?tting functions are also designed to be fast implementations of particular algorithms. To this end,nearly all of the actual MCMC sampling takes place in C++code that is called from within R.
The model?tting functions in MCMCpack have been written to be as similar as possible to the corresponding R functions for classical estimation of the models in question.This largely eliminates the need to learn a specialized model syntax for anyone who is even a novice user of R.For instance,to?t a linear regression model in R via ordinary least squares one could use the following syntax:
reg.out<-lm(y~x1+x2+x3,data=mydata)
Bayesian?tting of the same model with an improper uniform prior on the coe?cient vector and the default settings for the parameters governing the MCMC algorithm can be accomplished with MCMCpack using the following syntax:
Bayes.reg.out<-MCMCregress(y~x1+x2+x3,data=mydata)
4MCMCpack:Markov Chain Monte Carlo in R
It is our experience that such similarities greatly decrease the amount of time it takes to become a competent user of the MCMCpack package.
In addition,the MCMCpack model?tting functions are designed to be as similar to each other as possible.As a result,once a user becomes familiar with the basics of using one model ?tting function she can quickly move on to con?dently use other model?tting functions.For instance,to?t a probit model in which theβvector is assumed to have a multivariate normal prior with mean m and precision P one would use the MCMCpack syntax:
Bayes.probit.out<-MCMCprobit(y~x1+x2+x3,b0=m,B0=P,
data=mydata)
To?t a Poisson regression model in which theβvector is assumed to have a multivariate normal prior with mean m and precision P one would use the MCMCpack syntax:
Bayes.poisson.out<-MCMCpoisson(y~x1+x2+x3,b0=m,B0=P,
data=mydata)
Using R as the home environment for MCMCpack allows us to make use of a wide range of software for MCMC convergence assessment,posterior summarization,and data visualization (such as the coda package of Plummer et al.2006)that is already part of the R system. Output objects from all of the MCMCpack model?tting functions are formatted as coda mcmc objects,sometimes with additional attributes to allow for other types of analyses.This provides for a very easy-to-use interface between the two packages.For instance,to generate trace plots and density plots of the output in Bayes.poisson.out using the plot method, one would enter
plot(Bayes.poisson.out)
at the R command prompt.Similarly,to compute posterior means,standard deviations, quantiles,https://www.wendangku.net/doc/618588846.html,ing the summary method one would enter
summary(Bayes.poisson.out)
at the R command prompt.All other coda functions can used in conjunction with MCMCpack with similar ease.
https://www.wendangku.net/doc/618588846.html,parison with existing software
Perhaps the most widely used piece of software for applied Bayesian inference is the WinBUGS package,distributed by the MRC Biostatistics Unit at Cambridge(Spiegelhalter,Thomas, Best,and Lunn2003).WinBUGS is an important contribution,and is ideally suited for small models(both in the number of parameters and in the amount of data).WinBUGS is based on a set of general computational algorithms that can be used to estimate models speci?ed using a model de?nition syntax.WinBUGS oftentimes does not exploit model-speci?c information for estimation.All models are estimated using general purpose black-box algorithms.For the most part,sampling is done parameter-by-parameter.There are two other implementations of the BUGS language:JAGS(Plummer2003,2005)and OpenBUGS(Thomas,O’Hara,Ligges, and Sturtz2006).Furthermore,there are various R interfaces including R2WinBUGS(Sturtz,
Journal of Statistical Software5 Ligges,and Gelman2005),rjags(Plummer2009),or BRugs(Ligges,Thomas,Spiegelhalter,
Best,Lunn,Rice,and Sturtz2009).
While the BUGS language is useful for?tting many types of models,it has its limitations. The generic estimation engine used in WinBUGS is quite?exible,but because it does not exploit model-speci?c information,it is often ine?cient,and sometimes ine?ective.It is well known that WinBUGS has trouble with algorithms that require sampling from truncated distributions,such as the Albert and Chib(1993)algorithm for the binary probit model, and algorithms for?tting ordinal data models(Cowles1996;Johnson and Albert1999).The parameter-by-parameter approach is computationally ine?cient,both in terms of computa-tion time and clock time(Liu,Wong,and Kong1994).While we think MCMCpack has de?nite advantages over WinBUGS for many users,we emphasize that we view BUGS and MCMCpack as complementary tools for the applied researcher.In particular,the greater?ex-ibility of the BUGS language is perfect for users who need to build and?t custom probability models.
Other general purpose software packages with designs similar to that of MCMCpack exist. One is the BACC package(for Bayesian analysis,computation,and communication)discussed in Geweke(1999).The development of BACC was supported by the National Science Foun-dation,and the newest release was made available on2003-06-03.BACC is available as a computational engine for MATLAB,R,S-PLUS,GAUSS,and the UNIX or Windows console. BACC shares a design component with MCMCpack,as each algorithm is coded in compiled code for each individual model.Another general purpose software package is the recently released bayesm package by Rossi and McCulloch(2006).This package is geared to analyses in marketing and applied econometrics.Both BACC and bayesm use an interface where users provide appropriate scalars,vectors,and matrices for data and priors,and where posterior density samples are returned as matrices,requiring the user to perform additional computa-tion to summarize results,assess model?t,and the like.MCMCpack,on the other hand,has a more R-like interface.We also view these packages as complementary to MCMCpack,as they provide estimation engines for useful models.
In addition to these packages that o?er a comprehensive set of tools,there are a number of other available R packages that perform Bayesian inference for single classes of models.For a list of available packages,see the“Bayesian Inference”task view on the Comprehensive R Archive Network(Park,Martin,and Quinn2010).
2.3.An overview of MCMCpack features
In addition to the model?tting functions mentioned above,MCMCpack has numerous fea-tures aimed at both researchers and instructors.For the sake of space we do not comprehen-sively demonstrate this functionality here.Rather,we highlight three features that may be of interest.We encourage interested users to consult the documentation for further details. The MCMCmetrop1R function allows users to sample from a user-de?ned continuous density using a random walk Metropolis algorithm with a multivariate normal proposal distribution. This can be used to explore a posterior(or log-posterior)distribution,as well as any other target distribution of interest.See Gelman et al.(2003)and Robert and Casella(2004)for details of the random walk Metropolis algorithm.The sampler itself is coded in C++,so the iterations run quite https://www.wendangku.net/doc/618588846.html,ers only have to provide the target density as an R function. MCMCpack is sometimes used on large problems where parallel computation might a?ord
6MCMCpack:Markov Chain Monte Carlo in R
some advantages.While MCMCpack does not currently support parallelization within the Monte Carlo loop,for some problems it is useful to perform embarrassingly parallel simuations, e.g.,sampling from the posterior density of the same model with twenty-?ve simultaneous processes.Doing so requires the availability of a random number generator that provides in-dependent substreams of pseudo-random digits across processes.The default pseudo-random rumber generator in MCMCpack is the Mersenne twister(Matsumoto and Nishimura1998). MCMCpack provides a default seed,which can be changed in any model-?tting function. The package also provides the generator of L’Ecuyer,Simard,Chen,and Kelton(2002), which provides independent substreams suitable for parallel computation.Again,the seed can be changed by the user,who only needs to provide a unique substream number to en-sure independence.This generator works just?ne with the powerful snow package(Tierney, Rossini,Li,and Sevcikova2008).
To enhance the usability of MCMCpack in the classroom,the package also contains a number of“toy”instructional models.While the use of Monte Carlo simulation is not necessary to characterize the posterior density of these models,we?nd it useful to introduce Monte Carlo methods in contexts where analytical results are readily available.To wit,MCMCpack provides MCmultinomdirichlet(Monte Carlo simulation from a multinomial likelihood with a Dirichlet prior),MCnormalnormal(Monte Carlo simulation from a normal likelihood with known variance and a normal prior),and MCpoissongamma(Monte Carlo simulation from a Poisson likelihood with a gamma prior).
3.An example of Bayesian linear regression
In this section,we look at how to?t a Bayesian linear model using MCMCregress.We use Wilkerson’s(1999)analysis of“killer”amendments in US Congress as an example.A killer amendment causes a bill,that would pass absent the amendment,to fail.According to Wilkerson(1999),a typical example of a killer amendment occurred in1996, ...when the House was asked to approve an international shipbuilding trade agree-
ment that had taken more than a decade to negotiate.During the debate,a major-
ity party member,Herbert Bateman(R-VA),proposed adding30months of“tran-
sition assistance”for US shipbuilders to the agreement....Jennifer Dunn(R-WA)
argued“Some individuals argue that no agreement is better than this agreement.
In reality,if the Bateman amendment is adopted,that is exactly what we would
have:No agreement”(Congressional Record,June131996,H6321)....Ignoring
Dunn’s warning,the House passed the amendment(278-149)and voted for the
bill(325-100).Recognizing that the shipbuilding agreement was dead on arrival
...,the Senate never took it up.(pp.544–545)
Thus,in order to achieve the preferred?nal result(an international shipbuilding trade agree-ment),legislative majorities are often forced to vote against an amendement which seems closer to their preferences(the Bateman amendment)than an original bill.Doing so is a classic example of strategic voting;i.e.,voting for a less preferred choice at an early stage of a voting procedure to ultimately get a better outcome in the end.
Using a linear regression model estimated with ordinary least squares,Wilkerson shows that killer amendments are rare and hence legislative majorities seldom face situations that require
Journal of Statistical Software7 strategic voting.To investigate the existence of strategic voting,Wilkerson relies upon Poole and Rosenthal’s(1997)measures that capture the?t of actual roll call voting positions to estimated unidimensional legislative preferences(NOMINATE,see also Poole,Lewis,Lo, and Carroll2011).Among them,we focus on proportional reduction in error(PRE),which compares NOMINATE’s success at predicting roll call voting positions with the success rate of simply assuming that everyone votes with the majority:
PRE=(Minority vote?NOMINATE classi?cation errors)
(Minority vote)
.
A high value of PRE for a particular roll call implies that the NOMINATE is accurately classifying the votes on that roll call compared to the prediction based on the majority voting. Wilkerson identi?es?ve types of killer amendments and includes them as covariates:major weakening,strengthening,political cover,minor weakening,and new issue.Wilkerson also considers the sponsor.Amendments proposed by a member of the majority party are much more likely to be viewed as sincere attempts to improve the content of the legislation,rather than deliberate e?orts to sabotage https://www.wendangku.net/doc/618588846.html,stly,Wilkerson controls for the chamber of intro-duction,in case there are chamber e?ects on sincerity of voting.The model to be estimated is:
y i=x iβ+εiεi~N(0,σ2)
With priors:
b~N K(b0,B?10)σ2~IG(c0/2,d0/2)
To simulate from the posterior distribution of the model using MCMCregress we load Wilker-son’s data and call MCMCregress.
R>library("MCMCpack")
R>load("killamdt.rda")
R>wilkerson<-MCMCregress(APRE1~STRENGTH+COVER+WEAKMIN+
+NEWISSUE+SPONSOR+CHAMBER,data=killamdt,b0=0,
+B0=0.1,c0=2,d0=0.11,marginal.likelihood="Chib95")
Note that the output from the call to MCMCregress has been placed in an mcmc object called wilkerson.By default,MCMCregress uses a noninformative prior for the coe?cient param-eters.To facilitate model comparison,here we use a weakly informative prior.The prior for βis speci?ed using the b0and B0arguments,where b0is the prior mean and B0is the prior precision ofβ.The prior forσ2is governed by two parameters c0and d0.We use a weakly informative prior(c0=2and d0=0.11)with a mean equal to the marginal variance of the APRE1(0.11).
Before examining the results of our model,it is necessary to examine some diagnostics to assess whether the Markov chain has converged to its stationary distribution.The plot method for an mcmc object produces a trace plot and a density plot for each parameter.To save space,we only report these plots for?rst two parameters in Figure1.There are also a number of standard convergence diagnostics implemented in coda one should check before using the posterior density sample for inference.
The coda summary method can be used to summarize an the mcmc object produced by MCMCregress.
8MCMCpack:Markov Chain Monte Carlo in R
Figure1:Trace and density plots for two coe?cients in the MCMCregress example. R>summary(wilkerson)
Iterations=1001:11000
Thinning interval=1
Number of chains=1
Sample size per chain=10000
1.Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
(Intercept)0.726310.069340.00069340.0007085
STRENGTH-0.095730.158280.00158280.0015134
COVER-0.478220.151220.00151220.0015871
WEAKMIN-0.282460.079890.00079890.0008861
NEWISSUE-0.496850.110020.00110020.0009653
SPONSOR-0.149880.074250.00074250.0007096
CHAMBER0.189240.076660.00076660.0007562
Journal of Statistical Software9
Mean St.Dev Lower Upper
Constant0.72630.06930.59020.8640
STRENGTH?0.09570.1583?0.40370.2115
COVER?0.47820.1512?0.7734?0.1847
WEAKMIN?0.28250.0799?0.4395?0.1274
NEWISSUE?0.49680.1100?0.7128?0.2811
SPONSOR?0.14990.0743?0.2957?0.0050
CHAMBER0.18920.07670.03840.3422
σ20.07490.01390.05230.1069
Table1:Analysis of killer amendments(Wilkerson1999).The dependent variable is propor-tional reduction in error(PRE)when the minority vote is used as benchmark.Lower and Upper indicate central95percent Bayesian credible intervals.
sigma20.074930.013870.00013870.0001681
2.Quantiles for each variable:
2.5%25%50%75%97.5%
(Intercept)0.590170.680420.725960.772070.864040
STRENGTH-0.40370-0.20210-0.095700.011910.211530
COVER-0.77342-0.57890-0.48005-0.37593-0.184722
WEAKMIN-0.43952-0.33595-0.28221-0.22865-0.127368
NEWISSUE-0.71283-0.57028-0.49699-0.42332-0.281150
SPONSOR-0.29571-0.19903-0.14914-0.10004-0.005032
CHAMBER0.038410.138040.189120.240310.342175
sigma20.052320.065120.073430.082940.106865
One of main advantage of Bayesian methods is that posterior inference is quite straightfor-ward using MCMC output since the methods provide direct samples of parameters of inter-est.In this example,Wilkerson’s primary quantity of interest is in whether PRE is higher for minority-sponsored major weakening amendments than for majority-sponsored major weaken-ing amendments.Wilkerson considers a higher PRE for minority-sponsored major weakening amendments than for majority-sponsored major weakening amendments as major evidence for killer amendments.
We can answer this question using our posterior samples.In our data set,SPONSOR is coded 1if majority-sponsored amendments and WEAKMIN is coded1if minor weakening amend-ments,and0otherwise.Thus,the probability that PRE is higher for minority-sponsored major weakening amendments than for majority-sponsored major weakening amendments can be shown by comparing a posterior distribution of Constant with a posterior distribution of Constant+SPONSOR.
Figure2demonstrates that minority-sponsored major weakening amendments are better pre-dicted by NOMINATE scores than majority-sponsored major weakening amendments.This result is consistent with Wilkerson’s?nding that members of the majority party are more likely to give serious consideration to the substance and strategic implications of amendments
10MCMCpack :Markov Chain Monte Carlo in R
0.20.40.6
0.8 1.0
123456PRE D e n s i t y
majority sponsored major weakening
minority sponsored major weakening
Figure 2:PRE for various types of amendments and sponsors.
o?ered by fellow party members.
MCMCpack also provides a tool for Bayesian model comparison (for some of its models):Bayes factors.Bayes factors can be used to compare any model of the same data.Suppose that the observed data y could have been generated under one of two models A 1and A 2.A natural thing to ask from the Bayesian perspective is:“What is the posterior probability that A 1is true assuming either A 1or A 2is true?”Using Bayes theorem we can write:
Pr(A 1|y )=p (y |A 1)Pr(A 1)
p (y |A 1)Pr(A 1)+p (y |A 2)Pr(A 2)
It is instructive to look at the posterior odds in favor of one model (say A 1):
Pr(A 1|y )Pr(A 2|y )=p (y |A 1)p (y |A 2)×Pr(A 1)Pr(A 2)
.What this means is that if we want to move from the prior odds in favor of A 1to the posterior odds in favor of A 1we simply multiply the prior odds by:
B 12=p (y |A 1)
p (y |A 2)
which is called the Bayes factor for A 1relative to A 2.B 12provides a measure of the strength of evidence for A 1over A 2.For a thorough discussion of Bayes factors see Kass and Raftery (1995).
Journal of Statistical Software11
Model1Model2Model3Wilkerson
Model10.00?3.11?1.39?0.52
Model2 3.100.00 1.71 2.58
Model3 1.39?1.720.000.86
Wilkerson0.52?2.58?0.860.00
Table2:A matrix of Bayes factors(on the natural log scale)for the Wilkerson(1999)repli-cation.The dependent variable is PRE.
Mean St.Dev Lower Upper
(Intercept)0.78880.05910.67200.9042
STRENGTH?0.26810.1572?0.57650.0385
COVER?0.47410.1583?0.7849?0.1676
WEAKMIN?0.30590.0839?0.4687?0.1411
NEWISSUE?0.45180.1060?0.6607?0.2404
σ20.08480.01570.05930.1210
Table3:Killer amendments results from model2.Lower and Upper indicate central95 percent Bayesian credible intervals.
The key ingredients necessary for calculating Bayes factors(and posterior model probabilities) are the marginal likelihoods.There are several e?cient methods(Tierney and Kadane1986; Chib1995;Chib and Jeliazkov2001)for these calculations.MCMCregress computes marginal likelihoods using either the Laplace approximation(Tierney and Kadane1986)or Chib’s (1995)method to compute the marginal https://www.wendangku.net/doc/618588846.html,ing Chib’s method,we compare three alternative models with Wilkerson’s model using comparable prior distributions.
R>model1<-MCMCregress(APRE1~STRENGTH+COVER,data=killamdt,
+mcmc=10000,b0=0,B0=0.1,c0=2,d0=0.11,
+marginal.likelihood="Chib95")
R>model2<-MCMCregress(APRE1~STRENGTH+COVER+WEAKMIN+
+NEWISSUE,data=killamdt,mcmc=10000,b0=0,
+B0=0.1,c0=2,d0=0.11,marginal.likelihood="Chib95")
R>model3<-MCMCregress(APRE1~SPONSOR+CHAMBER,data=killamdt,
+mcmc=10000,b0=0,B0=0.1,c0=2,d0=0.11,
+marginal.likelihood="Chib95")
R>BF<-BayesFactor(model1,model2,model3,wilkerson)
R>summary(BF)
The output from BayesFactor shows a matrix of Bayes Factors for all pairs of models speci?ed in the argument.Table2shows that there is positive evidence or better to support model2 over all other models considered.That is,information about the sponsorship and the chamber of introduction does not help explain the variation of PRE given the information on?ve types of killer amendments.We report a posterior summary of model2in Table3for interested readers.
12MCMCpack:Markov Chain Monte Carlo in R
Framework Cycle phases
Goldstein Stagnation Rebirth Expansion War
1816183118491861
1873188518941911
1921193419411969
Level of con?ict Low con?ict High con?ict High con?ict Highest con?ict Wallerstein Victory Maturity Decline Ascent
1816185018731897
191319451967
Level of con?ict Scattered con?ict Lowest con?ict Competitiveness Highest con?ict Modelski and World Power Delegitimation Deconcentration Global war Thompson1816185018741915
19461974
Level of con?ict Lowest con?ict High con?ict Low con?ict Highest con?ict Gilpin UK hegemony No hegemony US hegemony
181518731945
Level of con?ict Low con?ict High con?ict Low con?ict
Table4:Starting dates of cycle phases and predicted level of con?ict.Source:Pollins(1996).
4.A Bayesian Poisson changepoint model example
We show how to?t a Poisson change-point model using MCMCpoissonChange in this section. Generally speaking,a change-point model is a time series model in which data are generated from multiple regimes or states.Main quantities of interests in the change-point analysis are the number,timing,and e?ects of regime changes.Chib(1998)reinterprets the change-point problem in the context of Markov mixture model in which observations are assumed to be drawn from latent states.The sampling algorithm of latent states in MCMCpoissonChange is based on Chib(1998).
Using conjugate prior distributions,we can write a Bayesian Poisson change-point model without covariates as follows:
y t~P oisson(λi),i=1,...,k
λi~G(c0,d0),i=1,...,k
p ii~B eta(α,β),i=1,...,k.
When there is no covariate,all parameters can be sampled using the Gibbs sampler as shown in Chib(1998).When users have covariates to model the expected count,MCMCpoissonChange uses the auxiliary mixture sampling method suggested by Fr¨u hwirth-Schnatter and Wagner (2006)to sample regression parameters.The details of the algorithm are discussed in Park (2010).
In this section,we employ MCMCpoissonChange to revisit the empirical test of cyclical theories of great-power war by Pollins(1996).Pollins identi?es four core frameworks in the study of cyclical patterns of international con?icts based on“(1)the relative importance of global political versus economic subsystems,(2)the dating of the phases they identify in the cycle, and(3)predictions regarding the level of con?ict within particular phases of their respective
Journal of Statistical Software 13
Time m i d a
18501900
19502000
01020304050607
Figure 3:Militarized interstate disputes data.
Mod.1
Mod.2Mod.3Mod.4Mod.5Mod.6Mod.7Mod.8Mod.10.0
?60.5?80.0?93.3?96.0?111.61?109.02?106.46Mod.260.5
0.0?19.5?32.9?35.6?51.12?48.53?45.96Mod.380.0
19.50.0?13.3?16.0?31.59?29.00?26.44Mod.493.3
32.913.30.0?2.7?18.26?15.67?13.11Mod.596.0
35.616.0 2.70.0?15.56?12.97?10.41Mod.6111.6
51.131.618.315.60.00 2.59 5.15Mod.7109.0
48.529.015.713.0?2.590.00 2.56Mod.8
106.546.026.413.110.4?5.15?2.560.00Table 5:A matrix of Bayes factors (on the natural log scale)for the Poisson change-point analysis of militarized interstate disputes between 1816and 1995.
cycles”(Pollins 1996,106).These four core frameworks are a theory of hegemony by Gilpin (1981),a theory of the world system by Wallerstein (1983),the long wave in the global economy by Goldstein (1988),and the world leadership cycle by Modelski and Thompson (1988).
These four theories give di?erent weights on the role of global economic order and global political order in explaining the patterns of international con?icts.For example,while the theory of Long Wave and the theory of world system argue that the growth and stagnation of global economy (Goldstein 1988)or economic hegemon in the system (Wallerstein 1983)generate di?erent levels of international con?icts,Modelski and Thompson (1988)and Gilpin (1981)consider the rise and decline of global leadership or a hegemonic power to play a central role in producing international con?icts.These four theories are further diverged in their predictions on the speci?c cycle phases of international con?icts as shown in Table 4.Following Pollins (1996),we use the annual number of military con?icts from the militarized interstate disputes (MID)data (Jones,Bremer,and Singer 1996)to detect the number and
14MCMCpack :Markov Chain Monte Carlo in R
Posterior Regime Probability
Time P r (S t = k |Y t )18501900
19502000
0.00.2
0.40.60.8
1.0q q q q q q q q q q q q q q q q q q q q q q q q q q
q
q q q
q
q
q
q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q
q q q q
q
q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q
q q q q q q q q q q q q q q q q q q q q q q q q q
q q q q q q q q q
q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q
q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q
q
q
q
q
q q State1State2State3State4State5State6State7Figure 4:Posterior probabilities of states:each circled line shows the time-varying probabil-ities of each regime.For example,the probability of Regime 1(the solid circled line)starts from 1in the beginning of the sample period and drops to 0around the 1840s.Estimates are from the Poisson change-point analysis of militarized interstate disputes from 1816to 1995.timing of cycle phases in international con?icts.
To implement MCMCpoissonChange ,a user needs to specify data,the number of states (m ),and hyperparameters of λ(c0and d0).Users must provide hyperparameters of λ(c0and d0)for model comparison using Bayes factors.Marginal likelihoods are not identi?ed under di?use priors.Thus,users should avoid di?use priors for model comparison.In this speci?c example,we set c0=13and d0=1since the mean of MID is about 13.The following R commands implement eight Poisson change-point models for the militarized interstate disputes data.R>
model1<-MCMCpoissonChange(mida,m =1,c0=13,d0=1,+
marginal.likelihood =c("Chib95"))R>
model2<-MCMCpoissonChange(mida,m =2,c0=13,d0=1,+
marginal.likelihood =c("Chib95"))R>
model3<-MCMCpoissonChange(mida,m =3,c0=13,d0=1,+
marginal.likelihood =c("Chib95"))R>
model4<-MCMCpoissonChange(mida,m =4,c0=13,d0=1,+
marginal.likelihood =c("Chib95"))R>
model5<-MCMCpoissonChange(mida,m =5,c0=13,d0=1,+
marginal.likelihood =c("Chib95"))R>
model6<-MCMCpoissonChange(mida,m =6,c0=13,d0=1,+
marginal.likelihood =c("Chib95"))R>
model7<-MCMCpoissonChange(mida,m =7,c0=13,d0=1,
+marginal.likelihood =c("Chib95"))
Journal of Statistical Software 15
18501900195020000.000.15
1850190019502000
0.00.6
18501900195020000.00.4
1850190019502000
0.00.61850190019502000
0.00.41850190019502000
0.00.2
Figure 5:Posterior density of regime change probability:each panel shows the posterior distribution of regime change probabilities.For example,the top panel shows the probability density of the regime change from Regime 1to Regime 2.
R>model8<-MCMCpoissonChange(mida,m =8,c0=13,d0=1,
+marginal.likelihood =c("Chib95"))
The results of the model comparison are summarized in Table 5.The ratio of marginal likeli-hoods shows that the sixth change-point model is favored over the https://www.wendangku.net/doc/618588846.html,ing Je?rey’s rule,which provides guidelines as to the strength of a Bayes factor,we can conclude that there is positive evidence for the sixth change-point model over the alternatives considered.
plotState()generates a plot of posterior state probabilities and plotChangepoint()plots of posterior probabilities of six change-points as shown in Figures 4and 5,respectively.Both graphs indicate that critical shifts in the level of international con?icts are concentrated in the ?rst half of the twentieth century and cyclical patterns predicted by the four theoretical models are not supported by the evidence.However,the picture provided by the Poisson change-point analysis is closest to Gilpin’s theory of hegemony which predicts high levels of con?icts in the absence of a world hegemon between 1873and 1945.
16MCMCpack:Markov Chain Monte Carlo in R
5.Conclusion
MCMCpack remains a work in progress.While it is di?cult to commit to speci?c future development e?orts,we hope to implement additional models others will?nd useful in their research and teaching.The ability to estimate marginal likelihoods and compute Bayes factors is a new feature in MCMCpack,and is only currently available for a handful of models.We hope to implement this for a wider range of models in the future.Of course,since MCMCpack is open source,others can build on this code-base to implement models of their choice.
In the next few years Jong Hee Park will develop a number of models to be used for Bayesian time series analysis,including change-point models and state space models.He intends to expand some of these models to include covariates,and to provide various methods to enhance usability for this class of models.Finally,while the the coda package provides a great deal of power for convergence checking and posterior density summarization,it does not provide adequate methods for formal model assessment,including prior and predictive checks,Bayes factors,and Bayesian model averaging(Hoeting,Madigan,Raftery,and Volinsky1999).Our hope is develop software to make these tools as easy to use as possible.
Acknowledgments
We have received a great deal of support for the MCMCpack project.First,we thank all of those in the R community who have used our software,found bugs,and provided sugges-tions,feedback,and patches.Second,we would like to thank all of the research assistants who have worked with us at Harvard University and Washington University over the years on this project:Lucy Barnes,Matthew Fasman,Ben Goodrich,Steve Haptonstahl,Kate Jensen,Laura Keys,Dan Pemstein,and Ellie Powell.Finally,we gratefully acknowledge ?nancial support from:the US National Science Foundation,Program in Methodology,Mea-surement,and Statistics,Grants SES-0350646and SES-0350613,the Institute for Quantitative Social Sciences(https://www.wendangku.net/doc/618588846.html,/)at Harvard University,the Weidenbaum Center on the Economy,Government,and Public Policy(https://www.wendangku.net/doc/618588846.html,/)at Washington University,and the Center for Empirical Research in the Law(https://www.wendangku.net/doc/618588846.html,/) at Washington University.Neither the National Science Foundation,Washington University, or Harvard University bear any responsibility for the content of this package.In addition, Quinn thanks the Center for Advanced Study in the Behavioral Sciences for its hospitality and support.
References
Albert JH,Chib S(1993).“Bayesian Analysis of Binary and Polychotomous Response Data.”Journal of the American Statistical Association,88,669–679.
Anderson B,Fox J,Frankin C,Gill J(2003).“The Times They R Achanging:The State of Advanced Methodology Curriculum at ICPSR.”The Political Methodologist,11,22–24. Becker RA,Chambers JM,Wilks AR(1988).The New S Language.Chapman&Hall,London. Bernardo J,Smith AFM(1994).Bayesian Theory.John Wiley&Sons,New York.
Journal of Statistical Software17 Besag J,Green P,Higdon D,Mengersen K(1995).“Bayesian Computation and Stochastic Systems.”Statistical Science,10,3–66.
Chib S(1995).“Marginal Likelihood From the Gibbs Output.”Journal of the American Statistical Association,90,1313–1321.
Chib S(1998).“Estimation and Comparison of Multiple Change-Point Models.”Journal of Econometrics,86(2),221–241.
Chib S,Jeliazkov I(2001).“Marginal Likelihood From the Metropolis Hastings Output.”Journal of the American Statistical Association,96,270–281.
Clinton J,Jackman S,Rivers D(2004).“The Statistical Analysis of Roll Call Data.”American Political Science Review,98,355–370.
Cowles MK(1996).“Accelerating Monte Carlo Markov Chain Convergence for Cumulative-Link Generalized Linear Models.”Statistics and Computing,6,101–110.
Fox J(2002).An R and S-PLUS Companion to Applied Regression.Sage Publications, Thousand Oaks,CA.
Fr¨u hwirth-Schnatter S,Wagner H(2006).“Auxiliary Mixture Sampling for Parameter-driven Models of Time Series of Small Counts with Applications to State Space Modelling.”Biometrika,93,827–841.
Gelfand AE,Smith AFM(1990).“Sampling-Based Approaches to Calculating Marginal Den-sities.”Journal of the American Statistical Association,85,398–409.
Gelman A,Carlin JB,Stern HS,Rubin DB(2003).Bayesian Data Analysis.2nd edition. Chapman&Hall/CRC,Boca Raton.
Gelman A,King G(1990).“Estimating the Electoral Consequences of Legislative Redistrict-ing.”Journal of the American Statistical Association,85,274–282.
Geweke J(1999).“Using Simulation Methods for Bayesian Econometric Models:Inference, Development,and Commnication.”Econometric Reviews,18,1–126.
Geweke J,Keane M(2000).“An Empirical Analysis of Income Dynamics Among Men in the PSID:1968-1989.”Journal of Econometrics,96,293–356.
Gill J(2002).Bayesian Methods:A Social and Behavioral Sciences Approach.Chapman& Hall,Boca Raton.
Gilpin R(1981).War and Change in World Politics.Cambridge University Press,Cambridge. Goldstein J(1988).Long Cycles:Prosperity and War in the Modern Age.Yale University Press,New Haven.
Hoeting JA,Madigan D,Raftery AE,Volinsky CT(1999).“Bayesian Model Averaging:A Tutorial.”Statistical Science,14,382–417.
Ho?PD,Raftery AE,Handcock MS(2002).“Latent Space Approaches to Social Network Analysis.”Journal of the American Statistical Association,97,1090–1098.
18MCMCpack:Markov Chain Monte Carlo in R
Ho?PD,Ward MD(2004).“Modeling Dependencies in International Relations Networks.”Political Analysis,12,160–175.
Je?reys SH(1998).Theory of Probability.3rd edition.Clarendon Press,Oxford.
Johnson VE,Albert JH(1999).Ordinal Data Modeling.Springer-Verlag,New York.
Jones DM,Bremer SA,Singer JD(1996).“Militarized Interstate Disputes,1816–1992:Ra-tionale,Coding Rules,and Empirical Patterns.”Con?ict Management and Peace Science, 15,163–213.
Kass RE,Raftery AE(1995).“Bayes Factors.”Journal of the American Statistical Association, 90,773–795.
King G,Rosen O,Tanner MA(1999).“Binomial-Beta Hierarchical Models for Ecological Inference.”Sociological Research and Methods,28,61–90.
L’Ecuyer P,Simard R,Chen EJ,Kelton WD(2002).“An Object-Oriented Random-Number Package with Many Long Streams and Substreams.”Operations Research,50,1073–1075.
Ligges U,Thomas A,Spiegelhalter D,Best N,Lunn D,Rice K,Sturtz S(2009).BRugs0.5: OpenBUGS and Its R/S-PLUS Interface BRugs.URL https://www.wendangku.net/doc/618588846.html,/ pub/RWin/src/contrib/.
Liu JS,Wong WH,Kong A(1994).“Covariance Structure of the Gibbs Sampler with Ap-plications to the Comparisons of Estimators and Augmentation Schemes.”Journal of the American Statistical Association,89,958–966.
Martin AD,Quinn KM(2002).“Dynamic Ideal Point Estimation via Markov Chain Monte Carlo for the U.S.Supreme Court,1953–1999.”Political Analysis,10,134–153.
Martin AD,Quinn KM,Park JH(2011).MCMCpack:Markov Chain Monte Carlo (MCMC)Package.R package version 1.0-10,URL https://www.wendangku.net/doc/618588846.html,/ package=MCMCpack.
Matsumoto M,Nishimura T(1998).“Mersenne Twister:A623-Dimensionally Equidistributed Uniform Pseudo-Random Number Generator.”ACM Transactions on Modeling and Com-puter Simulation,8,3–30.
Modelski G,Thompson WR(1988).Seapower in Global Politics.University of Washington Press,Seattle.
Park JH(2010).“Structural Change in U.S.Presidents’Use of Force.”American Journal of Political Science,54(3),766–782.
Park JH,Martin AD,Quinn KM(2010).“CRAN Task View:Bayesian Inference.”Ver-sion2010-12-14,URL https://www.wendangku.net/doc/618588846.html,/view=Bayesian.
Plummer M(2003).“JAGS:A Program for Analysis of Bayesian Graphical Models Using Gibbs Sampling.”In K Hornik,F Leisch,A Zeileis(eds.),Proceedings of the3rd Interna-tional Workshop on Distributed Statistical Computing,Vienna,Austria.ISSN1609-395X, URL http://www.ci.tuwien.ac.at/Conferences/DSC-2003/Proceedings/.
Journal of Statistical Software19 Plummer M(2005).JAGS:Just Another Gibbs Sampler.Version0.8,URL http://www-fis. iarc.fr/~martyn/software/jags/.
Plummer M(2009).rjags:Interface to the JAGS MCMC Library.R package version1.0.3, URL https://www.wendangku.net/doc/618588846.html,/package=rjags.
Plummer M,Best N,Cowles K,Vines K(2006).“coda:Convergence Diagnosis and Out-put Analysis for MCMC.”R News,6(1),7–11.URL https://www.wendangku.net/doc/618588846.html,/doc/ Rnews/.
Pollins BM(1996).“Global Political Order,Economic Change,and Armed Con?ict:Coe-volving Systems and the Use of Force.”American Political Science Review,90,103–117.
Poole K,Lewis J,Lo J,Carroll R(2011).“Scaling Roll Call Votes with wnominate in R.”Journal of Statistical Software,42(14),1–21.URL https://www.wendangku.net/doc/618588846.html,/v42/i14/.
Poole KT,Rosenthal H(1997).Congress:A Political-Economic History of Roll Call Voting. Oxford University Press,Oxford.
Quinn KM,Martin AD(2002).“An Integrated Computational Model of Electoral Competi-tion.”Statistical Science,17,405–419.
Quinn KM,Martin AD,Whitford AB(1999).“Voter Choice in Multi-Party Democracies: A Test of Competing Theories and Models.”American Journal of Political Science,43, 1231–1247.
R Development Core Team(2011a).An Introduction to R.Vienna,Austria.ISBN3-900051-12-7,URL https://www.wendangku.net/doc/618588846.html,/.
R Development Core Team(2011b).R:A Language and Environment for Statistical Comput-ing.R Foundation for Statistical Computing,Vienna,Austria.ISBN3-900051-07-0,URL https://www.wendangku.net/doc/618588846.html,.
R Development Core Team(2011c).R Language De?nition.Vienna,Austria.ISBN3-900051-13-5,URL https://www.wendangku.net/doc/618588846.html,/.
Robert C,Casella G(2004).Monte Carlo Statistical Methods.2nd edition.Springer-Verlag, New York.
Rossi P,McCulloch R(2006).bayesm:Bayesian Inference for Marketing/Micro-Econometrics.R package version 2.0-8,URL https://www.wendangku.net/doc/618588846.html,/package= bayesm.
Spiegelhalter DJ,Thomas A,Best NG,Lunn D(2003).WinBUGS Version1.4User Manual. MRC Biostatistics Unit,Cambridge.URL https://www.wendangku.net/doc/618588846.html,/bugs/.
Sturtz S,Ligges U,Gelman A(2005).“R2WinBUGS:A Package for Running WinBUGS from R.”Journal of Statistical Software,12(3),1–16.URL https://www.wendangku.net/doc/618588846.html,/ v12/i03/.
Tanner MA,Wong W(1987).“The Calculation of Posterior Distributions by Data Augmen-tation.”Journal of the American Statistical Association,82,528–550.
20MCMCpack:Markov Chain Monte Carlo in R
Thomas A,O’Hara B,Ligges U,Sturtz S(2006).“Making BUGS Open.”R News,6(1), 12–17.URL https://www.wendangku.net/doc/618588846.html,/doc/Rnews/.
Tierney L,Kadane JB(1986).“Accurate Approximation for Posterior Moments and Marginal Densities.”Journal of the American Statistical Association,81,82–86.
Tierney L,Rossini AJ,Li N,Sevcikova H(2008).snow:Simple Network of Workstations. R package version0.2-2,URL https://www.wendangku.net/doc/618588846.html,/package=snow.
Venables WN,Ripley BD(2000).S Programming.Springer-Verlag,New York.
Venables WN,Ripley BD(2002).Modern Applied Statistics with S.4th edition.Springer-Verlag,New York.
Wallerstein I(1983).“Three Instances of Hegemony in the History of the World Economy.”International Jounal of Comparative Sociology,24,100–108.
Western B(1998).“Causal Heterogeneity in Comparative Research:A Bayesian Hierarchical Modelling Approach.”American Journal of Political Science,42,1233–1259.
Wilkerson JD(1999).“‘Killer’Amendments in Congress.”American Political Science Review, 93,535–552.
A?liation:
Andrew D.Martin
Washington University School of Law
Campus Box1120
One Brookings Drive
St.Louis,MO63130,United States of America
E-mail:admartin@https://www.wendangku.net/doc/618588846.html,
URL:https://www.wendangku.net/doc/618588846.html,/
Kevin M.Quinn
UC Berkeley School of Law
490Simon#7200
University of California,Berkeley
Berkeley,CA94720-7200,United States of America
E-mail:kquinn@https://www.wendangku.net/doc/618588846.html,
URL:https://www.wendangku.net/doc/618588846.html,/kevinmquinn.htm
Jong Hee Park
Department of Political Science