A Bayesian Primer for the Organizational Sciences

An Introduction to Bayes’s Theorem

Thomas Bayes’s (1702-1761) original theorem is a mathematical equation built on simpler probability axioms (Jackman, 1999). Its basic form is

p ( A | B ) = p ( B | A ) p ( A ) p ( B ) ,

1

where p indicates probability, and the “|” symbol denotes a conditional probability, the probability of one event given that another event has occurred—for example, the probability that I will eat given that I am hungry. (Thus, the left side of the equation, p ( A | B ) , is the probability of A given that B has occurred.) Although this equation may seem complex, its utility is straightforward: It allows us to calculate any conditional probability of interest, provided that we can specify the terms on the right side of the equation. As previously stated, in statistics we are typically interested in estimating the parameter values that led to the observed data. Thus, our primary interest can be said to lie in one conditional probability above all others: p p a r a m e t e r | d a t a . That is, we are interested in the probability that a specific parameter estimate is the true population value given the data we have observed. By translating Bayes’s original equation into statistical terms, we can calculate this value:

p ( p a r a m e t e r | d a t a ) = p ( d a t a | p a r a m e t e r ) p ( p a r a m e t e r ) p ( d a t a ) .

2

The Logic of Using Functions

Using this version of Bayes’s theorem, one could hypothetically calculate the probability for any single estimate of a parameter. However, it would be much more useful if we could use this equation to simultaneously calculate the probability of each possible parameter value; obtaining the probability for each potential value would be extremely informative and allow one to make highly informed statistical inferences. Bayesian estimation is a means to this end and is therefore often described as a remarkably rich source of information on which to base inferences (Kruschke, 2013; Kruschke et al., 2012). Mathematically, this is accomplished through a revolutionary step: replacing the point probabilities in Bayes’s original formula with entire probability distributions that represent the entire range of potential parameter values. A probability distribution is simply a mathematical function: a relation that produces an output for a given input. For instance, f x = x 2 is a function where the relation between the inputs (the values on the x-axis) and outputs (the values on the y-axis) is the square of each x-value (i.e., it is a function of x). What makes a probability distribution unique is that its outputs on the y-axis are probabilities, and in Bayesian statistics, the inputs on the x-axis are the range of possible parameter values. Thus, using functions allows the analyst to assess the probability of each possible value simultaneously.

Users of frequentist methods are assumedly familiar with several probability distributions, including the univariate normal, Student’s t, and F-distributions, which are commonly used for testing the statistical significance of a parameter estimate. However, they are perhaps less familiar with the underlying mathematical functions that govern their graphical shape. For example, the form of the univariate normal distribution is governed by the following function,

f x = 1 2 π σ 2 e x p − ( x − μ ) 2 σ 2 2 ,

3

and the F-distribution by this function:

f x = ( d 1 x ) d 1 d 2 d 2 ( d 1 x + d 2 ) d 1 + d 2 x B ( d 1 2 , d 2 2 ) .

4

It is important to note that, although these functions appear complex and markedly different, one does not work with them directly when doing Bayesian estimation; as in frequentist procedures, the actual mathematics are performed by the software. However, it is important to understand conceptually that Bayesians use probability distributions for statistical modeling and that these all serve the same purpose: to return a measure of probability for each input value on the x-axis.

Because we can calculate the relative probability for all possible parameter values by using the functions that underlie probability distributions, Bayes’s formula originally defined for single probabilities can be reexpressed as,

f ( B | A ) = f ( A | B ) f ( B ) f ( A ) ,

5

or in statistical terms,

f ( p a r a m e t e r | d a t a ) = f ( d a t a | p a r a m e t e r ) f ( p a r a m e t e r ) f ( d a t a ) .

6

An added convenience is that when using functions Bayes’s formula reduces to

f ( p a r a m e t e r | d a t a ) ∝ f d a t a | p a r a m e t e r × f p a r a m e t e r , 1

7

where the “∝” symbol means “proportional to” and can be interpreted as an equals sign for the purposes of this article. This simplified version of Bayes’s theorem is the form that is explicitly used in Bayesian estimation; every Bayesian is highly familiar with it, and its three terms represent the posterior, likelihood, and prior, respectively. The posterior distribution on the left side of the equation is the ultimate goal of Bayesian estimation: the probability distribution of the parameters given the data that have been observed. The function underlying this probability distribution is simply the product of the two functions on the right that correspond to the two different sources of information that comprise these updated beliefs. The first term, f ( d a t a | p a r a m e t e r ) , is the likelihood, a function that captures the information provided by the data regarding the plausibility of different parameter values. The second term, f ( p a r a m e t e r ) , represents the prior distribution and is a function that formalizes one’s prior knowledge regarding the parameters.

The Posterior as a “Compromise”

