2.1 Basic Bayesian Theory

The Bayesian process of data analysis is characterized by three primary attributes: (1) a willingness to assign prior distributions to unknown parameters, (2) the use of Bayes’ rule to obtain a posterior distribution for unknown parameters and missing data conditioned on observable data, and (3) the description of inferences in probabilistic terms. The core philosophical foundation of Bayesian inference is the consideration of parameters as random quantities. A primary advantage of this approach is that there is no restriction to building complex models with multiple levels and many unknown parameters. Because model assumptions are much more overt in the Bayesian setup, readers can more accurately assess model quality.

The primary goal in the Bayesian approach is to characterize the posterior density of each of the parameters of interest through the computation and the marginalization of the joint posterior distribution of the unknown parameter vector, given the observed data and the model specification, $\pi (\mathit{\theta }|\mathbf{\text{y}})$, for vectors $\mathit{\theta }$ and data $\mathbf{\text{y}}$. Here the use of “ $\pi$” is a reminder that this is a posterior (post-data) probability distribution. The posterior distribution of the coefficients of interest is obtained by multiplying the prior by the likelihood function using proportionality:

$\begin{array}{r}\pi (\mathit{\theta }|\mathbf{\text{y}})\propto p(\mathit{\theta })L(\mathit{\theta }|\mathbf{\text{y}}),\end{array}$(2.1)

meaning that the posterior distribution is a compromise between prior information and data information. The fundamental principle and process of Bayesian inference starts with an unconditional prior distribution, $p(\mathit{\theta })$, that reflects prior belief about the parameters of interest, and multiplies it with the likelihood function, $L(\mathit{\theta }|\mathbf{\text{y}})$, bringing in the data information to produce the posterior distribution, $\pi (\mathit{\theta }|\mathbf{\text{y}})$. The posterior is now an updated distribution that reflects the most recent knowledge about some phenomenon of interest. Additionally, the result is a distribution, not a point estimate, so uncertainty is built into the results as distributional variability.

This idea of compromise is critical to Bayesian thinking. When there is very reliable prior information and sketchy data information, then we should rely more on the specification of the prior. On the other hand, when there is a large amount of reliable data, then we should favor this type of information over any prior. A key feature of Bayesian inference is that this compromise is mathematically built into the process for any Bayesian result. So we get an important balance based on the amount of contribution from each of the two information sources. In addition, a posterior distribution can be treated as a prior distribution when another set of data is observed, and this process can be repeated as many times as new data arrive to obtain new posteriors that are conditional on all previous sets of observed data, as if they had arrived at the same time. This is a process called “Bayesian updating” or “Bayesian learning.”

The Bayesian updating process is illustrated in Figure 1 where a gamma distribution prior is applied to a Poisson likelihood function with synthetic data ($n=20$) generated by the R command rpois(20,5) or the Python command np.random.poisson(5,20), for a single parameter $\theta$. This is a conjugate setup, meaning that the form of the prior flows down to the posterior, albeit with different parameterization. Conjugacy is a special, and convenient, case of Bayesian modeling that has simple mathematical properties for determining the posterior distribution. In this case the prior and posterior parameterizations are

$\begin{array}{rl}p(\theta |4,1)& =1^4\mathrm{\Gamma }(4{)}^{-1}{\theta }^{4-1}\text{exp}[-1\times \theta ]\\ \pi (\theta |109,21)& ={21}^{109}\mathrm{\Gamma }(109{)}^{-21}{\theta }^{109-1}\text{exp}[-21\times \theta ],\end{array}$

using the rate version of the gamma PDF:

$\mathcal{G}(x|\alpha ,\beta )=\frac{{\beta }^{\alpha }}{\mathrm{\Gamma }(\alpha )}{x}^{\alpha -1}\text{exp}[-x\beta ],0\le x<\mathrm{\infty },0\le \alpha ,\beta .$(2.2)

Suppose the threshold at 4 on the $x$-axis is substantively important. In the prior distribution, 0.566 of the density is to the left of 4, whereas 0.005 of the posterior density is to the left of 4. So, even though there are only 20 data values, the likelihood function changes the prior substantially.

Figure 1 Poisson–Gamma update.

In many Bayesian models, priors themselves have prior parameters (hyperpriors), so there can be additional terms in the product but no substantial increase in complexity:

$\pi (\theta |\mathbf{\text{y}})\propto p(\theta |\gamma )p(\gamma )L(\theta |\mathbf{\text{y}}),$(2.3)

where $\gamma$ is a higher-level parametric assignment with prior distribution $p(\gamma )$. Proportionality as used in (2.3) is typically employed to simplify the notation, since in the Bayesian setup, $p(\mathbf{\text{y}})$ is observed and therefore is assumed fixed and independent of $p(\theta )$. In addition, Bayesian model specification can be naturally hierarchical, since parameters and hyperparameters can be conduits to additional specifications at group levels.

2.2 Interval Estimation

Highest posterior density (HPD) regions are the Bayesian equivalent of confidence intervals in which the region describes the area over the coefficient’s support with $(1-\alpha )$% of the density with the highest probability. Like frequentist confidence intervals, an HPD region that does not contain zero implies that the coefficient estimate is deemed to be reliable, but instead of being $(1-\alpha )$% “confident” that the interval covers the true parameter value over multiple replications of the exact same experiment, an HPD interval provides a $(1-\alpha )$% probability that the true parameter is in the interval. Sometimes HPD regions are relatively wide, indicating that the coefficient varies considerably and is less likely to be reliable, and sometimes these regions are quite narrow, indicating greater certainty about the central location of the parameter. Highest posterior density regions do not have to be contiguous, and contiguous HPDs (no breaks on the $x$-axis) are usually called credible intervals.

2.3 The Bayes Factor

Formal hypothesis testing can also be easily performed in Bayesian inference. Start with ${\mathrm{\Theta }}_0$ and ${\mathrm{\Theta }}_1$ representing competing hypotheses about the location of some $\theta$ parameter to be estimated. These two statements about $\theta$ form a partition of the sample space, meaning that $\mathrm{\Theta }={\mathrm{\Theta }}_0\cup {\mathrm{\Theta }}_1,{\mathrm{\Theta }}_0\cap {\mathrm{\Theta }}_1=\mathrm{\varnothing }$, where “ $\mathrm{\varnothing }$” denotes the empty or null set. Now prior probabilities are assigned to each of the two outcomes: $p_0=p(\theta \in {\mathrm{\Theta }}_0)$ and $p_1=p(\theta \in {\mathrm{\Theta }}_1)$. Since the priors are different (but the form of the likelihood function is the same), there are two competing posterior distributions: ${\pi }_0=p(\theta \in {\mathrm{\Theta }}_0|\mathbf{\text{X}})$ and ${\pi }_1=p(\theta \in {\mathrm{\Theta }}_1|\mathbf{\text{X}})$. We can now define the prior odds, $p_0/p_1$, and the posterior odds, ${\pi }_0/{\pi }_1$. This leads to a prior-to-posterior ratio of $(p_0/p_1)/({\pi }_0/{\pi }_1)$, which is called the Bayes factor, interpreted as odds favoring $H_0$ versus $H_1$ given the observed data, and it leads naturally to Bayesian hypothesis testing. Another interpretation is that the Bayes factor is the ratio of improvement in odds by accounting for the data.

More formally, start with the two competing models, $M_0:f_0(\mathbf{\text{y}}|{\theta }_0)$ and $M_1:f_1(\mathbf{\text{y}}|{\theta }_1)$. Specify corresponding parameter priors: $p_0({\theta }_0)$, $p_1({\theta }_1)$, and model priors: $p(M_0)$ and $p(M_1)$. Note that $p(\mathbf{\text{y}}|M_i)={\int }_{{\theta }_i}{f}_i(\mathbf{\text{x}}|{\theta }_i)p_i({\theta }_i)d{\theta }_i$. Therefore:

$\begin{array}{rl}\underset{\text{posterior odds}}{\underbrace{\frac{p(M_0|\mathbf{\text{y}})}{p(M_1|\mathbf{\text{y}})}}}& =\underset{\text{prior odds/data}}{\underbrace{\frac{p(M_0)/p(\mathbf{\text{y}})}{p(M_1)/p(\mathbf{\text{y}})}}}\times \underset{\text{Bayes factor}}{\underbrace{\frac{{\int }_{{\theta }_0}{f}_0(\mathbf{\text{y}}|{\theta }_0)p_0({\theta }_0)d{\theta }_0}{{\int }_{{\theta }_1}{f}_1(\mathbf{\text{y}}|{\theta }_1)p_1({\theta }_1)d{\theta }_1}}}.\\ \text{posterior odds ratio}& =\text{prior odds ratio}\times \text{integrated likelihood ratio.}\end{array}$(2.4)

Now rearranging and canceling out $p(\mathbf{\text{y}})$ gives:

$\begin{array}{rl}B(\mathbf{\text{y}})& =\frac{p(M_0|\mathbf{\text{y}})/p(M_0)}{p(M_1|\mathbf{\text{y}})/p(M_1)}\\ & =\text{``posterior to prior odds ratio.''}\end{array}$(2.5)

The Bayes factor here is analogous to a likelihood ratio test, except that it includes prior information. Also, the tested hypotheses do not need to be “nested” as does the likelihood ratio (one model is a subset of the other). Regretfully, though, there is no distributional test for the Bayes factor. Since a result near one is obviously inconclusive, what we want is a number that is decisively large, like 10 or greater, or decisively small, like 0.1 or less. There are some published thresholds contained in Jeffreys (Reference Jeffreys1998) and Kass and Raftery (Reference Kass and Raftery1995), but they are completely arbitrary.

It is important to know that the Bayes factor as evidence is completely different than using p-values (or similar devices) in non-Bayesian inference. For example, consider testing whether a coin is fair or biased at a specific value: $H_0:\theta =0.5$ versus $H_1:\theta =0.7$. Assign equal priors: $p(H_0)=p(H_1)=\frac{1}{2}$, so that the prior odds ratio is just one, which is commonly done. We flip the coin 10 times and get 7 heads, which is exactly at the alternative hypothesis. The Bayes factor between “model 0 supported” and “minimal evidence against model 0” is calculated by

$B(\mathbf{\text{y}})=\frac{p(M_0|\mathbf{\text{y}})/p(M_0)}{p(M_1|\mathbf{\text{y}})/p(M_1)}=\frac{p(y=7|\theta =0.5)}{p(y=7|\theta =0.7)}=\frac{\left(\genfrac{}{}{0em}{}{10}{7}\right)(0.5{)}^7(0.5{)}^3}{\left(\genfrac{}{}{0em}{}{10}{7}\right)(0.7{)}^7(0.3{)}^3}=0.44,$