As noted above, the function that underlies the posterior distribution is simply the product of those that underlie the likelihood and prior (i.e., the two contributing sources of information: data and prior knowledge). Therefore, it is important to understand that one only needs to specify a particular probability distribution for the likelihood and prior, whose underlying functions, when multiplied together, are the posterior. In other words, one does not specify a probability distribution for the posterior directly; it naturally results from the combination of these “two sources.” It is therefore useful to think of the Bayesian posterior as “a point that represents a compromise between the prior information and the data, and the compromise is controlled to a greater extent by the data as the sample size increases” (Gelman et al., 2004, p. 36). That is, even if we start out with some very strong prior beliefs, heavily weighing particular parameter values over others, a large amount of contrary data will overwhelm them, and the resulting posterior distribution will be closer in form to the likelihood. It is especially helpful to understand this relationship graphically. Figure 1 depicts how the shape of the posterior distribution for a parameter (here, a probability) changes under three different combinations of likelihoods and priors. The upper panel illustrates a case where the analyst has specified strong prior beliefs regarding the probability of certain parameter values: Values near 1 are heavily favored. However, the likelihood function—based on how the data “speak” to the probability of the different parameter values—strongly contradicts these prior beliefs. Because the posterior can be viewed as a “compromise” between the two functions, its form and location reflect the influence of both. However, it is ultimately much closer in form to the likelihood, as the large number of observations in the likelihood caused it to be weighted more heavily in the analysis relative to the prior.

OPEN IN VIEWER

Figure 1.

Three examples of the likelihood and prior combining to form the posterior distribution. Note that in the bottom panel, the posterior is identical to the likelihood function, as the flat noninformative prior exerts no influence in the analysis.

The center panel illustrates how the same data (i.e., the likelihood) interacts with a prior distribution that is notably more similar. The resulting posterior is very similar to that of the first plot due to the influence of the data, but importantly, is less spread out (i.e., it has a reduced variance). This results from the fact that, relative to the first plot, the likelihood and prior contain information that is much more consistent, manifesting in less uncertainty (i.e., spread) within the posterior distribution.

Finally, the bottom panel displays a “flat,” noninformative prior distribution that weighs each possible parameter value equally so that the posterior is entirely dominated by the data in the likelihood function. In fact, because the prior distribution asserts no influence in the analysis, the posterior is identical to the likelihood and is completely overlapped by it in the figure.

As can be seen, it is the product (i.e., combination) of the likelihood and prior that determines the underlying function and corresponding graphical shape of the posterior—the goal of Bayesian estimation. In fact, because in a mathematical sense they are the posterior, specifying these terms essentially is the process of Bayesian modeling. It is therefore crucial to understand how these functions are mathematically constructed—the process of how specific probability distributions are selected for these terms in Bayes’s theorem. Depending on the specific analysis at hand, the distributions might be those cited above (the normal or F-distribution) or other, perhaps less familiar, forms (e.g., the Poisson, binomial, or multivariate normal distributions). Specifying an appropriate distribution for the likelihood and prior in Bayes’s theorem (Equation 7) will be a central topic in the forthcoming sections, where each of these two sources is described in turn.

Addendum: Probability Versus Probability Density

It is important to note that in the plots in Figure 1 and in Bayesian estimation in general, the y-axis does not usually represent probability per se. This is because for continuous probability distributions, the y-axis represents probability density, and not probability proper. The reason for this is that a continuous distribution contains an infinite number of continuous values on the x-axis. For instance, the normal distribution ranges from –∞ to +∞, but given any two particular values (e.g., 0 and 1) there are an infinitude of numbers between them (.1, .01, .001, etc.). This fact makes the true probability of any one specific value infinitesimally small. Thus, probability density is used as the metric of probability in these distributions which is defined as the amount of probability within a particular interval divided by its width. In other words, it is the “density” of probability within a specific interval along the x-axis. This is why the continuous probability distributions seen in Figure 1 allow for raw y-axis values greater than 1 (i.e., greater than unity). However, for discrete probability distributions the value on the y-axis is, in fact, probability properly conceived. Intuitively, a discrete distribution is one whose x-axis contains a number of distinct values (e.g., integers only) which allows each value to be allotted its own value of probability. Though the probability versus probability density distinction may seem to be a hindrance and add a layer of unnecessary complexity, the practical consequences are actually minimal; regardless of whether continuous or discrete distributions are used, the prior, likelihood, and posterior still provide accurate measures of relative probability, which is what the analyst is interested in—that is, which parameter values are the most likely relative to the others (and by how much) given the observed data. For a more thorough discussion regarding this distinction see Kruschke (2011, chap. 3).

Selecting a Distribution for the Likelihood

Because the likelihood (the first term in Bayes’s theorem) is a function that models the data, the first step in its construction is selecting a probability distribution that is appropriate to the dependent variable(s) being modeled. That is, one must select a distribution that plausibly represents the underlying population that generated the observed data. For instance, if one is undertaking a multivariate analysis with multiple dependent variables, then a multivariate distribution must be chosen for the likelihood (e.g., the multivariate normal or multinomial distributions), as multivariate observations cannot come from a univariate distribution. Similarly, if the dependent variable is categorical, then a discrete probability distribution is most appropriate, rather than one that is continuous. Here, the goal is to select a specific distribution that is most accurately able to model the population from which the sample came.

Evidently, what should determine this selection is the actual nature of the values being modeled. Therefore, when initially specifying the distribution for the data in the likelihood function, it is always useful to first ask oneself several questions, such as, “Are my observations univariate or multivariate?” “Are they continuous or discrete?” and “Are my data restricted to certain values?” (e.g., integers, nonnegative values). A useful heuristic is to view this selection as, in part, a process of elimination, filtering out distributions that are clearly not suitable as population models for the data. Furthermore, data often come in one of several canonical forms (e.g., continuous or dichotomous values), and there are standard choices for what distribution to select for a particular data type. Accordingly, Table 1 lists different common forms of data, the distributions that are conventionally chosen as their population model for the likelihood, and the common statistical analyses associated with each. Appendix A provides intuitive conceptual descriptions of many different distributions and how they are conventionally used in Bayesian estimation.