which provides no strong evidence for either hypothesis. Now calculate the p-value in the standard way:

$\begin{array}{rl}p& =p(\text{observed data or more extreme}|H_0)=p(x=\{7,8,9,10\}|\theta =0.5)\\ & =\sum _{x=7}^{10}\left(\genfrac{}{}{0em}{}{10}{x}\right)(0.5{)}^x(0.5{)}^{10-x}=0.172,\text{fail to reject}{H}_0.\end{array}$

Notice that the manner in which we make these conclusions is very different. The surprising result is that we cannot reject $H_0:\theta =0.5$ with the p-value as evidence even though the data landed right on the alternative of $H_1:\theta =0.7$. This is because the two hypotheses are relatively close and p-values are poor tools for asserting scientific evidence. Bayes factors are very compelling because they provide a direct way to compare competing hypotheses that include both prior and data information, but there are some shortcomings. First, they are known to be sensitive to the choice of priors, although an overwhelming number of authors set the prior probabilities to 0.5 each. Second, there is no ability to have directional hypotheses in the conventional sense. However, direction can be built into the hypothesis structure by specifying regions such as $\theta <0$ versus $\theta \ge 0$. Finally, improper priors (those that do not sum or integrate to a finite value) do not work because the infinity values do not cancel out in the ratio, although there are several rather tortured work-arounds; see Gill (Reference Gill2014), chapter 6.

2.4 Using Simulation

The challenge before us now is to use these Bayesian principles to develop regression-style specification of interest to social scientists. However, many of the important and popular forms require Markov chain Monte Carlo (MCMC) so that the computer performs the difficult or impossible tasks for us in estimating these forms and getting marginal posterior distributions for each of the unknown parameters. One major impediment to wide usage of Bayesian regression specifications is the hurdle of learning MCMC. We break it down into three stages here: (1) understanding the general principle of Monte Carlo methods for making the computer do a lot of repetitive tasks, (2) defining what a Markov chain does mathematically in an accessible way, and (3) combining these two principles to produce a powerful tool that lets Bayesians model nearly anything. Our objective in the next three sections is to demystify MCMC tools and make them widely usable to a large class of researchers. We do this by introducing Monte Carlo simulation in the next section, then in Section 4 we provide a gentle introduction to the mathematics of Markov chains and then combine these principles to produce the powerful MCMC estimation engine in Section 5.

3.1 Random Number Generation with Uses

Surprisingly, all random number generation starts with producing $[0:1]$ continuous uniform random variables within computer algorithms. Then there is a large set of transformations that provide any distributional form of interest. These uniform random numbers are actually not random at all. They are pseudo-random numbers that are produced deterministically by the computer with a defined period before repeating the sequence identically. Fortunately, due to work mostly done in the 1980s, the length of this period is so large that we do not have to care. We can generate different sequences by setting a “seed” to start a new process. This is illustrated in the following code boxes in R and Python where we set this random seed before generating five $\mathcal{U}(0:1)$ random numbers, change the seed to produce a different set, and then change the seed back to get the original set again.

Box 1R Code for Random Uniforms

Code 1.01

Code 1.01

 set.seed(1) 
 runif(5) 
 [1] 0.2655087 0.3721239 0.5728534 0.9082078 0.2016819 
 set.seed(2) 
 runif(5) 
 [1] 0.0693609 0.8177752 0.9426217 0.2693818 0.1693481 
 set.seed(1) 
 runif(5) 
 [1] 0.2655087 0.3721239 0.5728534 0.9082078 0.2016819

Box 2Python Code for Random Uniforms

Code 1.02

Code 1.02

 import numpy as np 
 np.random.seed(1) 
 np.random.uniform(low = 0.0, high = 1.0, size = 5) 
 array([4.17022005e-01, 7.20324493e-01, 1.14374817e-04, 
 3.02332573e-01, 1.46755891e-01]) 
 np.random.seed(2) 
 np.random.uniform(low = 0.0, high = 1.0, size = 5) 
 array([0.4359949 , 0.02592623, 0.54966248, 0.43532239, 
 0.4203678 ])

The Monte Carlo principle is actually much more general than this, even though generating random numbers according to some distribution is very important. Suppose we could generate random numbers in a way that does work for us, and often you can get the right answer faster with simulation than by analytics. For instance, what is the 90 percent confidence interval for a random variable distributed $N(6,3)$? This is obviously an easy analytical calculation but we can also do this quickly with two lines of code as shown in what follows (one line if we want to be clever). We simply generate 1,000,000 values from the desired normal distribution, sort them, and find the desired interval end-points empirically as if these were actual data. The large number of generated values (it can be larger) assures us a high level of accuracy.

Box 3R Code for a Confidence Interval

Code 1.03

Code 1.03

 n.data <- rnorm(1000000,6,3) 
 sort(n.data)[c(50000,950000)] 
 [1] 1.063709 10.944484 
 mean(n.data) + c(-1,1)*qnorm(0.95)*sd(n.data) 
 [1] 1.064769 10.943997

Box 4Python Code for a Confidence Interval

Code 1.04

Code 1.04

 import numpy as np 
 import scipy.stats 
 n_data = np.random.normal(loc=6, scale=3, size=1000000) 
 [sorted(n_data)[50000], sorted(n_data)[950000] ] 
 [1.0671053614550638, 10.926188882182167] 
 [np.mean(n_data) 
 + scipy.stats.norm.ppf(0.05)*np.std(n_data), 
 np.mean(n_data) 
 + scipy.stats.norm.ppf(0.95)*np.std(n_data)] 
 [1.066819319906851, 10.927139948501146]

Here R and Python have slightly different answers because the default random number generators are different. Notice also that the two methods also have slightly different answers and this is due to rounding in the different calculations. There is also an error term introduced by using Monte Carlo methods instead of an analytical calculation: simulation error. The good news, however, is that this type of error decreases at a rate of $1/\sqrt{n}$, where $n$ is the number of simulated values used. This means that the user controls the simulation error and the only cost is additional user time on increasingly fast computers.

3.2 Basic Monte Carlo Integration

Another use of Monte Carlo simulation is doing integral calculus that is hard or impossible analytically. Suppose $f(y)$ is a function that is difficult to integrate over a specific range but we can draw values from this form. For example, the expected value of $y$ distributed $f(y)$ restricted to being between $a$ and $b$ is given by the integral:

$\begin{array}{r}I[a,b]={\int }_a^{b}yf(y)dy.\end{array}$(3.1)

To get this integral by simulation, first generate a large number ($n$) of simulated values of $y$ from $f(y)$. Then we discard any values below $a$ and above $b$ to get the sample $n_{trunc}$ (e.g. left after truncation). The estimate of the integral is merely

$\begin{array}{r}\hat{I}[a,b]=\frac{1}{n_{trunc}}\sum _{n_{trunc}}{y}_i.\end{array}$(3.2)

So we substitute a mean of the empirical sample for the analytical integral calculation. From the Strong Law of Large Numbers, the calculation of $\hat{I}[a,b]$ converges with probability one to the actual true answer. As in the preceding example the simulation error of $\hat{I}[a,b]$ goes to zero proportional to $1/\sqrt{n}$. In this case the expression of the variance of the simulation estimate is

$\text{Var}(\hat{I}[a,b])=\frac{1}{n_{trunc}(n_{trunc}-1)}\sum _{n_{trunc}}(y_{i,trunc}-\hat{I}[a,b]{)}^2$(3.3)

Note that the researcher fully controls the simulation size, $n$, and therefore the simulation accuracy of the estimate.

As a specific example, suppose we want to integrate a $\mathcal{N}(15,3)$ over the range $[12:13]$:

$I[12,13]={\int }_{12}^{13}(2\pi {\sigma }^2{)}^{-\frac{1}{2}}\text{exp}[-\frac{1}{{\sigma }^2}(y-\bar{y}{)}^2]dy.$

This is easily calculated by the following:

$\hat{I}[12,13]=\frac{1}{n_{trunc}}\sum _{n_{trunc}}{y}_i.$

The calculation is given in the following R and Python code boxes.

Box 5R Code for a Normal Integral

Code 1.05

Code 1.05

 norm.sample <- rnorm(mean=15,sd=3,1000000) 
 mean(norm.sample[norm.sample>12 & norm.sample <13]) 
 [1] 12.52342

Box 6Python Code for a Normal Integral

Code 1.06

Code 1.06

 norm_sample = rng.normal(loc=15, scale=3, size=1000000) 
 np.mean(norm_sample[(norm_sample >= 12) 
 & (norm_sample <= 13)]) 
 12.524258697969533

In a Bayesian context using 3.1 we can replace arbitrary $f(y)$ with a posterior statement $\pi (\theta |y)$ and $y$ in the preceding with $f(\theta )$. So $I[a,b]$ is the posterior expectation of $\theta$, given by:

$|\mathbb{E}|[\theta |y]=\int \pi (\theta |x)f(\theta )d\theta \approx \frac{1}{n}\sum _{i=1}^{n}f({\theta }_i).$(3.4)

We can see here that if it is possible to generate values from $f(\theta )$, then we can calculate this expected value or any other related statistic. This is actually the last step of the MCMC process, which we will see in Section 4.2 on Markov Chain Theory.

3.3 Rejection Sampling

Suppose that we lack the ability to produce samples from some known distribution but we still have a difficult analytical problem that suggests Monte Carlo simulation. In some circumstances we can then use rejection sampling to integrate as done earlier. Rejection sampling is actually used in two primary contexts: obtaining dimensional integral quantities, and generating random variables. In addition, rejection sampling can be particularly useful in determining the normalizing factor for nonnormalized posterior distributions. The very general steps are as follows:

• find some other distribution that is convenient to sample from, where it must enclose the target distribution.

• compare generated values with the target distribution and decide which values to keep and which to reject.

• the integral of interest is the normalized ratio of these values.

The easier case is when the target function has bounded support: a defined, noninfinite support of some PDF or other function. Consider a function of interest $f(y)$ such that it can be evaluated for any value over the support of $y$: $S_y=[a,b]$. We want to put a “box” around the function and we have the two sides of the rectangle from $a$ and $b$. Since our interest is almost always centered on probability functions, the bottom of the box is the x-axis at zero. To determine the top of the box we can use various mathematical tools like analytical root finding, numerical techniques such as gridding, Newton–Raphson, EM, etc. to find $\tilde{y}$. So we have defined a rectangle that bounds all possible values for the pair $(y,f(y))$. Other geometric forms work too. The objective is now to sample the two-dimensional random variable in this rectangle by independently sampling uniformly over $[a,b]$ and $[0,\tilde{y}]$ then pairing to produce a point in 2-space. Finally all we have to do is count the proportion of points that fall inside the region of interest. The value of the integral, which is the area under the curve, is the ratio of the generated points under the curve compared to the total number of points scaled by the known size of the box:

$\frac{\text{count of points under the curve}}{\text{total number of generated points}}\times \text{dimension of the box}\underset{n\to \mathrm{\infty }}{⟶}{\int }_a^{b}f(y)dx.$

As with all Monte Carlo simulations, the error rate goes to zero proportional to $1/\sqrt{n}$.

A curve of interest is plotted in Figure 2 with 500 points generated uniform randomly inside the rectangle. Normally many more points are generated to ensure good Monte Carlo accuracy (low error), but the number is limited here for graphical purposes. There are 99 accepted points (marked with “+”) under the curve, and 401 rejected points (marked with “x”) above the curve. and 401 rejected red points above the curve. The rectangle has sides of length $[16.25,0.2197324]$ by construction. So the area under the curve (the integral) is calculated by:

$\begin{array}{rl}I(a,b)& =\text{proportion accepted}\times \text{size of box}\\ & =\frac{99}{500}\times (16.25\times 0.2197324)\\ & =0.706989\end{array}$

This was straightforward because we were able to easily put a box around the function of interest. The key feature of this example is that while we could not integrate the function, we were able to insert values of $y$ and get values of $f(y)$.

Figure 2 Example of Rejection Sampling

A more difficult version of this problem is when the target function, $f(y)$, has unbounded support, like say a gamma distribution with support $[0:\mathrm{\infty }]$. In this case we need to identify a majorizing function that fully envelopes the target in a way that a rectangle cannot, and can still have values drawn from. The rectangle was convenient because the size is simply height times width. Typically a chosen majorizing function, $g(y)$, is a PDF chosen such that $g(y)\ge f(y)$ for the full support of $f(y)$. To make sure that this property holds, it is always necessary to multiply the candidate generating function by a constant: $g^{\prime }(y)=kg(y)$. The steps are slightly different since we do not have the same height on the y-axis over the support. The steps are as follows:

• uniformly randomly generate many points on the support of $y$ (x-axis), $n$.

• for all $i=1,\dots n$, points calculate $a_i=f(y_i)/g^{\prime }(y_i)$.

• for all $i=1,\dots n$, generate a unique $\mathcal{U}(0,1)$ value, $u_i$.

• accept points if $u_i<a_i$.

• otherwise reject.

The setup is critical here. If the majorizing function is highly dissimilar than the target function, then a lot of values will be rejected and the algorithm will be inefficient. Rejection sampling is actually a fundamental part of one of the two important MCMC procedures.

Now we can state the Fundamental Principle of Monte Carlo Simulation. If we have many draws from a distribution of interest (or other functional forms), then we can use these draws to calculate any value of interest, with a known error rate. A key word in that definition is “any,” meaning that we can empirically calculate means, variances, intervals, quantiles, higher-order moments, and more. Thus, this process replaces analytical work for hard problems but does not deter us from getting values or summaries of interest.

Monte Carlo analysis has been around since the infancy of programmable computers in the early 1950s. Physicists such as Ulam, Von Neumann, Metropolis, and Teller were exploring subatomic properties developing nuclear weapons, but in some cases could not design the relevant physical experiments. So they turned to this new device called the computer and simulated the properties of interest. One important paper that is still cited and relevant to us in Section 5, is Metropolis et al. (Reference Metropolis, Rosenbluth, Rosenbluth, Teller and Teller1953). In this case the authors were interested in the paths of particles inside a box when the box is being squeezed down to a singularity.

General statistical interest in Monte Carlo methods in statistics started roughly in the 1980s when researchers were first having personal computers on their office desk. As these computers became more common and more powerful, Monte Carlo simulation became much more common because it made difficult calculations easier. If one were to look into the office of a research statistician in the 1970s, she would almost certainly be deriving properties analytically. Today, that same person is doing that same work on the same sorts of problems but mainly with various simulation tools. We now turn to the theoretical basics of Markov chains.

4.1 A Modeling Exercise

As a motivating example for our efforts to learn MCMC for Bayesian estimation, consider a dataset used in Gill and Casella (2009) looking at a puzzle in bureaucratic politics. These data contain every federal political appointee to full-time positions requiring Senate confirmation from November 1964 through December 1984 (collected by Mackenzie and Light, ICPSR Study Number 8458, Spring 1987). Their survey queries various aspects of the Senate confirmation process, acclamation to running an agency or program, and relationships with other functions of government. But these researchers needed to preserve anonymity of the interviewed so they embargoed some variables and randomly sampled 1,500 down to 512. This provides some methodological challenges, but we will ignore them here. So every respondent in the data ran a US federal agency or some major component of it. The puzzle is that even though all of them have to nominated by the president, underwent grilling on Capitol Hill, moved to an expensive real estate market, hired a lawyer or advisor for the process, often left a high-paying corporate position, and uprooted a family,  $\dots$ they stayed in the position on average for only about two years. Fortunately, there is a convenient outcome variable that may explain this phenomenon: stress, which can be considered as a surrogate measure for self-perceived effectiveness and job-satisfaction recorded as a five-point scale from “not stressful at all” to “very stressful.” The emphasis in this Element is on Bayesian regression models, so consider the explanatory variables in Table 1.

This is a rich dataset and there are many other variables to explore. Given that the chosen outcome variable is an ordinal scale, an “ordered choice model” makes sense. This name actually comes from economics, but applies to any regression outcome that is ordinal (i.e. not necessarily a “choice”). The model starts with assuming that there is a latent dimension underlying these observed ordinal categories. In a perfect world we could put some sort of a meter up to people’s heads and measure the exact level of stress in that person. In reality we have to live with some ordinal survey response, with five categories in this case. However, the model specification first assumes that we have the exact level measurement on the utility metric and then assumes that there are categories imposed on this latent scale from the survey wording. The setup looks like this:

${\mathbf{\text{U}}}_i:{\theta }_0\underset{c=1}{⟸⟹}{\theta }_1\underset{c=2}{⟸⟹}{\theta }_2\underset{c=3}{⟸⟹}{\theta }_3\dots {\theta }_{C-1}\underset{c=C}{⟸⟹}{\theta }_C$

Now the vector of (unseen) utilities across individuals in the sample, $\mathbf{\text{U}}$, is determined by a linear additive specification of explanatory variables: $\mathbf{\text{U}}=-{\mathbf{\text{X}}}^{\prime }\mathit{\beta }+\mathit{ϵ}$ as if we could see the latent dimension for real, where $\mathit{\beta }=[{\beta }_1,{\beta }_2,\dots ,{\beta }_k]$ does not depend on the ${\theta }_j$ cutpoints between categories, and the residuals have some distribution that we have not assigned yet, $\mathit{ϵ}\sim F_{\mathit{ϵ}}$, stated in cumulative distribution terms (CDF). These cutpoints between categories also need to be estimated on this latent scale and are sometimes called “fences” or “thresholds.” The key to understanding this model is to see that it initially assumes we know something that we cannot know: the true latent dimension. Another key point to note is that the model is initially given all in cumulative terms. This means that the probability that individual $i$ in the sample is observed to be in category $r$ or lower is:

$P({\mathbf{\text{Y}}}_i\le r|{\mathbf{\text{X}}}_i)=P({\mathbf{\text{U}}}_i\le {\theta }_r)=P(\mathit{ϵ}\le {\theta }_r+{{\mathbf{\text{X}}}_i}^{\prime }\mathit{\beta })=F_{{\mathit{ϵ}}_i}({\theta }_r+{{\mathbf{\text{X}}}_i}^{\prime }\mathit{\beta }),$(4.1)

which is differently signed than the popular MASS package in R, which states $logitP(Y<=k|x)={\zeta }_k-\eta$ from the help page. It is very irritating that the sign in front of ${\mathbf{\text{X}}}^{\prime }\mathit{\beta }$ varies by author and software, but we just have to live with that. One convenient way to figure out which sign is being used is to check the effect of an “obvious” explanatory variable. So, for example, in voting models we know that people get more conservative as they get older and richer. This is not a perfect recommendation but it helps to make up for technical omissions by some authors. Now specifying a logistic distributional assumption on the residuals produces a logistic cumulative specification for the whole sample:

$F_{\mathit{ϵ}}({\theta }_r+{\mathbf{\text{X}}}^{\prime }\mathit{\beta })=P(\mathbf{\text{Y}}\le r|\mathbf{\text{X}})=[1+\text{exp}(-{\theta }_r-{\mathbf{\text{X}}}^{\prime }\mathit{\beta }){]}^{-1}$(4.2)

However, a probit (cumulative normal) distribution is another popular choice.

Of course, this work is about Bayesian specifications for social science models so we need to specify prior distributions for the $\mathit{\beta }$ parameters. However, this differs from simpler regression models in that we also need to specify prior distributions for the $\theta$ cutpoint parameter as well. So, let’s specify the most simple possible and nearly conjugate forms as possible to keep this extended example as direct as possible:

$\begin{array}{r}p({\mathit{\beta }}_k)\sim \mathcal{N}({\mu }_{{\beta }_k},{\sigma }_{\beta }^2),k=1,\dots ,11\text{for the explanatory variables.}\\ p({\theta }_j)\sim \mathcal{N}(0,{\sigma }_{\theta }^2),j=1,\dots ,4\text{for the latent variable thresholds.}\end{array}$

Here is the problem; all this produces a posterior distribution according to

$\begin{array}{rl}\pi (\mathit{\beta },\mathit{\theta }|\mathbf{\text{X}},\mathbf{\text{Y}})& \propto L(\mathit{\beta },\mathit{\theta }|\mathbf{\text{X}},\mathbf{\text{Y}})p(\mathit{\theta })p(\mathit{\beta })p(b|\tau )p(\tau ).\\ & \propto \prod _{i=1}^n\prod _{j=1}^{C-1}\prod _{k=1}^p{\left[\mathrm{\Lambda }({\theta }_j-{{\mathbf{\text{X}}}_i}^{\prime }\mathit{\beta })-\mathrm{\Lambda }({\theta }_{j-1}-{{\mathbf{\text{X}}}_i}^{\prime }\mathit{\beta })\right]}^{z_{ij}}\\ & \times \text{exp}\left(-\frac{({\mathit{\beta }}_k-{\mu }_{{\beta }_k}{)}^2}{2{\sigma }_{\beta }^2}\right)-\frac{{\theta }_j^2}{2{\sigma }_{\theta }^2}-\frac{b_i^2}{2{\tau }^2}\end{array}$(4.3)

where $\mathrm{\Lambda }(\cdot )$ denotes the cumulative density function (CDF) of the logistic distribution, meaning this is actually a simplified version. What is this? It is a big joint posterior probability distribution with 16 dimensions, from the 11 explanatory variable parameters, the 4 threshold parameters, and the density (probability) over them. While this form tells us everything we need to know about the joint distribution of the model parameters, we need to marginalize (integrate) it for every one of these parameters to create an informative regression table. This may be very difficult or even impossible. We do not know. But wait! This is a standard model in the toolbox of most practicing quantitative social scientists and we just showed that specifying the most simple Bayesian version of this model leads to an intractable mathematical problem. This example is provided because it highlights the key problem of practicing Bayesian modelers prior to 1990. This was an insurmountable task for twentieth-century Bayesians that was solved in 1990. It turns out that there was a tool hiding in statistical physics and image restoration that completely solved this issue. Employing Markov chain Monte Carlo (Gibbs sampling in this case) allows us to report the standard regression table output in Table 2. Here we see a mixture of results quality. Two reliable effects are interesting: pre-hearing confirmation preparation lowers stress and greater work hours increase stress. Of course, these are nuanced qualitative effects that warrant further investigation.

Table 2Bureaucratic Politics Example Results

Table 2Long description

The table includes means, standard errors, and 95% highest posterior density (HPD) intervals for explanatory variables and threshold intercepts. HPD intervals are displayed as horizontal lines with a dashed vertical line at zero for reference. Several intervals lie entirely to one side of zero (e.g., *Career.Exec-Compet*, *Career.Exec-Liaison/Bur*, and *President Orientation*), indicating stronger evidence for nonzero effects. The estimated standard deviation of the latent scale is reported as 6.04 (SD = 1.20).

We know that scientific advances do not proceed in a linear fashion (Kuhn, Reference Kuhn1997), and this is a classic example of that. A review paper in 1990 (Gelfand & Smith, Reference Gelfand and Smith1990) changed all of this story. Now we turn our attention to this technology and the basic theory to understand it, starting with the mathematics of Markov chains.

4.2 Markov Chain Theory

The first step in this process is to understand some basic random variable definitions that occur over time. This is a case where the mathematical language does not reflect the simplicity of the principles, so it is important not to be intimidated by the vocabulary.

A stochastic process is a consecutive set of random quantities, $\theta$ s, defined on some known state space, $\mathrm{\Theta }$, indexed so that the order is known: $\{{\theta }^{[t]}:t\in T\}$. So $\theta$ is a generic random variable in the standard statistical sense except that it is observed in a serial manner rather than just all at once in a standard sampling scheme. Here $T$ just indexes the time of arrival and typically, but not necessarily, $T$ is the set of positive integers implying consecutive, even-spaced time intervals: $\{{\theta }^{[t=0]},{\theta }^{[t=1]},{\theta }^{[t=2]},\dots \}$. This state space is either discrete or continuous depending on how the variable of interest is measured. So, while the term “stochastic process” sounds very mathematical and intimidating it is really just the simple idea that random variables can have scheduled arrival times.

A Markov chain is a stochastic process with the property that any specified state in the series, ${\theta }^{[t]}$, is dependent only the previous value of the chain, ${\theta }^{[t-1]}$. Therefore, values are conditionally independent of all other previous values: ${\theta }^{[0]},{\theta }^{[1]},\dots ,{\theta }^{[t-2]}$. More formally:

$P({\theta }^{[t]}\in A|{\theta }^{[0]},{\theta }^{[1]},\dots ,{\theta }^{[t-2]},{\theta }^{[t-1]})=P({\theta }^{[t]}\in A|{\theta }^{[t-1]}),$(4.4)

where this is the probability that $\theta$ lands in some region $A$, which is any identified subset (an event or range of events) on the complete state space defined by $\mathrm{\Theta }$. We will use this $A$ notation extensively. Colloquially we can say:

Anyone who has ever owned a house cat or been around them much will certainly feel that cats have a very Markovian quality to their behavior. Again, this is a simple idea that unfortunately gets wrapped up in intimidating mathematical language. Necessarily Markov chains are based on a probabilistic process, a “kernel,” that makes a decision about the next draw of $\theta$ at time $T+1$ based only on the current existing draw of $\theta$ at time $T$. This leads to the standard analogy that the Markov chain is traveling through the state space defined by the support $\mathbf{\Theta }$. So we think of a serial path of the chain like an ant wandering around some topology on the ground. Consider Figure 3 where a 2-dimensional Markov chain is started at the point indicated by a larger dot (all Markov chains need an arbitrary starting point in the sample space), and the first step indicated by the bold black line segment segment. In this case we would call this a random walk Markov chain because we have not imposed any conditions on its behavior other than a random normal draw added to the current position.

Figure 3 Sample Markov Chain Path

So why is the Markovian property useful? This “short-term” memory property is very handy because when we impose some specific mathematical restrictions on the chain and allow it to explore a specific distributional form, it eventually finds the region of the state space with the highest density, and it will wander around there producing an empirical sample from it in proportion to the probability. If this distribution is the posterior from model, then we can use these “empirical” values as legitimate posterior sample values. Thus, difficult posterior calculations can be done with MCMC by letting the chain wander around “sufficiently long,” thus producing summary statistics from recorded values. For example, if we wanted marginal draws from the big awkward joint distribution in (4.3) then the Markov chain explores the multidimensional state space recording where it has been for us. This means that we can take the empirical values for each dimension (model parameter) and summarize them any way we want.

To clarify this strategy, consider the output of a Markov chain with 5 parameters run for 500,000 iterations. This means that it has visited 500,000 locations in the state space, and each 5-dimensional location is a row of a stored matrix. To get the marginal distribution analytically by contrast for just ${\theta }_1$, we would have to perform the integration in four steps:

$p({\theta }_1|\mathbf{\text{X}})={\int }_{{\theta }_2}{\int }_{{\theta }_3}{\int }_{{\theta }_4}{\int }_{{\theta }_5}p({\theta }_1,{\theta }_2,{\theta }_3,{\theta }_4,{\theta }_5|\mathbf{\text{X}})d{\theta }_{2}d{\theta }_{3}d{\theta }_{4}d{\theta }_5.$(4.5)

That’s a lot of work! And that’s for one of the dimensions; it needs to be done for all 5 of them. Recall from the structure of the joint distribution in (4.3) that this can be really hard or even impossible for humans with some model specifications. However, once we have run the Markov chain sufficiently long (more on this shortly) and stored the values in matrix form, then we get an empirical distribution for each of the parameters simply by looking down the columns. For instance, rather than do the integration process in (4.5), we would just pull out the values in column one to produce a vector of draws for the first parameter only. Doing the integration process gives an exact marginal distribution in mathematical terms, $p({\theta }_1|\mathbf{\text{X}})$. But wait! If we have half-a-million draws from this marginal distribution, then we actually have all of the information we need to know for scientific purposes. That is, we can take the mean and standard deviation of the vector to produce the information needed for a regression table. We can also calculate any other desired statistic or summary: quantiles, different measures of centrality and variance, probabilities over ranges, HPD intervals, and so on. Furthermore, given the Law of Large Numbers we know that we are very nearly precisely correct, given a very big number of draws like 500,000.

How does the Markov chain decide to move? The Markov chain moves based on a probability statement so we need to define how probability enters the decision. First, define the transition kernel, $K$, as a very general, and for the moment unspecific, mechanism for describing the probability of moving to some other specified state based on the current chain location. Here $K(\theta ,A)$ is a defined probability measure for all $\theta$ points in the state space to the set (subspace) $A\in \mathrm{\Theta }$. Therefore $K(\theta ,A)$ maps potential transition events to their probability of occurrence. Note that, in the current discussion, $\theta$ is a multidimensional point in the sample space and not a coefficient to be estimated. This is slightly confusing given the preceding discussion, but this notation is consistent with the literature.

Let us start with a simple discrete form of a two-dimensional state space, $\mathrm{\Theta }$, with 64 $\theta$ locations: the chess board. Discrete state spaces in this discussion are much easier to conceptualize so we will restrict ourselves to this form for now. It is also easy with this example to define various $A$ subspaces, like, say, the four squares in the top right-hand side of this chess board. The most flexible and powerful piece in chess is the Queen who can move any distance on the board vertically, horizontally, and diagonally. Now imagine a Super-Queen who can move from any square to any other square on the chessboard in one move, arbitrarily for from ${\theta }_i$ to ${\theta }_j$, where the transition kernel is uniform selection of any of the 64 squares. Thus we have 64 moves from some ${\theta }_i$ to another ${\theta }_j$, including moving to where we already are: ${\theta }_i\to {\theta }_i$. The latter scenario is included because the current location is always a valid destination for a Markov chain.

When the state space is discrete, $K$ is a matrix mapping for $k$ discrete elements in $A$, where each cell defines the probability of a state transition from the first term to all possible states. So for a rectangular 2-dimensional state space like the chess board ($k=64$), the transition matrix giving the probabilities of moving from any square to any other square for our Super-Queen is:

$K=\left[\begin{array}{lllll}p({\theta }_1,{\theta }_1)& p({\theta }_1,{\theta }_2)& \dots & p({\theta }_1,{\theta }_{k-1})& p({\theta }_1,{\theta }_k)\\ p({\theta }_2,{\theta }_1)& p({\theta }_2,{\theta }_2)& \dots & p({\theta }_2,{\theta }_{k-1})& p({\theta }_2,{\theta }_k)\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ p({\theta }_k,{\theta }_1)& p({\theta }_k,{\theta }_2)& \dots & p({\theta }_k,{\theta }_{k-1})& p({\theta }_k,{\theta }_k)\end{array}\right]$

where the row indicates where the chain is at this period. Thus the first row gives the 64 probabilities of starting at ${\theta }_1$ (say the top leftmost corner if we are mapping the chess board like reading a paragraph) to all of the 64 possible destinations starting with $p({\theta }_1,{\theta }_1)$ for returning to the current location using up a unit of time index, and proceeding to $p({\theta }_1,{\theta }_{64})$ for moving from the top left corner to the bottom right corner. It is important that each matrix row is a well-behaved probability function, $p({\theta }_i,{\theta }_j)\ge 0,\mathrm{\forall }i,j\in A$ in the sense of the Kolmogorov Axioms. Rows of $K$ sum to one and define a conditional probability mass function (PMF) since they are all specified for the same starting value as the condition and cover each possible destination in the state space: for row $i$: $\sum _{j=1}^{k}p({\theta }_i,{\theta }_j)i=1$. As a side note for now, when the state space is continuous, then $K$ is a conditional PDF: $K=f(\theta |{\theta }_i)$.