OPEN IN VIEWER

Table 1.

Standard Distributions for Modeling Different Types of Data in the Likelihood Function.

Type of Dependent Variable	Standard Distribution	Continuous or Discrete	Associated Analyses
Continuous	Normal	Continuous	Univariate analyses (e.g., linear and curvilinear regression; t test; ANOVA; multilevel models)
Continuous	Student’s t	Continuous	Robust univariate analyses
Dichotomous	Bernoulli or binomial	Discrete; discrete	Logistic regression
Ordinal	Normal or categorical	Continuous; discrete	Ordinal regression
Count	Poisson	Discrete	Poisson regression
Multivariate continuous	Multivariate normal	Continuous	Multivariate analyses (e.g., multivariate regression, MANOVA, latent variable models)
Multivariate continuous	Multivariate Student’s t	Continuous	Robust multivariate analyses
Multivariate categorical	Multinomial	Discrete	Multivariate analyses (e.g., latent variable models)

Note: Although these represent conventional selections for modeling the data, the reader should also be aware that alternative models may be used to address the same research question, just as in the frequentist paradigm.

Example of a Likelihood Function: Estimating a Population Mean

To make the process of modeling the data in the likelihood function more concrete, suppose that we wish to perform Bayesian estimation on a particular parameter, say, the mean number of past deviant work behaviors performed by job applicants. Suppose further that we have collected data from 100 random applicants, where each observation represents their number of prior deviant behaviors. Thus, each observation is an integer, and we are therefore modeling count data (i.e., integer data). Examining the nature of this variable using the preceding questions (e.g., “Is it continuous or discrete?”), we know that we require a probability distribution that is univariate, discrete, and ideally bound at 0 (because no individual can have a negative number of deviant behaviors). Looking at Table 1, the standard distribution for modeling count data is the Poisson distribution, and from the descriptions in Appendix A, we know that it is a probability distribution that fits each one of these desired characteristics. The following mathematical function underlies its form:

p x = e − λ λ x x ! .

8

In this function, x represents an observed data value, like every probability distribution, and the variable λ is a parameter—the mean of the distribution and the parameter of interest. Substituting a data value for x and leaving λ as a free parameter will cause this function to return the probability of this particular observation under each possible value of λ, as discussed in the preceding section. For instance, an observation of 3 (as a value of x) will not be very probable given a population mean (λ) of 100. Intuitively, this observation is most probable given values of λ nearest to 3, while the probability of this observation decreases as the value of λ becomes further removed. To illustrate this concept, the top panel of Figure 2 displays a Poisson likelihood function when the observed data value is 3. In this way, the data “speak” as to which specific parameter values are more or less likely. Finally, now that a parameter has been specified in the model, Bayes’s general theorem can be rewritten to reflect our current working example:

f ( λ | d a t a ) ∝ f d a t a | λ × f λ .

9

OPEN IN VIEWER

Figure 2.

Three Poisson likelihood functions with different numbers of observed values that make different values of the parameter, λ (the mean), more or less likely. The dashed line represents the point in the function that corresponds to the most likely value of λ given the data (3, 4, and 1.86, respectively).

Constructing a Full Likelihood

This is how a single observation is formalized into knowledge regarding a parameter. First, a probability distribution is selected that represents a plausible population model that generated the observation. Then, the observed value is substituted for x within its underlying function to yield the probability of that specific value under the range of possible parameter values. However, more generally we wish to know how probable the data set is as a whole relative to the range of possible parameter values; the function constructed from a single observation will not be very informative. If a single observation is modeled by a probability distribution, how then do we combine the information contained in the entire data set? The answer to this question lies in basic probability theory, where a fundamental axiom is that the probability of two independent events occurring is merely the product of their individual probabilities: For example, the probability of events A and B both occurring is p e v e n t A × p ( e v e n t B ) . This rule can be extended to the functions underlying probability distributions as well. That is, we can insert each independent data value within the particular function we have selected and then multiply them all together, just as when working with the point probabilities above. For instance, if we have selected a Poisson distribution as a suitable population model for the data, and there are only two observations in the data set, say, 3 and 5, then the full likelihood function for λ will be the following equation:

f d a t a | λ = e − λ λ 3 3 ! × e − λ λ 5 5 ! .

10

This equation is nothing more than two Poisson distribution functions multiplied together with each observation substituted as x. However, in actual research, data sets normally contain more than two observations, and so a Poisson likelihood for an entire data set can be expressed more generally as

f d a t a | λ = e − λ λ x 1 x 1 ! × e − λ λ x 2 x 2 ! × … × e − λ λ x N x N ! ,

11

where x₁ , x₂ , and x_N represent the individual observations, and λ is the population mean that is left as the free parameter on the x-axis. Placing this within the context of our example, suppose that the first three individuals in our fictitious data set (N = 100) reported 3, 5, and 4 prior instances of deviant workplace behavior, respectively. In constructing the likelihood function, we would substitute each data value as x in the particular function we have chosen (for us, the one underlying the Poisson distribution) and then multiply each together:

f d a t a | λ = e − λ λ 3 3 ! × e − λ λ 5 5 ! × e − λ λ 4 4 ! × … e − λ λ x 100 x 100 ! .

12

Equation 12 represents the full likelihood function for our example, with the remaining 97 observations excluded for brevity. Appendix B shows how Bayes’s theorem has been updated now that the likelihood has been fully defined.

The total number of observed values also influences the likelihood function. The center panel of Figure 2 displays a Poisson likelihood that results from having only three observed values (x = 3, 5, 4) in the data set. Notice that relative to the top panel, which displays the likelihood function from a single observation (i.e., 3), when more observations are included the spread (i.e., width) of the likelihood function is reduced. Thus, with accumulating consistent information, our uncertainty is reduced regarding which parameter values led to the collected data; each observation is more or less probable under different parameter values (here, λ, the mean) and contributes to the final form of the likelihood function. The bottom panel of Figure 2 displays a Poisson likelihood with four additional observations: 0, 0, 1, 0. Because these values are lower than the previous observations, the likelihood has shifted over to the left, reflecting the fact that a lower value of λ is more likely to lead to this entire set of data. The spread (variance) is also significantly reduced, as with accumulating information, λ values greater than 4 are increasingly unlikely.

Summarizing the preceding discussion, the following represents the general process of constructing a likelihood function. First, one assesses the nature of the dependent variable(s) and selects a probability distribution that appropriately reflects these characteristics to serve as a model for the population that generated each observation. Then, each individual observation takes the place of x within the function underlying the chosen distribution and each function is multiplied together to yield the full likelihood. Although our example focused on a Poisson likelihood, the very same principles apply for the construction of all likelihood functions and can be simply extrapolated accordingly (e.g., if a normal distribution is chosen, then the likelihood will be comprised of N normal distribution functions multiplied together). Once these steps have been accomplished, the first term in Bayes’s theorem, f ( d a t a | p a r a m e t e r ) , will have been specified. However, before the function underlying the posterior distribution is fully defined (the final goal of Bayesian estimation), the other term in Bayes’s theorem, the prior distribution— f ( p a r a m e t e r ) —must specified as well. Similar to the likelihood function, it too will be modeled by a specific probability distribution, but does not provide information regarding the parameters based on data. Rather, its source of information is the analyst’s prior beliefs.

Selecting a Prior Distribution

Irrespective of how much information is included in the prior, the first step in specifying this part of the Bayesian model is to select a probability distribution that appropriately reflects the properties of the parameter to be modeled—in much the same way the data are modeled by a distribution in the likelihood function. Returning to our example, the parameter of interest is λ, the mean number of deviant behaviors. What are the desirable characteristics of a distribution for modeling this parameter? The questions used for assessing the nature of the data can also be used here: For example, “Is this parameter univariate or multivariate?” “Continuous or discrete?” “Bound to a certain range?” To start, this particular parameter is univariate and continuous, as means are not limited to integers. We also know that this particular mean cannot be negative, as no individual can have a negative number of deviant behaviors. Selecting a prior distribution for this parameter is made easier by the fact that, just like when modeling data, there are conventional choices in Bayesian statistics for modeling parameters. Table 2 lists typical prior distributions for the parameters in different likelihood functions. Looking at this table, the most common prior distribution for modeling λ, the mean of the Poisson distribution, is the gamma distribution. This distribution is univariate, continuous, and bound at zero—all desirable characteristics for modeling this particular parameter.

OPEN IN VIEWER

Table 2.

Common Selections for Prior Distributions and Their Hyperparameters.

Likelihood Distribution	Associated Parameters in Likelihood	Standard Prior Distributions for Each Parameter	Hyperparametersa Associated With Each Prior
Normal	1. μ: location parameter (mean) 2. σ2: scale parameter (variance)	1. Normal 2. Inverse gamma	1. μ: location parameter (mean); σ2: scale parameter (variance) 2. α: shape parameter; β: scale parameter
Student’s t	1. μ: location parameter (mean) 2. σ2: scale parameter (variance) 3. ν: degrees of freedom	1. Normal 2. Inverse gamma 3. Uniform	1. μ: location parameter (mean); σ2: scale parameter (variance) 2. α: shape parameter; β: scale parameter 3. a: location parameter; b: scale parameter
Bernoulli	1. p: probability of event (mean)	1. Beta	1. α: shape parameter; β: shape parameter
Binomial	1. p: probability of event (mean) 2. n: number of total trials	1. Beta 2. None (fixed from data)b	1. α: shape parameter; β: shape parameter 2. None
Poisson	1. λ: both a location and scale parameter (mean and variance)	1. Gamma	1. α: shape parameter; β: rate parameter

Note: These priors represent many of the classic choices for modeling the parameters of likelihood functions. However, there are additional distributions that can be used, and the selection should always depend on the particular analytic context. For additional information on these distributions and others, see Gelman, Carlin, Stern, and Rubin (2004, Appendix A) and Lynch (2007, pp. 25-33).

^aHyperparameters are the parameters of prior distributions.

^bThe n parameter representing the number of total trials does not require a prior distribution to quantify uncertainty because it is given by the data.

Types of Parameters