Consider a very simple mathematical example to get us started. Define a two-dimensional state space: a discrete vote choice between two political parties, a commercial purchase decision between two brands, a choice by a nation to go to war with another nation, and so on. In this setup, voters, consumers, predators, viruses, and so on. who normally select ${\theta }_1$ have a 70% chance of continuing to do so, and voters/consumers/predators/viruses/and so on who normally select ${\theta }_2$ have an 80% chance of continuing to do so. This provides the transition matrix $K$:

$\begin{array}{rl}& \stackrel{\text{next period}}{\overbrace{\begin{array}{rr}{\theta }_1& {\theta }_2\end{array}}}\\ \text{current period}& \left\{\begin{array}{l}\begin{array}{r}{\theta }_1\\ {\theta }_2\end{array}\left[\begin{array}{rr}0.7& 0.3\\ 0.2& 0.8\end{array}\right].\end{array}\right.\end{array}$

As a running hypothetical case, suppose we are doing market research for a major toothpaste brand and ${\theta }_1$ designates a purchase of Crest, and ${\theta }_2$ designates a purchase of Colgate. So from the preceding transition matrix, in any given time period, 70% of the Crest buyers stay loyal whereas 80% of the Colgate buyers remain loyal. This is just a mapping from one time period to the next since it defines a Markovian process where the past before that does not matter. All Markov chains require a starting point, in this context set by a human to

$S_0=\left[\begin{array}{rr}0.5& 0.5\end{array}\right],$

meaning that we initially recruit a sample that is one-half Crest buyers and one-half Colgate buyers. To go from this initial state to the first state, we simply postmultiply the starting state by the transition matrix:

$S_1=\left[\begin{array}{rr}0.5& 0.5\end{array}\right]\left[\begin{array}{rr}0.7& 0.3\\ 0.2& 0.8\end{array}\right]=\left[\begin{array}{rr}0.45& 0.55\end{array}\right]=S_1.$

So, observing our sample’s purchase at time period $T=1$, we see that 45% bought Crest and 55% bought Colgate. This series of observed purchases continues multiplicatively as long as we want to study our sample:

$\begin{array}{rl}\text{2nd state:}{S}_2=\left[\begin{array}{rr}0.45& 0.55\end{array}\right]\left[\begin{array}{rr}0.7& 0.3\\ 0.2& 0.8\end{array}\right]& =\left[\begin{array}{rr}0.425& 0.575\end{array}\right]\\ \text{3rd state:}{S}_3=\left[\begin{array}{rr}0.425& 0.575\end{array}\right]\left[\begin{array}{rr}0.7& 0.3\\ 0.2& 0.8\end{array}\right]& =\left[\begin{array}{rr}0.4125& 0.5875\end{array}\right]\\ \text{4th state:}{S}_4=\left[\begin{array}{rr}0.4125& 0.5875\end{array}\right]\left[\begin{array}{rr}0.7& 0.3\\ 0.2& 0.8\end{array}\right]& =\left[\begin{array}{rr}0.40625& 0.59375\end{array}\right].\end{array}$

If we continued this process a few more times, the resulting state of $[0.4:0.6]$ would emerge (try it!) and would never change no matter how many times we performed this operation. This is called the steady state, the stationary distribution, or the ergodic state (the latter for reasons that we will see later on page 31). Markov chains are not guaranteed to have this property and those that do not are called transient Markov chains.

In this section we describe another simple discrete Markov chain simplified from Bayer and Diaconis (Reference Bayer and Diaconis1992), and also described in Gill (Reference Gill2014). Suppose we have the following algorithm for shuffling a deck of cards: take the top card from the deck and with equal probability place it somewhere in the deck, then repeat many times. The objective (stationary distribution) is a uniformly random distribution in the deck: for any given order to the stack, the probability of any one card occupying any one position is 1/52. Is this a stochastic process? Sure, there is a set of discrete outcomes which are the ordering of the deck that change with each step and have a state space that is $52!$ big. Is this a Markov chain? Yes, clearly, because if we want to know the probability that the deck is in some specific ordering after the next step we only need to know the current status of the deck to compute the probabilities of the 52 possible arrangements at the next state. These are presumably all $\frac{1}{52}$ if we adhere to the condition of uniform random replacement presented earlier. So the previous states of the deck do not help us in any way. As long as we know the current state of the deck we have all the information required to use the transition kernel to move. What is the limiting distribution? For this, we will use simulation as a convenience as noted many times here.

For simplicity and without loss of generality we will use a deck with only $n=3$ cards. This is the same problem as using $n=52$ cards but is a lot easier to describe and visualize. With three numbered cards we have a sample space of size $3!$ given by $\{[1,2,3],[1,3,2],[2,1,3],[2,3,1],[3,1,2],[3,2,1]\}$. Given the uniformly random replacement construction of the algorithm, the transition kernel matrix is:

$K=\left[\begin{array}{cccccc}1/3& 0& 1/3& 1/3& 0& 0\\ 0& 1/3& 0& 0& 1/3& 1/3\\ 1/3& 1/3& 1/3& 0& 0& 0\\ 0& 0& 0& 1/3& 1/3& 1/3\\ 1/3& 1/3& 0& 0& 1/3& 0\\ 0& 0& 1/3& 1/3& 0& 1/3\end{array}\right]$

We will use the starting point $[1,2,3]$, which is analogous to a new deck of cards after one takes off the plastic wrapping and discarding the jokers. We start by defining this matrix and a function to implement the shuffling in both R and Python.

Box 7R Code for Card Shuffling

Code 1.07

Code 1.07

 ( K <- matrix(c(1/3,0,1/3,0,1/3,0,0,1/3,1/3,0,1/3,0,1/3, 
 0,1/3,0,0,1/3, 1/3,0,0,1/3,0,1/3,0,1/3,0,1/3,1/3,0, 
 0,1/3,0,1/3,0,1/3),nrow=6) ) 
 [,1] [,2] [,3] [,4] [,5] [,6] 
 [1,] 0.3333333 0.0000000 0.3333333 0.3333333 0.0000000 0.0000000 
 [2,] 0.0000000 0.3333333 0.0000000 0.0000000 0.3333333 0.3333333 
 [3,] 0.3333333 0.3333333 0.3333333 0.0000000 0.0000000 0.0000000 
 [4,] 0.0000000 0.0000000 0.0000000 0.3333333 0.3333333 0.3333333 
 [5,] 0.3333333 0.3333333 0.0000000 0.0000000 0.3333333 0.0000000 
 [6,] 0.0000000 0.0000000 0.3333333 0.3333333 0.0000000 0.3333333 
 shuffle <- function(kernel,iterations) { 
 deck <- c(1,0,0,0,0,0) 
 for (i in 1:iterations) { 
 print(deck) 
 deck <- deck %*%kernel 
 } 
 } 
 shuffle(K,20)

Box 8Python Code for Card Shuffling

Code 1.08

Code 1.08

 import numpy as np 
 K = np.array([ 
 [1/3, 0, 1/3, 1/3, 0, 0], 
 [0, 1/3, 0, 0, 1/3, 1/3], 
 [1/3, 1/3, 1/3, 0, 0, 0], 
 [0, 0, 0, 1/3, 1/3, 1/3], 
 [1/3, 1/3, 0, 0, 1/3, 0], 
 [0, 0, 1/3, 1/3, 0, 1/3] 
 ]) 
 def shuffle(kernel, iterations): 
 deck = np.array([1, 0, 0, 0, 0, 0]) 
 for i in range(iterations): 
 print(deck) 
 deck = np.dot(deck, kernel) 
 shuffle(K, 20)

Running this in either language provides the output probabilities at times in Table 3. Since each value is the probability of being in one of the six states at the time indicated by the iteration number it is easy to see that this Markov chain reaches its stationary distribution at about $T=15$ (subject to a little rounding of course). It also shows that the shuffling algorithm provides a fully shuffled deck since every one of the six outcomes has equal probability at this point.

Table 3Card Shuffling Probabilities

#	${\mathit{\theta }}_{\mathbf{\text{1}}}$	${\mathit{\theta }}_{\mathbf{\text{2}}}$	${\mathit{\theta }}_{\mathbf{\text{3}}}$	${\mathit{\theta }}_{\mathbf{\text{4}}}$	${\mathit{\theta }}_{\mathbf{\text{5}}}$	${\mathit{\theta }}_{\mathbf{\text{6}}}$
0	1.000000	0.000000	0.000000	0.000000	0.000000	0.000000
1	0.333333	0.000000	0.333333	0.333333	0.000000	0.000000
2	0.222222	0.111111	0.222222	0.222222	0.111111	0.111111
3	0.185185	0.148148	0.185185	0.185185	0.148148	0.148148
4	0.172839	0.160493	0.172839	0.172839	0.160493	0.160493
5	0.168724	0.164609	0.168724	0.168724	0.164609	0.164609
6	0.167352	0.165980	0.167352	0.167352	0.165980	0.165980
7	0.166895	0.166438	0.166895	0.166895	0.166438	0.166438
8	0.166742	0.166590	0.166742	0.166742	0.166590	0.166590
9	0.166692	0.166641	0.166692	0.166692	0.166641	0.166641
10	0.166675	0.166658	0.166675	0.166675	0.166658	0.166658
11	0.166669	0.166663	0.166669	0.166669	0.166663	0.166663
12	0.166667	0.166665	0.166667	0.166667	0.166665	0.166665
13	0.166667	0.166666	r0.166667	0.166667	0.166666	0.166666
14	0.166666	0.166666	0.166668	0.166666	0.166666	0.166666
15	0.166667	0.166666	0.166667	0.166667	0.166666	0.166666
16	0.166667	0.166667	0.166667	0.166667	0.166667	0.166667
17	0.166667	0.166667	0.166667	0.166667	0.166667	0.166667
18	0.166667	0.166667	0.166667	0.166667	0.166667	0.166667
19	0.166667	0.166667	0.166667	0.166667	0.166667	0.166667

This example also highlights a very big and important point about Markov chains. The numbers in Table 3 after convergence are marginal probabilities of being in any given state in an arbitrary distant time and for this ergodic Markov chain they will never change. But wait! The Markovian process is a process that produces values at the next iteration conditional on the value at the current iteration. “Conditional” is the operative word here and this is how all Markov chains work. So the contrast is between the conditional distribution that produces the next serially dependent value versus the unconditional future probabilities. Thinking about this another way, in the transition kernel the probability of being at ${\theta }_1$ and moving to ${\theta }_6$ is zero (we cannot go from $[1,2,3]$ to $[3,2,1]$ in one step), but in the stationary state for some large $T$ the probability of being in any state is $\frac{1}{6}$ because we do not stipulate the state at time $T-1$.

Now formalize the idea of the marginal distribution of the chain from applying the transition kernel at some distant time point $T=m$ giving marginal probabilities that do not change if the chain has reached its stationary distribution. Because this probability structure does not change anymore, we can simply label it as $p$. Since time point $m$ is far from the starting point at time point $0$, then we know that $p$ has been applied many times.

Now we need a set of conditions that assure us that a given Markov chain will eventually reach its stationary distribution, where the marginal distribution is constant and permanent. These mostly have fancy sounding mathematical names and our experience in teaching shows that they can be a little off-putting, but they are in fact fairly direct concepts. Each of these important terms is in bold for later reference.

A Markov chain is said to be homogeneous at step $t$ if the transition probabilities at this step do not depend on the value of $m$. More colloquially, elderly Markov chains behave exactly the same as young Markov chains. Now with non-homogeneous Markov chains the transition matrix is not constant but a function of time. This type are useful for modeling probabilistic decay processes, but not for our purposes here.

Recall that Markov chains operate on states, which are subspaces of the state space, $\mathrm{\Theta }$, as in our example of subspaces on of the chess board. These can be individual discrete units or collections of them. The specific state is absorbing if, once the chain enters this state, it cannot leave: $p(A,A^c)=0$. The state is transient if the probability of the chain not returning to this state is nonzero: $1-p(A,A)>0$, meaning that there is a finite number of visits in infinite time for the chain. State $A$ is closed to state $B$ if a Markov chain on $A$ cannot reach $B$: $p(A,B)=0$. State $A$ is closed in the general sense if it is absorbing.

The state $A$ is irreducible if for all substates ${\theta }_i$ and ${\theta }_j$ in this state, the two substates “communicate,” meaning that every reached point or collection of points (necessarily a subspace in the continuous case) can be reached from every other reached point or collection of points:

$p({\theta }_i,{\theta }_j)\ne 0,\mathrm{\forall }{\theta }_i,{\theta }_j\in A.$

So being in any one part of $A$ does not prevent the Markov chain from reaching any other part of $A$.

The next definition is a “collector” term that combines some of these properties. If a state is (1) closed, (2) discrete, or continuous but finitely bounded, and (3) irreducible, then this state and all subspaces within this subspace are called recurrent. Markov chains operating on recurrent state spaces are labeled as recurrent themselves. A formal definition is:

Less formally: (1) a recurrent Markov chain run for an infinite amount of time moves into a recurrent state, (2) stays there forever, and (3) visits every subspace an infinite number of times in the limit. This sounds like an odd treatment of infinity, but mathematicians remind us that there are different “flavors” of infinity so it works. Also, a Markov chain is positive recurrent if the mean time to return to $A$ is bounded, Otherwise it is called null recurrent. Additionally, recurrence is an absolute condition in that irreducible chains are either recurrent or transient. For continuous state spaces that are unbounded, like, say, the support of the normal distribution, there is a stricter variant of recurrence called Harris recurrence where the probability of an infinite number of visits to every possible substate of every possible $A$ is one.

Recurrence is an abstract mathematical concept, but it helps to visualize the chess board again. Suppose $A$ is the topmost four squares on the right-hand side. If we run the Super-Queen for an infinite number of moves, then it will visit $A$ an infinite number of times and 1 out of 16 times as well. Now consider a different version of $A$: the strip of 8 vertical squares on the far left-hand side. The abstract notion of infinity is odd indeed. If the Super-Queen is run under recurrence for an infinite number of moves, it will visit $A$ $1/8$ th of the time, which is also infinity. In addition, the union of a set of recurrent states (nonempty, and bounded or countable) or Harris recurrent states is a new state that is closed and also irreducible. So the infinity property and the recurrence associated with it also works for the first two columns on the left-hand side consisting of 16 squares.

Define $\pi (\theta )$ as the stationary distribution of the Markov chain for $\theta$ on the state space $A$. Recall that $p({\theta }_i,{\theta }_j)$ is the probability that the chain will move from ${\theta }_i$ to ${\theta }_j$ at some arbitrary step $t$, and ${\pi }^t(\theta )$ is the corresponding marginal distribution. The stationary distribution satisfies

$\begin{array}{rl}\sum _{{\theta }_i}{\pi }^t({\theta }_i)p({\theta }_i,{\theta }_j)={\pi }^{t+1}({\theta }_j)& \text{Discrete case, and}\\ \int {\pi }^t({\theta }_i)p({\theta }_i,{\theta }_j)d{\theta }_i={\pi }^{t+1}({\theta }_j)& \text{Continuous case}.\end{array}$

meaning that application of the kernel to the current distribution probabilities returns the same probabilities as seen in the two simple numerical examples presented earlier. This means that the marginal distribution remains fixed when the chain reaches the stationary distribution and we might as well drop the superscript designation for iteration number and just use $\pi (\theta )$. In shorthand: $\pi =\pi p$. Once the chain reaches its stationary (absorbing) distribution, it stays in this distribution and moves about, or “mixes,” throughout the subspace according to the marginal distribution $\pi (\theta )$ forever. Notice that recurrence gives the range restriction property whereas stationarity gives constancy of the probability structure that dictates movement.

The period of a Markov chain is the length of time required to repeat an identical cycle of chain values. Actually, periodicity is bad and we want an aperiodic Markov chain: the repeating length is the trivial value of one. The problem with periodicity is that it destroys the probabilistic property of the transition kernel. Caveat: all Markov chains run on computers are actually periodic. John Von Neumann’s famous quote on this is: “Anyone who considers arithmetical methods of producing random digits is, of course, in a state of sin” (Von Neumann, Reference Von Neumann1963). Fortunately with modern pseudo-random number generators on computers the period is so large that we need not be concerned as previously noted.

Now we are ready for the big result that all of these conditions have been building up to. If a chain is positive Harris recurrent (unbounded continuous case) or recurrent (discrete/bounded continuous case), and aperiodic, then we call it ergodic. Ergodic Markov chains have the property:

$\underset{n\to \mathrm{\infty }}{\text{lim}}{K}^n({\theta }_i,{\theta }_j)=\pi ({\theta }_j),$

for all ${\theta }_i$, and ${\theta }_j$ in the subspace. What this means is that once a specified chain has reached its ergodic state, sample values behave as if they were produced independently by the posterior of interest from the model. The Ergodic Theorem is the equivalent of the Strong Law of Large Numbers but for Markov chains, since it states that any specified function of the posterior distribution can be estimated with samples from a Markov chain in its ergodic state because averages of sample values give strongly consistent parameter estimates. Suppose ${\theta }_{i+1},\dots ,{\theta }_{i+n}$ are $n$ (not necessarily consecutive) values from a Markov chain that has reached its ergodic distribution before time $i$, a statistic of interest, $h(\theta )$, can be calculated empirically:

$\hat{h}({\theta }_i)=\frac{1}{n}\sum _{j=i+1}^{i+n}h({\theta }_j)=\int h(\theta )d\pi (\theta ),$(4.6)

where the summation starts at $i+1$ since we assume that the Markov chain is in its stationary distribution by time $i$, and we consider $n$ values thereafter. This is the Ergodic Theorem. If $h$ is the identity function, then this is just the calculation of the expected value of the random variable $\theta$. Notice that once we take this sum we have discarded the order in which the samples arrived. For finite quantities this sum converges almost surely:

$p[\hat{h}({\theta }_i)\to h(\theta ),\text{as}n\to \mathrm{\infty }]=1$(4.7)

Even though Markov chain values, by their very definition, have serial dependence, the mean of the chain values provides a strongly consistent estimate of the true parameter after convergence. We get an important theoretical result from the Central Limit Theorem. For a given empirical estimator $\hat{h}({\theta }_i)$ with bounded limiting variance, we get:

$\sqrt{n}\frac{\hat{h}({\theta }_i)-h(\theta )}{\sqrt{\text{Var}(\hat{h}({\theta }_i))}}\underset{n\to \mathrm{\infty }}{⟶}\mathcal{N}(0,1)$

since $\hat{h}({\theta }_i)=\frac{1}{n}\sum _{j=i+1}^{i+n}h({\theta }_j)$ is really a mean. The Central Limit Theorem really simply means that nature loves the normal distribution. Here because we are looking at averages taken with large numbers of draws it shows up naturally (pun intended), which is very cool.

What this section provides is the necessary background to understand Markov chains. We will next move on to Markov chain Monte Carlo (MCMC) as a Bayesian estimation procedure. Unlike some standard estimation procedures, either in closed form or through maximization, it is important to understand how Markov chains work since it is easier to apply MCMC incorrectly, which produces misguided results. The mathematics of this section may seem abstract at first, but they really describe a reasonably simple process whereby a stochastic process moves around a state space under a set of very particular rules.

5.1 The Gibbs Sampler

In Section 4 we saw that the mechanism that dictates how a Markov chain moves is a transition kernel ($K$). For discrete state spaces, these were shown to be probability matrices and for continuous state spaces, we saw that these were merely probability density functions (PDFs). In general practice, the transition kernel is an algorithm: either the Gibbs sampler or the Metropolis–Hastings rejection procedure, or possibly a hybrid of the two. We now introduce the Gibbs sampler and save the short discussion of Metropolis–Hastings for later.

Gibbs sampling was introduced in a paper by Geman and Geman (Reference Geman and Geman1984) that appeared in the IEEE Transactions on Pattern Analysis and Machine Intelligence, and subsequently went ignored for years (at least by statisticians) until it was discovered by Gelfand and Smith (Reference Gelfand and Smith1990). This is a great example of Stigler’s Law of Eponymy that no discovery in science is ever named for its actual inventor since “Gibbs” comes from the distribution used in the Geman and Geman (Reference Geman and Geman1984) paper not the authors. The Gelfand and Smith publication fundamentally changed statistics and is easily one of the most important contributions to statistics even though it was essentially a review paper. To reflect on the importance, one merely needs to try and read the Geman and Geman paper. So suddenly a huge class of Bayesian models that were difficult or impossible to estimate suddenly became accessible. For a fascinating history of this era see Robert and Casella (Reference Robert and Casella2011).

The Gibbs sampler is a transition kernel created by a series of full conditional distributions that are sampled in cycles. So it is a Markov chain updating scheme based on conditional probability statements where every parameter has its own specification. Suppose that the limiting distribution of interest is $\pi (\mathit{\theta })$ where $\mathit{\theta }$ is a $k$-length vector of coefficients to estimate. The objective is to produce a Markov chain that cycles through these conditional statements moving toward and then around this distribution. For the moment, consider a joint posterior that we wish to marginalize with $k=3$ that is given by $\pi ({\theta }_1,{\theta }_2,{\theta }_3|\mathbf{\text{X}},\mathbf{\text{y}})$ with data $\mathbf{\text{X}}$ and $\mathbf{\text{y}}$ in the regression context. The full algorithm at time $t$ is given by drawing three random values according to the following expressions:

$\begin{array}{rl}{\theta }_1^{[t]}& \sim p\left({\theta }_1|{\theta }_2^{[t-1]},{\theta }_3^{[t-1]},\mathbf{\text{X}},\mathbf{\text{y}}\right),\\ {\theta }_2^{[t]}& \sim p\left({\theta }_2|{\theta }_1^{[t]},{\theta }_3^{[t-1]},\mathbf{\text{X}},\mathbf{\text{y}}\right),\text{and}\\ {\theta }_3^{[t]}& \sim p\left({\theta }_3|{\theta }_1^{[t]},{\theta }_2^{[t]},\mathbf{\text{X}},\mathbf{\text{y}}\right).\end{array}$

This is in total a single step of the Gibbs sampler. Clearly this is a Markov chain since the decision to move at time $t$ is only conditional on information at time $t-1$ in total. Importantly, notice that each draw is done with the most recent information from the other parameters: ${\theta }_1$ uses information strictly at time $t-1$ since it is drawn first, ${\theta }_2$ uses information at time $t$ from ${\theta }_1$ but information at time $t-1$ for ${\theta }_3$, and ${\theta }_3$ uses information completely from the current step since it is sampled last. Actually, the order of this within-step sampling does not matter, meaning 1,2,3 are in arbitrary order, but it is critical that no matter what the order is, each sub-sampling step is done conditional on the most “modern” version of the other parameters.

Consider data that are normally distributed: $y_1,\dots ,y_n\sim \mathcal{N}(\mu ,{\sigma }^2)$, meaning that we have two parameters to estimate from the joint posterior:

$\pi (\mu ,{\sigma }^2|\mathbf{\text{y}})\propto L(\mu ,{\sigma }^2|\mathbf{\text{y}})p(\mu )p({\sigma }^2),$(5.1)

with researcher-specified normal and gamma priors:

$p(\mu )\sim \mathcal{N}(0,s^2)p({\sigma }^2)\sim \mathcal{G}(\alpha ,\beta ).$(5.2)

The full conditional distribution for $\mu$ is given by:

$p(\mu |{\sigma }^2,\mathbf{\text{y}})\propto \mathcal{N}\left(\frac{{\sigma }^{-2}}{n{\sigma }^{-2}+s^{-2}}\sum _{i=1}^n{y}_i,\frac{1}{n{\sigma }^{-2}+s^{-2}}\right).$(5.3)

This notation here uses a variance, ${\sigma }^2$, although in Bayesian specifications it is common to use the precision: $\tau ={\sigma }^{-2}$, meaning high variance is low precision. The full conditional distribution for ${\sigma }^2$ is given by the following formula:

$p({\sigma }^2|\mu ,\mathbf{\text{y}})\propto \mathcal{G}\left(\alpha +\frac{n}{2},\beta +\frac{1}{2}\sum _{i=1}^n(y_i-\mu {)}^2\right).$(5.4)

Notice that in both cases the Bayesian principle that posterior information is a compromise between prior information ($\alpha ,\beta ,{\sigma }^2$) and data information ($n,\mathbf{\text{y}}$) is apparent, and in both cases the data dominate in the limit. So now we have the complete recipe for the Gibbs sampler where we alternate drawing from the full conditional distributions for $\mu$ and ${\sigma }^2$.

As an example, Table 4 contains the percent of the total US vote on “others” (non-Democratic and non-Republican) in the 2020 presidential election by state and the District of Columbia. Some might call these “independents” but there is a literature in political science that finds more nuances. These data are unimodal and almost symmetric and are thus a good candidate for modeling with a normal distribution. Running the Gibbs sampler 20,000 iterations and dispensing with the first 10,000 draws to let the Markov chain reach its stationary distribution (we will discuss “burn-in” in Section 8). This is specified in the following code which returns an empirical mean of 0.347 for $\mu$ and 0.008 for ${\sigma }^2$. Interestingly, this model specification for 51 cases is fairly sensitive to the priors selected. This is often not the case with decent sized social science data, so we encourage readers to download the code and data for experiential purposes on this type of situation.

Box 9R Code for 2020 Election Percents

Code 1.09

Code 1.09

 elect <- read.table(election.2020.dat",header=TRUE) 
 norm.gibbs <- function(y,s2.inv,alpha,beta,n.sims) { 
 mcmc.out <- matrix(NA,nrow=n.sims,ncol=2) 
 mcmc.out[1,] <- c(0,1) 
 n <- length(y) 
 for (i in 2:nrow(mcmc.out)) { 
 sig2.inv <- mcmc.out[(i-1),2] 
 mcmc.out[i,1] <- rnorm(1, 
 (sig2.inv/(n*sig2.inv + s2.inv))*sum(y), 
 1/(n*sig2.inv + s2.inv)) 
 mcmc.out[i,2] <- rgamma(1, shape=alpha+n/2, 
 rate=beta + 0.5*sum(y-mcmc.out[i,1])^2) 
 } 
 return(mcmc.out) 
 } 
 gibbs.vals <- norm.gibbs(elect$Vote.pct,2,2,12,20000) 
 apply(gibbs.vals[10001:20000,],2,mean)

Box 10Python Code for 2020 Election Percents

Code 1.010

Code 1.010

 import pandas as pd 
 import numpy as np 
 from scipy.stats import norm, gamma 
 elect = pd.read_csv("election.2020.csv", header=0) 
 def norm_gibbs(y, s2_inv, alpha, beta, n_sims): 
 mcmc_out = np.zeros((n_sims, 2)) 
 mcmc_out[0, :] = [0, 1] 
 n = len(y) 
 for i in range(1, n_sims): 
 sig2_inv = mcmc_out[i-1, 1] 
 mu_mean = (sig2_inv / (n * sig2_inv + s2_inv)) * 
 np.sum(y) 
 mu_var = 1 / (n * sig2_inv + s2_inv) 
 mcmc_out[i, 0] = norm.rvs(loc=mu_mean, 
 scale=np.sqrt(mu_var)) 
 shape = alpha + n / 2 
 rate = beta + 0.5 * np.sum((y - 
 mcmc_out[i, 0])**2) 
 mcmc_out[i, 1] = gamma.rvs(a=shape, scale=1/rate) 
 return mcmc_out 
 gibbs_vals = 
 norm_gibbs(elect[’Vote.pct’], 2, 2, 12, 20000) 
 mean_vals = np.mean(gibbs_vals, axis=0) 
 print(mean_vals)

Table 42020 Independent Vote Percent

Table 4Long description

The table lists independent vote percentages in the 2020 election for all U.S. states and Districtic of Columbia. The top five are Alaska 4.39%, Utah 4.22%, Wyoming 3.51%, Vermont 3.24%, and Washington 3.26%. The bottom five are Florida 0.92%, Pennsylvania 1.15%, Georgia 1.25%, New Jersey 1.27%, and Mississippi 1.34%.

5.2 Model Fitting with MCMC Software

The MCMC software implementation JAGS is based on the language BUGS that was developed at Cambridge Biostatistics Unit and released in 1991 (Gilks, Thomas, & Spiegelhalter, Reference Gilks, Thomas and Spiegelhalter1994). It is a flexible and powerful platform developed by Martyn Plummer and is also distributed freely. The syntax is deliberately relatively R-like, but unlike R there are much fewer functions to work with because it is not intended to be an all-purpose statistical environment. JAGS can be run stand-alone but it is much easier to call it from R using the package rjags or in Python using the package PyJAGS.

There are five general steps to running the model and producing results:

1. 1. Specify the distributional features of the model, and the quantities to be estimated in a text editor and import into R as a function.

2. 2. Specify vectors for starting points and names of coefficients.

3. 3. Run the model with the rjags function.

4. 4. Assess convergence using the diagnostics in the R package superdiag.

5. 5. Summarize empirically by variable as if the samples were data.

Most of the work is in the first step, but a nice feature of the software, especially the ones we focus on in this book, is that the model is specified very much as one would write it out on paper as shown in the following example. Interestingly, unlike other programming languages, statements are not processed serially, they constitute a full specification. Still, there is a structure to writing the code.

The default MCMC transition kernel in JAGS is the Gibbs sampler and therefore the first step in the software identifies the full conditional distributions for each variable in the model. This is done by the software not by the user. If the software cannot determine the full conditional distribution for a given parameter, then it automatically switches to the Metropolis–Hastings algorithm for that dimension only, making the result a hybrid sampler. These statements show the beauty of how this software is constructed: the user codes a statistical model and JAGS encodes the full conditional distributions for Gibbs sampling.

Since the JAGS syntax is the same for R and Python, it should be fairly easy to adapt the R JAGS code to Python and vice versa. PyJAGS provides a Python interface to JAGS, which uses the exact same JAGS syntax for model specification and Python code to specify data and run sampler. Therefore, instead of repeating the JAGS code for Python, we want to introduce PyMC for Python implementation of the Bayesian regression models, which is a popular probabilistic programming library for building Bayesian models in Python. Developed as an open-source project, PyMC offers a powerful and flexible platform for Bayesian statistical modeling. Also PyMC is built on top of Theano (now maintained as Aesara), a Python library that is for working with multi-dimensional arrays including computing gradients using automatic differentiation. This makes it particularly adept at handling more sophisticated models that involve large datasets and complex posterior distributions. The syntax of PyMC is intentionally Pythonic, making it accessible to users familiar with Python and conducive to integration into larger Python applications. Additionally, PyMC shares similarities with JAGS in its expressive model specification, making it relatively easy for R JAGS users to transition to PyMC. It should be noted that PyMC uses advanced MCMC algorithms such as the No-U-Turn Sampler (NUTS), a variant of Hamiltonian Monte Carlo (HMC), which is more efficient for high-dimensional parameter spaces.

5.2.1 Murder Sentencing in Georgia Data

In a seminal study, David C. Baldus, Charles Pulaski, and George Woodworth (e.g. the Baldus study) looked at the potential disparity in the imposition of the death sentence in Georgia based on the race of the murder victim and the race of the defendant (Baldus, Pulaski, & Woodworth, Reference Baldus, Pulaski and Woodworth1983). This is actually two studies, the second one examining about 762 cases with a murder conviction in Georgia from March 1973 to December 1979. The data contains 160 covariates, including legal background, crime description, and demographics. From the 1970 US Census 1,187,149 out of 4,589,575, or about 26%, of Georgia residents were Black. The death penalty was imposed during this time according to:

22% cases of Black defendant, White victim

8% cases of White defendant and White victim

1% of cases of Black defendant and Black victim

3% of cases of White defendant and Black victim.

The Baldus study was cited in the US Supreme Court case McClesky v. Kemp (1987), in which a Black defendant (McClesky) was sentenced to death for killing a White police officer in Georgia. The central argument was that the sentence violated the Equal Protection clause of the Fourteenth Amendment, since statistically he was more likely to get the death penalty since the victim was White. The Supreme Court rejected McClesky’s argument (5-4) on the grounds that statistical trends did not effectively prove the existence of discrimination among the jury who decided his particular case. Supreme Court McClesky was executed in 1991.

The outcome variable is sentence, which is dichotomous according to: 0= life sentence, and 1=death penalty given. Consider the general statement of the regression structure in (1.1) where we will now specify a link function that is appropriate for 0/1 outcomes. This is a preview of the regression models in Sections 6 and 7. Of the 160 possible variables we choose the 9 in Table 5 as explanatory variables for this example. We matched 287 out of 454 White cases to the existing 287 Black cases using Mahalanobis distance matching with Match in R (Sekhon, Reference Sekhon2008). Since the analysis is done in R (by calling JAGS using the rjags package) or in Python, we will store these data in a list that looks like the following code boxes. In these boxes we only show a representative fraction of the data for space purposes. Notice also that we include data size definitions first in the list as good programming practice. For R the equal sign is required here not the “<-” notation.

Box 11R Data List for the Baldus Example

Code 1.011

Code 1.011

 baldus.list <- list(NUM.CASES = 358, NUM.VARS = 10, 
 sentence = c( 0, 1, 1, 0, 0, 0, 0, 0, 0, 0,...), 
 race = c( 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...), 
 educatn = c( 2, 2, 3, 1, 2, 1, 1, 3, 2, 2,...), 
 SES = c( 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...),...

Box 12Python Data Format for the Baldus Example

Code 1.012

Code 1.012

 NUM_CASES = len(baldus_data) 
 NUM_VARS = 10 
 sentence = baldus_data[’sentence’].values 
 race = baldus_data[’race’].values 
 educatn = baldus_data[’educatn’].values 
 SES = baldus_data[’SES’].values 
 num_chld = baldus_data[ʼnum.chld’].values 
 defense = baldus_data[’defense’].values 
 vict_sex = baldus_data[’vict.sex’].values 
 pr_incrc = baldus_data[’pr.incrc’].values 
 aggrevat = baldus_data[’aggrevat’].values 
 mitcir = baldus_data[’mitcir’].values

Generally, it is most convenient to write actual code in an ASCII text editor to create a file that the rjags or PyJAGS package functions can find. We will now break down the components of the logit regression model that models the dichotomous sentencing outcome. The outside structure of the model is simply a function definition in R with no terms passed-in:

Box 13R Outer Shell the Baldus Example

Code 1.013

Code 1.013

 baldus.model <- function() { 
 : 
 : 
 }

Table 5Explanatory Variables for the Death Penalty Model

Table 5Long description

Each row defines a predictor and its coding: race (0 = White, 1 = Black), education level (1 - 3, 1 = middle school or lower, 3 = high school degree), socioeconomic status (0 = not low wage, 1 = low wage), number of children (1–9+), defense type (1 = retained, 2 = appointed), victim’s sex (1 = male, 2 = female), prior incarceration record, presence of aggravated methods, and presence of mitigating circumstances.

In PyMC, we use a model context to define model components, including the priors, likelihood, and the sampler. Model context is initiated using a with statement. This context manager encapsulates the entire model, ensuring that all the variables and definitions are scoped correctly within this block. The as model part will store this model context as an object called model in the environment, which can be called later.

Box 14Python Model Context for the Baldus Example

Code 1.014

Code 1.014

 with pm.Model() as model:

Now it is our responsibility to tell JAGS to loop through the data and perform a statistical procedure. This looks very much like a model written out by hand. First, we sum up the linear statistical component $\mathbf{\text{X}}\mathit{\beta }$ by multiplying each variable per case by the corresponding coefficient to be estimated. The requirement in rjags is to put the link function on the left-hand side of the equation rather than specify the inverse link function on the right-hand side of the equation as most text books do: logit(p[i]). Here p[i] is the outcome variable for the $i$ th case. Notice that this is the exact structure identified in (1.1), where the link function is now the logit function: $g(p)=\text{log}\left[p/(1-p)\right]$. Finally inside the cases loop we connect the outcome variable to the distributional assumption, which in this case is the Bernoulli PMF. This part of the code including the looping is shown in the following R code box. Just a reminder here that we will cover linear and nonlinear regression in much greater detail in the subsequent two sections.

Box 15R Statistical Specification for the Baldus Example

Code 1.015

Code 1.015

 for (i in 1:NUM.CASES) { 
 logit(p[i]) <- beta[1] + beta[2]*race[i] 
 + beta[3]*educatn[i] + beta[4]*SES[i] 
 + beta[5]*num.chld[i] + beta[6]*defense[i] 
 + beta[7]*vict.sex[i] + beta[8]*pr.incrc[i] 
 + beta[9]*aggrevat[i] + beta[10]*mitcir[i] 
 sentence[i]   dbern(p[i]) 
 }

In PyMC, the model context will handle looping through the data and perform the statistical procedure. The model specification also closely mirrors the conceptual framework of the model. Unlike JAGS, where the link function is explicitly specified in the model structure, PyMC handles this through the sigmoid function on the right-hand side of the equation.

Box 16Python  Statistical Specification for the Baldus Example

Code 1.016

Code 1.016

 logit_p = (beta[0] + beta[1]*race + 
 beta[2]*educatn + beta[3]*SES + 
 beta[4]*num_chld + beta[5]*defense + 
 beta[6]*vict_sex + beta[7]*pr_incrc + 
 beta[8]*aggrevat + beta[9]*mitcir) 
 p = pm.math.sigmoid(logit_p) 
 sentence_obs = pm.Bernoulli(’sentence_obs’, 
 p=p, 
 observed=sentence)

Finally, since this is a Bayesian model, we need to specify the prior distributions. This is done in a very R-like syntax such as g.coef   dgamma(alpha, beta), for each variable to be estimated. Here we could have done this 7 times for the seven beta coefficients, but it is more elegant to do it in a loop if they all get the same priors. This is why the betas are indexed. So now we add the final statement, which is outside the NUM.CASES loop in the NUM.VARS loop to finish the code.

Box 17Full JAGS Code for Death Penalty Model

Code 1.017

Code 1.017

 baldus.model <- function() { 
 for (i in 1:NUM.CASES) { 
 logit(p[i]) <- beta[1] + beta[2]*race[i] 
 + beta[3]*educatn[i] + beta[4]*SES[i] 
 + beta[5]*num.chld[i] + beta[6]*defense[i] 
 + beta[7]*vict.sex[i] + beta[8]*pr.incrc[i] 
 + beta[9]*aggrevat[i] + beta[10]*mitcir[i] 
 sentence[i]   dbern(p[i]) 
 } 
 for (j in 1:NUM.VARS) beta[j]   dnorm(0.0,0.01) 
 }

Just as with JAGS, we start by defining the model, specifying priors, likelihoods, and, if necessary, transformations directly using the PyMC library. It also provides distribution functions to specify the prior distributions. In our case, we use the Normal() function to draw normal distributions with specified mu (mean) and sigma (standard deviation) for the priors of beta. The Normal() function also allows you to draw multiple values at a time by using the shape argument. One major difference between is that in PyMC the model specification and sampler need to be put in the same model context. PyMC code is also highly expressive just like rjags.

Box 18Full PyMC Code for Death Penalty Model

Code 1.018

Code 1.018

 with pm.Model() as model: 
 beta = pm.Normal(’beta’, mu=0, 
 sigma=10, shape=NUM_VARS) 
 logit_p = (beta[0] + beta[1]*race + 
 beta[2]*educatn + beta[3]*SES + 
 beta[4]*num_chld + beta[5]*defense + 
 beta[6]*vict_sex + beta[7]*pr_incrc + 
 beta[8]*aggrevat + beta[9]*mitcir) 
 p = pm.math.sigmoid(logit_p) 
 sentence_obs = pm.Bernoulli(’sentence_obs’, 
 p=p, 
 observed=sentence) 
 baldus_out = pm.sample(draws=50000, 
 tune=20000)

All JAGS models have essentially this structure above: loop through the data statements to provide the statistical model, separately specify the prior distributions on all of the unknown parameters in the model, and add any additional features as appropriate like longitudinal and hierarchical components. Now we need to identify a vector of starting points, and it is convenient to create a list of the coefficient names so as not to clutter up the model call. In this case since we indexed beta only a single character object is required. Calling rjags means specifying where the model, the data, the parameters, and the initial values are located. It is also required to specify the number of iterations. The full production of marginal posterior values involves two step shown below: create the model object, and then run this model object to draw samples. Note that in PyMC we don’t need a separate sampling procedure. The sampler is already defined in the model context and will run once we execute the code above.

Box 19R Code for Calling the Death Penalty Model

Code 1.019

Code 1.019

 baldus.inits <- list("beta" = rep(1,10)) 
 baldus.params <- c("beta") 
 baldus.jags.obj <- jags.model( 
 file="./BOOKS/Book.Bayes.IV/baldus.model.jags", 
 data=baldus.list, inits=baldus.inits, n.chains = 1) 
 baldus.out1 <- coda.samples(baldus.jags.obj, 
 variable.names=baldus.params, n.iter=50000)

This produces the following output when we summarize the columns as data. These are deliberately portrayed as if from a standard logit regression although Bayesians usually prefer more information.

	Coef	StdErr
(Intercept)	-3.44669686	1.11392418
race	0.87493310	0.28946874
educatn	-0.01799630	0.19581014
SES	-1.69381203	0.29300638
num.chld	0.01006399	0.08053441
defense	1.18060470	0.33752817
vict.sex	0.91479574	0.34767625
pr.incrc	0.22999805	0.32048919
aggrevat	1.00660730	0.42226770
mitcir	-1.09579542	0.28368242

We see here five of the eight coefficients imply reliable effects with these point estimations: posterior means and standard errors. Consider the first two posterior coefficient distributions shown in Figure 5. For the race variable on the left (White=0, Black=1) we see that 0.0012 of the density (pictured in red) is below zero indicating that there is a 0.9982 probability that this effect is positive meaning that Black defendants are reliably more likely to get the death penalty than White defendants, given the model. For the education variable (middle school=1, some high school = 2, high school =3) the story is quite different: 0.5413 of the density is below zero meaning that we can make no reliable conclusion about the effect of this variable as there is a near equal probability of its effect being positive and negative. So posterior illustrations can be much more informative than standard regression tables.

Figure 5 Posterior Density for ${\beta }_1$ and ${\beta }_2$