However, the gamma distribution, like any probability distribution, can take on many different shapes. Though every distribution has its own specific, unchanging characteristics (whether it includes negative values, is continuous or discrete, etc.), there are variables in its underlying function that set its particular form. For instance, the univariate normal distribution function has two variables that define the specific form it will take: the mean (μ) and variance (σ²). Although it is always symmetrical, univariate, and continuous, varying either one of these values alters the ultimate form of the distribution: The mean determines the location of the center, and the variance determines its spread. These variables that determine the particular form of a probability distribution are themselves referred to as parameters. Although each distribution has its own unique parameters, they often perform similar functions and can usually be classified into one of three broad categories: (1) location parameters that determine where the distribution is located on the x-axis (e.g., μ, the mean of the normal distribution), (2) scale and rate parameters that define the spread of the distribution (i.e., increasing a scale parameter such as σ², the variance of the normal distribution, increases its width, while increasing a rate parameter, such as β of the gamma distribution, decreases its width), and (3) shape parameters that affect neither the location nor the width directly, but alter the shape of the distribution in a particular way (e.g., by altering its kurtosis or the width of its tails).

That said, Bayesian statistics requires not only a working knowledge of the different distributions available for modeling, but also how they are parameterized (i.e., their underlying parameters). Becoming familiar with the different types of distributions and their parameters is a necessary yet rewarding challenge for the student of Bayesian analysis. To help the beginner through this process, Appendix A contains conceptual descriptions of many distributions and their parameters, along with whether or not they are classified as a particular type (viz., location, scale/rate, or shape). Furthermore, Appendix C provides R code that enables the user to plot a number of distributions under different values of their parameters. Experimenting with various values promotes an intuitive understanding of how different parameter values affect the final graphical form of a distribution.

Hyperparameters: Defining Prior Knowledge

Returning to our example, we have decided that the gamma distribution is an appropriate choice to model the parameter of interest: the mean number of prior deviant behaviors by current job applicants. Its underlying algebraic function is

f x = β α α − 1 ! x α − 1 e − β x ,

13

where α and β are the distribution’s parameters that define its form. The parameter α is a shape parameter that controls its skewness (as α increases the gamma distribution decreases in skew and approaches normality, as explained in Appendix A), and β is a rate parameter (i.e., as its value increases, the width of the distribution decreases). Because these are higher-level parameters—that is, parameters of the prior distribution that models our beliefs regarding a parameter—they are conventionally referred to as hyperparameters. This terminology can be somewhat confusing for fledgling Bayesian analysts: Parameters in the likelihood function (e.g., the mean number of deviant work behaviors) are modeled by prior distributions whose forms are themselves determined by additional, higher-order parameters. ² However, the distinction is made because Bayesian models are hierarchical by nature, and using this terminology allows one to distinguish the parameters found at different levels. Figure 3 depicts the hierarchy of our current Bayesian model. At the bottom level, the data is modeled by the Poisson distribution which has a single parameter, λ, the mean. In turn, our prior beliefs regarding this parameter are modeled by the gamma distribution, which has its own parameters: the hyperparameters α and β.

OPEN IN VIEWER

Figure 3.

A hierarchical illustration of the current Bayesian model in the style of Kruschke (2011).

As varying either one of these variables will alter the ultimate form of this distribution, we wish to select specific values that will make it congruent to our prior beliefs regarding λ. That is, if we have prior knowledge, we may wish to make use of this by weighing certain parameter values more strongly at the outset of the analysis and fix the values of the hyperparameters accordingly. A suggested method of doing this is by plotting multiple distributions with the parameters varied each time. This allows one to see how the form of the distribution changes and promotes an intuitive understanding of the intimate relationship between its shape and its parameter values. Appendix C contains further R code for plotting overlaid distributions for comparison purposes.

We now return to our running example of estimating the population mean (λ) of the number of past deviant work behaviors for job applicants. Despite this being merely a hypothetical example, previous research could potentially inform the analysis. Dilchert, Ones, Davis, and Rostow (2007) reported the frequency of officially documented deviant work behaviors by individuals applying for police officer jobs. Approximately 77% of individuals had no officially documented cases, 9% had one case, 5% had two, and 7% had three or more. Other studies have found similar general rates of deviant behaviors (e.g., Bolton, Becker, & Barber, 2010; Mount, Ilies, & Johnson, 2006). This information can be used to define our prior beliefs regarding which parameter values are more probable relative to others. Based on this information, the most likely values for λ are quite low (the mean of the Dilchert et al., 2007, sample was around 0.4), so we would like to weigh these values more heavily. Weighing low values relative to higher ones entails that the distribution be significantly right-skewed. We therefore initially set the shape parameter that determines the skewness of the gamma distribution, α, to a very low value, say, 1. (As noted in Appendix A, the gamma distribution is increasingly skewed as this parameter is reduced.) We also would like the distribution to be somewhat wide to capture our uncertainty regarding this parameter; we have previous evidence for a low mean, but weighing certain values too heavily can influence the shape of the posterior more strongly than we would like. We therefore initially set the rate parameter β to another low value: 0.5. The form of this distribution is shown in the top panel of Figure 4. As can be seen, values below 2 are weighted very strongly with the probability density of larger values decreasing quickly, accurately reflecting the prior research findings as well as reflecting some of our subjective uncertainty regarding the true population value.

OPEN IN VIEWER

Figure 4.

An informative gamma prior distribution based on previous research findings and a noninformative gamma prior distribution for the same hypothetical example.

Our prior beliefs have now been formalized by this probability distribution: It strongly favors low values, while those that are farther away are increasingly less probable. As stated earlier, a common practice in Bayesian modeling is to use a noninformative priors to let the data fully dominate the form of the posterior distribution (e.g., the bottom panel of Figure 1). These are much simpler to specify because (a) they do not require precise calibration of the hyperparameters to reflect specific prior knowledge and (b) there are conventional noninformative priors that are used to model various parameters (e.g., the uniform distribution discussed in Appendix A). Often, specifying a noninformative prior simply requires increasing the spread of the distribution so that all plausible parameter values are weighted equally. This can often be done by either increasing the value of a scale parameter or, equivalently, by decreasing the value of a rate parameter (which, again, has the inverse effect of decreasing a scale parameter). For example, if one has chosen a normal distribution as a prior distribution for a parameter (e.g., a mean or regression slope), then setting its variance to a very large number (e.g., 10,000) will result in a very wide distribution where all plausible values are weighted equally. This either reflects the analyst’s total lack of knowledge regarding the potential parameter values or their desire to let the data alone determine the form of the posterior. Returning to our example, a noninformative gamma distribution is commonly specified by significantly reducing both the shape and rate hyperparameters to 0.001 (see Spiegelhalter, Thomas, Best, & Lunn, 2003). This noninformative gamma prior is shown in the bottom panel of Figure 4. Although this distribution seems to “spike” as it approaches zero, it is essentially flat for all plausible values of λ.

Now that both functions in Bayes’s theorem—the likelihood and prior—have been fully specified by different probability distributions, the posterior distribution of our example is defined. Again, the posterior distribution provides us with our conditional probability of interest, f ( p a r a m e t e r | d a t a ) , the probability of all potential parameter values given the observed data, and is itself a function—one that is the product of the likelihood function and the function underlying the prior distribution. (See Appendix B for the full mathematical function that represents the posterior in our example.) Furthermore, because the posterior represents a “compromise” between these two terms, their combination will shift our beliefs away from the prior to the extent that is supported by the data, a highly attractive feature of Bayesian estimation.

Summary

We close this section by offering a general summary of the process of specifying a prior distribution. First, each parameter in the model must be defined by a prior distribution that reflects the amount of previous knowledge the analyst wishes to incorporate within the analysis. Selecting an appropriate distribution to model parameter is a similar process to modeling data in the likelihood function: The nature of the values to be modeled is the primary focus (e.g., discrete or continuous, bound at certain intervals). Once a suitable distribution has been selected, the parameters of this prior distribution (i.e., its hyperparameters) must be specified to reflect the amount of desired information included within the prior. Often, flat, noninformative priors are used for simplicity and to forestall potential criticisms regarding this “subjective” component (Howson & Urbach, 1991, p. 372; see also Gelman, 2008). However, if reliable information is available, this can and should be used to increase the precision of the analysis.

Informal Methods

Assessment of Markov chain algorithm performance (viz., convergence and mixing) is accomplished through both formal and informal methods. It is important to note that there is no single, definitive method that guarantees that proper MCMC sampling has occurred (Lynch, 2007, p. 135). However, a combination of both formal and informal methods can provide sufficient evidence for proper MCMC performance. Despite being more subjective, informal methods need not be any less rigorous than those that are formal. In fact, these are often the best methods for diagnosing good or poor chain performance. Primarily, this entails examining trace plots. A trace plot displays the value of each sampled parameter value (y-axis) at each step in the chain (x-axis) and can provide evidence for both convergence and mixing. Convergence is evinced when the chains appear to “settle” in the same general region, indicating that they have reached the central regions (i.e., the bulk) of the distribution. Mixing is indicated when the chains exhibit enough variability (i.e., vertical movement) around this region to sample properly from the entire posterior. Therefore, each chain ought to not only consistently sample around a general region (its bulk), but also display variability around it. Figure 5 depicts two trace plots of different parameters whose posteriors have each been sampled by two chains started at different values. The two chains in the upper plot (represented by the two lines) have clearly not reached convergence, as they are completely isolated; converged chains will always overlap, as they sample values from the same general region of the posterior. In contrast, the chains in the lower plot appear to have converged almost immediately, as both have settled around the same general region and are consistently sampling comparable values. They also demonstrate good mixing as there is adequate variability around this region. To ensure proper estimation, thousands of samples should be taken regardless of how quickly the chains seeming to converge. Note that these plots only show the first 400 values sampled by the chains. This should emphasize how variable convergence rates can be across different models and even for parameters within the same model. (The two parameters displayed here are from the same Bayesian model.)

OPEN IN VIEWER

Figure 5.

An example of trace plots generated by WinBUGS that show the value of the sampled parameter value (y-axis) at each step in the chain (x-axis). Both plots contain two Markov chains each represented by a different line.

Formal Methods

Formal methods for diagnosing convergence and mixing involve the calculation of statistics. When multiple chains are used, the Brooks–Gelman–Rubin convergence diagnostic (BGR; also called the potential scale reduction factor; Gelman & Rubin, 1992) is a common method for assessing whether or not the Markov chains have reached convergence (i.e., the bulk of the posterior). Essentially, this statistic partitions the variance of the sampled parameter values into its between and within-chains variance. Significant between-chain variance is undesirable because it indicates that the different chains are sampling notably different parameter values and implies that at least one chain has not yet reached convergence. When the chains have converged, there should be very little between-chain variance because they will be sampling similar values. Conversely, within-chain variance is desirable, as it indicates mixing—that is, there is variability in the samples drawn by the chains to represent the entire posterior. The BGR convergence diagnostic is simply a ratio of the total variance (both within and between-chains) to the within-chains variance of the sampled parameter values. Thus, when this ratio approaches 1 (e.g., < 1.10), this indicates that there is little between-chain variance and provides evidence that the chains have converged.

A formal method for determining if proper mixing has occurred, and also if enough samples have been collected, is to calculate the Monte Carlo error for each parameter. This statistic assesses the accuracy of the summary statistics that would be calculated from the current set of samples (Spiegelhalter et al., 2003). A general rule is that the number of collected samples is sufficient if the Monte Carlo error is below at least 5% of the standard deviation of the set of samples drawn for each parameter (Spiegelhalter et al., 2003).

The aforementioned methods represent the primary tools for assessing MCMC performance; identifying convergence and mixing (or lack thereof) are often accomplished through these diagnostics. However, other supplementary methods exist, and for a more additional information regarding MCMC diagnostics the reader is directed to Lynch (2007, chap. 6).

Model Specification and MCMC Diagnostics

When reporting a Bayesian analysis, there are several components that must be communicated (Anderson, Link, Johnson, & Burnham, 2001; Kruschke, 2011, chap. 23): One must first describe the model for the data, the choice of priors, and the process and diagnostics of the MCMC sampling (see Table 3). In our analyses, the dependent variable was dichotomous and modeled in the likelihood function by the Bernoulli distribution. The intercept and slope parameters were given noninformative priors due to the lack of previous work in this subject area. They were modeled by the univariate normal distribution with a mean (μ) of 0 and variance (σ²) of 10,000 to ensure a very wide distribution so that all plausible parameter values were weighted equally at the outset of the analysis. The first 10,000 MCMC samples were discarded in the burn-in period, and posterior inferences were made after an additional 20,000 samples were collected. The Monte Carlo error for each parameter was less than 5% of its standard deviation (range: 0.01-1.28%), indicating that enough samples had been collected for accurate estimation. Both regression models converged very quickly, as evidenced by trace plots and a sufficiently low BGR convergence diagnostic. (For all parameters it was < 1.07 by the end of the burn-in period.) Finally, proper mixing was also evidenced by trace plots and the low Monte Carlo error for each parameter.

The Current Results

Before making Bayesian inferences, it is useful to first visually assess the posterior distribution approximated by the set of MCMC samples so that appropriate and theoretically relevant descriptions can be made. We therefore begin our analysis by looking at the posterior distributions of the three regression slopes that were statistically significant in the frequentist analysis: those relating hobo syndrome to contractual forces, behavioral forces, and openness to experience, displayed in Figure 7. The BugsXLA tool for posterior plotting juxtaposes a parameter’s prior distribution with the approximated posterior, showing how one’s initial beliefs have been updated by the observed data (i.e., the “Bayesian process” discussed in Part 1). Because our analyses used flat, noninformative priors, the posterior effectively represents the likelihood function—that is, how the data differentially weigh candidate parameter values. It is also important to note that the y-axis in these plots does not display probability proper because these distributions only approximate the posterior through the set of samples drawn by the Markov chains. In other words, it is a histogram that, given proper convergence and mixing, will accurately represent the posterior, the “relative frequencies” of the y-axis indicating how probable a particular parameter value is relative to the others (i.e., it is functionally equivalent to measures of relative probability).

OPEN IN VIEWER

Figure 7.

BugsXLA posterior plots of the contractual, behavioral, and openness slopes. The black line at the bottom of the distribution is the “flat” noninformative prior, and the darker area constitutes the region 0.90-1.10, which was defined as practically equivalent to a null effect (the ROPE).

Because we are interested in the presence of practically significant effects, we defined the region of practical equivalence (the ROPE) as the range 0.90 to 1.10. That is, for a predictor to be “practically significant,” a standardized one-unit increase must make the odds of hobo syndrome increase or decrease by at least 10% (i.e., an odds ratio greater than 1.10 or less than 0.90, respectively). The approximate region of the ROPE is depicted by the darker shade in the posterior plots. BugsXLA has a function to retrieve the probability that the parameter value is above or below a specific value in the distribution. Once the ROPE has been defined, one can simply use this function to yield the probability that the parameter lies outside or within this region.

In deciding which point estimates to use, the summary statistics of Table 4 indicate that the means and medians for each parameter are virtually identical (i.e., the posteriors are approximately normal), so using any point estimate will yield a similar conclusion. Just as in the frequentist analysis, the contractual slope parameter was the strongest predictor of hobo syndrome. It had a posterior mean of 3.79, a 95% credible interval (CI) of [2.19, 6.82], and 100% of its distribution was located above the ROPE. Similarly, the regression slope for behavioral forces had a mean value of 1.60, 95% CI [1.09, 2.38], and 97.33% of the distribution was located above the ROPE, indicating that low perceived costs of leaving are positively related to hobo syndrome. Examination of the openness posterior (mean: 1.19, 95% CI [.86, 1.64]) also indicated a positive relationship, with 67.80% of the posterior located above the ROPE. Despite the fact that this posterior is more ambiguous, the majority of the evidence still supports a positive relationship, even if the parameter value is most likely smaller than the previous two. Had only 50% been located above the ROPE, then it would have been prudent to abstain from making substantial conclusions regarding the presence of an effect. Again, because the posterior contains such a large amount of information, the analyst must be able to make flexible assessments that are appropriate given the nature of the present evidence.

Expectedly, these posterior estimates were similar to the maximum likelihood estimates of the frequentist models because, when noninformative priors are used, Bayesian and frequentist estimators often coincide (Bayarri & Berger, 2004, p. 63; Gelman et al., 2004, p. 363). In research settings where the use of prior knowledge is justified, Bayesian estimation can differ markedly from maximum likelihood estimation; the former can be more precise as the latter estimation procedure is unable to incorporate previous findings. Moreover, even in the current situation where noninformative priors were used, the posterior speaks directly on the probability of all possible parameter values and yields point estimates that are unsuppressed by NHST.

In the frequentist regression, the vast majority of predictors (six out of nine) failed to reach statistical significance at α = .05. This leads to ambiguous results (“approaching significance”) and an inability to make inferences. However, a considerable advantage of Bayesian estimation is that even when estimates are more ambiguous (which is often the case in research), inferences can still be made, and these estimates can inform subsequent research. The posteriors of the three “marginally significant” slopes of the frequentist analysis (p > .05 but < .10), are presented in Figure 8.

OPEN IN VIEWER

Figure 8.

BugsXLA posterior plots of the affective, alternative, and constituent slopes.

Here, the slope relating affective forces to hoboism had a posterior mean of 0.64, 95% CI [0.41, 1.00], and 93.75% of its posterior lay below the ROPE, indicating the existence of a practically significant negative relationship. The alternative forces slope had a posterior mean of 1.40, 95% CI [0.97, 2.08], and over 89.78% of its distribution supported a positive relationship. Finally, the constituent forces slope had a mean of 1.48, 95% CI [1.00, 2.25], and over 93.56% of its posterior was also located above the ROPE. Though these point estimates were similar across both paradigms, unconstrained by NHST, Bayesian data analysis allows for substantive inferences that are essential to understanding the topic of interest. In contrast to the original predictions, hobos were predicted by higher levels of positive affect toward the organization. Thus, these results suggest that hobos do not quit because they are averse to the organization. In fact, they may actually view the organization more favorably because they do not feel restricted by it, an explanation consistent with the very large contractual forces slope (i.e., that they feel little contractual obligation to stay). Practically, this is a crucial insight because it suggests a fundamental limit to the extent to which organizations are responsible for the chronic turnover exhibited by this group—that is, it may not be a function of job satisfaction.

The analysis also indicated that hobos are sensitive to the job market; they are more likely to believe outside alternatives are available and perceive them as desirable (alternative forces). They also display reduced personal attachment to individuals or groups within the organization (constituent forces). Their patterns of frequent turnover and positive attitudes toward quitting perhaps inhibit motivation to create lasting workplace relationships. Alternatively, this finding may indicate features of their general personality and social functioning: Hobos may feel low attachment to others and groups in general, organizations being one particular manifestation of this pattern of interpersonal behavior. In any case, this is a novel research finding that was suppressed using conventional frequentist practices. Through Bayesian estimation, these “marginally significant” effects have provided novel insights into the underlying cognitive and affective processes of job hobos.

As for the “nonsignificant” slopes of the frequentist analysis, the distributions of the calculative forces, normative forces, and impulsivity regression slopes are displayed in Figure 9. The posterior mean of the calculative slope was 1.26, 95% CI [0.71, 2.27], and 67.62% of its posterior was located above the ROPE. Thus, much like the openness posterior, the majority of the evidence supports a practically significant positive relationship. The posterior mean of the normative slope was 1.18, 95% CI [0.81, 1.73], and 64.82% of its posterior was located above the ROPE, also supporting a positive relationship. For the impulsivity slope, the posterior mean was 0.96, 95% CI [0.70, 1.32], and its posterior was located largely within the ROPE (46.32%), with the remaining probability split on either side (below 32.84%; above 20.84%).

OPEN IN VIEWER

Figure 9.

BugsXLA posterior plots of the calculative, normative, and impulsivity slopes.

As expected, these slopes were not as large as those that had lower p values in the frequentist analysis. The data here were more ambiguous, and we are inevitably less confident when drawing conclusions about these parameters. Though nonnegligible percentages of the posterior were located within the ROPE, they still provide support for two practically significant effect sizes. Certainly, further research is warranted, but tentative inferences can still be made.

From our results, it is probable that hobos feel less external pressure from individuals outside the organization (e.g., family members) to retain their job (normative forces). This lends further support to the idea that hobos have weaker personal attachments in general. Alternatively, they may indeed have these relationships but simply assign lower importance to these perceived obligations. With some confidence, we can also state that hobos perceive lower probabilities of attaining long-term career goals at their current organization (calculative forces). The weight of the evidence (46.32%) supported that impulsivity is not practically related to hobo classification. However, there is still some ambiguity as this percentage is not particularly large, and 32.84% of the posterior supported a negative relationship. What is more important is that our uncertainty in this set of analyses is not a detriment that is systematically ignored, but organized information that supports tentative inferences and informs future research in an unbiased way.