12.1.1 Region of practical equivalence

A region of practical equivalence ($\mathrm{R}$OPE) indicates a small range of parameter values that are considered to be practically equivalent to the null value for purposes of the particular application. For example, if we wonder whether a coin is fair, for purposes of determining which team will kick off at a football game, then we want to know if the underlying bias in the coin is reasonably close to 0.50, and we don't really care if the true bias is 0.473 or 0.528, because those values are practically equivalent to 0.50 for our application. Thus, the $\mathrm{R}$OPE on the bias might go from 0.45 to 0.55. As another example, if we are assessing the efficacy of a drug versus a placebo, we might only consider using the drug if it improves the probability of cure by at least 5 percentage points. Thus, the $\mathrm{R}$OPE on the difference of cure probabilities could have limits of ±0.05. There will be more discussion later of how to set the $\mathrm{R}$OPE limits.

Once a $\mathrm{R}$OPE is set, we make a decision to reject a null value according the following rule:

For example, suppose that we want to know whether a coin is fair, and we establish a $\mathrm{R}$OPE that goes from 0.45 to 0.55. We flip the coin 500 times and observe 325 heads. If the prior is uniform, the posterior has a 95% HDI from 0.608 to 0.691, which falls completely outside the $\mathrm{R}$OPE. Therefore we declare that the null value of 0.5 is rejected for practical purposes. Notice that when the HDI excludes the$R$OPE, we are not rejecting all values within the$R$OPE; we are rejecting only the null value.

Because the $\mathrm{R}$OPE and HDI can overlap in different ways, there are different decisions that can be made. In particular, we can decide to “accept” a null value:

With this decision rule, a null value of a parameter can be accepted only when there is sufficient precision in the estimate of the parameter. For example, suppose that we want to know whether a coin is fair, and we establish a $\mathrm{R}$OPE that goes from 0.45 to 0.55. We flip the coin 1,000 times and observe 490 heads. If the prior is uniform, the posterior has a 95% HDI from 0.459 to 0.521, which falls completely within the $\mathrm{R}$OPE. Therefore we declare that the null value of 0.5 is confirmed for practical purposes, because all of the most credible values are practically equivalent to the null value.

In principle, though rarely in practice, a situation could arise in which a highly precise HDI does not include the null value but still falls entirely within a wide $\mathrm{R}$OPE. According to the decision rule, we would “accept” the null value despite it having low credibility according to the posterior distribution. This strange situation highlights the difference between the posterior distribution and the decision rule. The decision rule for accepting the null value says merely that the most credible values are practically equivalent to the null value according to the chosen$R$OPE, not necessarily that the null value has high credibility. If this situation actually arises, it could be a sign that the $\mathrm{R}$OPE is too large and has been ill-conceived or is outdated because the available data are so much more precise than the $\mathrm{R}$OPE.

When the HDI and $\mathrm{R}$OPE overlap, with the $\mathrm{R}$OPE not completely containing the HDI, then neither of the above decision rules is satisfied, and we withhold a decision. This means merely that the current data are insufficient to yield a clear decision one way or the other, according to the stated decision criteria. The posterior distribution provides complete information about the credible parameter values, regardless of the subsequent decision procedure. There are other types of decisions that could be declared, if the HDI and $\mathrm{R}$OPE overlap in different ways. We will not be pursuing those other types of decisions here, but they can be useful in some applied situations. For further discussion of the $\mathrm{R}$OPE, under somewhat different appellations of “range of equivalence” and “indifference zone,” see, for example, Carlin and Louis (2009), Freedman, Lowe, and Macaskill (1984, Hobbs and Carlin (2008), and Spiegelhalter, Freedman, and Parmar (1994).

Aside from the intuitive appeal of using a $\mathrm{R}$OPE to declare practical equivalence, there are sound logical reasons from the broader perspective of scientific method. Serlin and Lapsley (1985, 1993) pointed out that using a $\mathrm{R}$OPE to affirm a predicted value is essential for scientific progress, and is a solution to Meehl's paradox (e.g., Meehl, 1967, 1978, 1997). Meehl started with the premise that all theories must be wrong, in the sense that they must oversimplify some aspect of any realistic scenario. The magnitude of discrepancy between theory and reality might be small, but there must be some discrepancy. Therefore, as measurement precision improves (e.g., with collection of more data), the probability of detecting the discrepancy and disproving the theory must increase. This is how it should be: More precise data should be a more challenging test of the theory. But the logic of null hypothesis testing paradoxically yields the opposite result. In null hypothesis testing, all that it takes to “confirm” a (non-null) theory is to reject the null value. Because the null hypothesis is certainly wrong, at least slightly, increased precision implies increased probability of rejecting the null, which means increased probability of “confirming” the theory. What is needed instead is a way to affirm substantive theories, not a way to disconfirm straw-man null hypotheses. Serlin and Lapsley (1985, 1993) showed that by using a $\mathrm{R}$OPE around the predicted value of a theory, the theory can be affirmed. Crucially, as the precision of data increases, and the width of the $\mathrm{R}$OPE decreases, then the theory is tested more stringently.²

How is the size of the $\text{R}$OPE determined? In some domains such as medicine, expert clinicians can be interviewed, and their opinions can be translated into a reasonable consensus regarding how big of an effect is useful or important for the application. Serlin and Lapsley (1993, p. 211) said, “Admittedly, specifying [$\mathrm{R}$OPE limits] is difficult. … The width of the [$\mathrm{R}$OPE] depends on the state of the art of the theory and of the best measuring device available. It depends on the state of the art of the theory …[because] a historical look at one's research program or an examination of a competing research program will help determine how accurately one's theory should predict in order that it be competitive with other theories.” In other words, the limits of the $\mathrm{R}$OPE depend on the practical purpose of the $\mathrm{R}$OPE. If the purpose is to assess the equivalence of drug-treatment outcomes, then the $\mathrm{R}$OPE limits depend on the real-world costs and benefits of the treatment and the ability to measure the outcome. If the purpose is to affirm a scientific theory, then the $\mathrm{R}$OPE limits depend on what other theories need to be distinguished.

The $\mathrm{R}$OPE limits, by definition, cannot be uniquely “correct,” but instead are established by practical aims, bearing in mind that wider $\mathrm{R}$OPEs yield more decisions to accept the $\mathrm{R}$OPEd value and fewer decision to reject the $\mathrm{R}$OPEd value. In many situations, the exact limit of the $\mathrm{R}$OPE can be left indeterminate or tacit, so that the audience of the analysis can use whatever $\mathrm{R}$OPE is appropriate at the time, as competing theories and measuring devices evolve. When the HDI is far from the $\mathrm{R}$OPEd value, the exact $\mathrm{R}$OPE is inconsequential because the $\mathrm{R}$OPEd value would be rejected for any reasonable $\mathrm{R}$OPE. When the HDI is very narrow and overlaps the target value, the HDI might again fall within any reasonable $\mathrm{R}$OPE, again rendering the exact $\mathrm{R}$OPE inconsequential. When, however, the HDI is only moderately narrow and near the target value, the analysis can report how much of the posterior falls within a $\mathrm{R}$OPE as a function of different $\mathrm{R}$OPE widths. An example is shown at this book's blog at http://doingbayesiandataanalysis.blogspot.com/2013/08/how-much-of-bayesian-posterior.html.

It is important to be clear that any discrete decision about rejecting or accepting a null value does not exhaustively capture our knowledge about the parameter value. Our knowledge about the parameter value is described by the full posterior distribution. When making a binary decision, we have merely compressed all that rich detail into a single bit of information. The broader goal of Bayesian analysis is conveying an informative summary of the posterior, and where the value of interest falls within that posterior. Reporting the limits of an HDI region is more informative than reporting the declaration of a reject/accept decision. By reporting the HDI and other summary information about the posterior, different readers can apply different $\mathrm{R}$OPEs to decide for themselves whether a parameter is practically equivalent to a null value. The decision procedure is separate from the Bayesian inference. The Bayesian part of the analysis is deriving the posterior distribution. The decision procedure uses the posterior distribution, but does not itself use Bayes’ rule.

In applications when the Bayesian posterior is approximated with an MCMC sample, it is important to remember the instability of the HDI limits. Recall the discussion accompanying Figure 7.13, p. 185, which indicated that the standard deviation of a 95% HDI limit for a normal distribution, across repeated runs with an MCMC sample that has an effective sample size (ESS) of 10,000, is roughly 5% of the standard deviation of parameter posterior. Thus, if the MCMC HDI limit is very near the $\mathrm{R}$OPE limit, be cautious in your interpretation because the HDI limit has instability due to MCMC randomness. Analytically derived HDI limits do not suffer this problem, of course.

It may be tempting to try to adopt the use of the $\mathrm{R}$OPE in NHST because an NHST confidence interval (CI) has some properties analogous to the Bayesian posterior HDI. The analogous decision rule in NHST would be to accept a null hypothesis if a 95% CI falls completely inside the $\mathrm{R}$OPE. This approach goes by the name of equivalence testing in NHST (e.g., Rogers, Howard, & Vessey, 1993; Westlake, (1976, 1981)). While the spirit of the approach is laudable, it has two main problems. One problem is technical: CIs can be difficult to determine. For example, with the goal of equivalence testing for two proportions, Dunnett and Gent (1977) devoted an entire article to various methods that approximate the CI for that particular case. With modern Bayesian methods, on the other hand, HDIs can be computed seamlessly for arbitrary complex models; indeed the case of two proportions was addressed in Section 7.4.5, p. 176, and the more complex case of hundreds of grouped proportions was addressed in Section 9.5.1, p. 253. The second problem with frequentist equivalence testing is foundational: A CI does not actually indicate the most credible parameter values. In a Bayesian approach, the 95% HDI actually includes the 95% of parameter values that are most credible. Therefore, when the 95% HDI falls within the $\mathrm{R}$OPE, we can conclude that 95% of the credible parameter values are practically equivalent to the null value. But a 95% CI from NHST says nothing directly about the credibility of parameter values. Crucially, even if a 95% CI falls within the $\mathrm{R}$OPE, a change of stopping or testing intention will change the CI and the CI may no longer fall within the $\mathrm{R}$OPE. For example, if the two groups being compared are intended to be compared to other groups, then the 95% CI is much wider and may no longer fall inside the $\mathrm{R}$OPE. This dependency of equivalence testing on the set of tests to be conducted is pointed out by Rogers et al. (1993, p. 562). The Bayesian HDI, on the other hand, is not affected by the intended tests.

12.1.2 Some examples

The left side of Figure 9.14 (p. 256) shows the difference of batting abilities between pitchers and catchers. The 95% HDI goes from −0.132 to −0.0994. If we use a $\mathrm{R}$OPE from −0.05 to +0.05 (somewhat arbitrarily), the HDI falls far outside the $\mathrm{R}$OPE, and we reject the null, even taking into account the MCMC instability of the HDI limits.

The right side of Figure 9.14 (p. 256) shows the difference of batting abilities between catchers and first basemen. The 95% HDI goes from −0.0289 to 0.0 (essentially). With any non-zero $\mathrm{R}$OPE, and taking into account the MCMC instability of the HDI limits, we would not want to declare that a difference of zero is rejected, instead saying that the difference is only marginally non-zero. The posterior gives the full information, indicating that there is a suggestion of difference, but the difference is small relative to the uncertainty of its estimate.

The right side of Figure 9.15 (p. 257) shows the difference of batting abilities between two individual players who provided lots of data. The 95% HDI goes from −0.0405 to +0.0368. This falls completely inside a $\mathrm{R}$OPE of −0.05 to +0.05 (even taking MCMC instability into account). Thus, we could declare that a difference of zero is accepted for practical purposes. This example also illustrates that it can take a lot of data to achieve a narrow HDI; in this case both batters had approximately 600 opportunities at bat.

The plotPost function, in the utilities that accompany this book, has options for displaying a null value and $\mathrm{R}$OPE. Details are provided in Section 8.2.5.1 (p. 205), and an example graph is shown in Figure 8.4 (p. 205). In particular, the null value is specified with the compVal argument (which stands for comparison value), and the $\mathrm{R}$OPE is specified as a two-element vector with the ROPE argument. The resulting graph displays the percentage of the posterior on either side of the comparison value, and the percentage of the posterior within the $\mathrm{R}$OPE and on either side of the $\mathrm{R}$OPE limits.

12.1.2.1 Differences of correlated parameters

It is important to understand that the marginal distributions of two parameters do not reveal whether or not the two parameter values are different. Figure 12.1, in its left quartet, shows a case in which the posterior distribution for two parameter values has a strong positive correlation. Two of the panels show the marginal distributions of the single parameters. Those two marginal distributions suggest that there is a lot of overlap between the two parameters values. Does this overlap imply that we should not believe that they are very different? No! The histogram of the differences shows that the true difference between parameters is credibly greater than zero, with a difference of zero falling well outside the 95% HDI, even taking into account MCMC sampling instability and a small $\mathrm{R}$OPE. The upper left panel shows why: The credible values of the two parameters are highly correlated, such that when one parameter value is large, the other parameter value is also large. Because of this high correlation, the points in the joint distribution fall almost all on one side of the line of equality.

Figure 12.1 shows, in its right quartet, a complementary case. Here, the marginal distributions of the single parameters are exactly the same as before: Compare the histograms of the marginal distributions in the two quartets. Despite the fact that the marginal distributions are the same as before, the bottom right panel reveals that the difference of parameter values now straddles zero, with a difference of zero firmly in the midst of the HDI. The plot of the joint distribution shows why: Credible values of the two parameter are negatively correlated, such that when one parameter value is large, the other parameter value is small. The negative correlation causes the joint distribution to straddle the line of equality.

In summary, the marginal distributions of two parameters do not indicate the relationship between the parameter values. The joint distribution of the two parameters might have positive or negative correlation (or even a non-linear dependency), and therefore the difference of the parameter values should be explicitly examined.

12.1.2.2 Why hdi and not equal-tailed interval?

I have advocated using the HDI as the summary credible interval for the posterior distribution, also used in the decision rule along with a $\mathrm{R}$OPE. The reason for using the HDI is that it is very intuitively meaningful: All the values inside the HDI have higher probability density (i.e., credibility) than any value outside the HDI. The HDI therefore includes the most credible values.

Some other authors and software use an equal-tailed interval (ETI) instead of an HDI. A 95% ETI has 2.5% of the distribution on either side of its limits. It indicates the 2.5th percentile and the 97.5th percentile. One reason for using an ETI is that it is easy to compute.

In symmetric distributions, the ETI and HDI are the same, but not in skewed distributions. Figure 12.2 shows an example of a skewed distribution with its 95% HDI and 95% ETI marked. (It is a gamma distribution, so its HDI and ETI are easily computed to high accuracy.) Notice on the right there is a region, marked by an arrow, that is outside the HDI but inside the ETI. On the left there is another region marked by an arrow, that is inside the HDI but outside the ETI. The ETI has the strange property that parameter values in the region marked by the right arrow are included in the ETI, even though they have lower credibility than parameter values in the region marked by the left arrow that are excluded from the ETI. This property seems undesirable as a summary of the credible values in a distribution.

The strange property of the ETI also leads to weirdness when using it as a decision tool. If a null value and $\mathrm{R}$OPE were in the region marked by the right arrow, it would be rejected by the HDI, but not by the ETI. Which decision makes more sense? I think the decision by HDI makes more sense, because it is saying that the values outside its limits have low credibility. But the decision by ETI says that values in this region are not rejected, even though they have low credibility. The complementary conflict happens in the region marked by the left arrow. If a null value and $\mathrm{R}$OPE overlap that region, the decision by HDI would be not to reject, but the decision by ETI would be to reject. Again, I think the decision by HDI makes more sense, because these values have high credibility, even though they are in the extreme tail of the distribution.

Proponents of using the ETI point out that the ETI limits are invariant under nonlinear transformations of the parameter. The ETI limits of the transformed parameter are just the transformed limits of the original scale. This is not the case for HDI limits (in general). This property is handy when the parameters are arbitrarily scaled in abstract model derivations, or in some applied models for which parameters might be nonlinearly transformed for different purposes. But in most applications, the parameters are meaningfully defined on the canonical scale of the data, and the HDI has meaning relative to that scale. Nevertheless, it is important to recognize that if the scale of the parameter is nonlinearly transformed, the HDI limits will change relative to the percentiles of the distribution.

12.2.1 Is a coin fair or not?

For the null hypothesis, the prior distribution is a “spike” at the null value. The prior probability is zero for all values of θ other than the null value. The probability of the data for the null hypothesis is

$p\left(z,\text{?}N|M_{\text{null}}\right)={\theta }_{\text{null}}^z{\left(1-{\theta }_{\text{null}}\right)}^{\left(N-z\right)}$

where M_null denotes the null model, that is, the null hypothesis. For the alternative hypothesis, we assume a broad beta distribution. Recall from Footnote 5 on p. 132 that for a single coin with a beta prior, the marginal likelihood is

$p\text{?}\left(z,\text{?}N|M_{\text{alt}}\right)=B\text{?}\left(z+a_{\text{alt}},\text{?}N-z+b_{\text{alt}}\right)/B\text{?}\left(a_{\text{alt}},\text{?}{b}_{\text{alt}}\right)$

This equation was expressed as functions in $\mathrm{R}$ immediately after Equation 10.6 on p. 270. Combining Equations 12.2 and 12.1 we get the Bayes’ factor for the alternative hypothesis relative to the null hypothesis:

$\frac{p\text{?}\left(z,\text{?}N|M_{\text{alt}}\right)}{p\text{?}\left(z,\text{?}N|M_{\text{null}}\right)}=\frac{B\text{?}\left(z+a_{\text{alt}},\text{?}N-z+b_{\text{alt}}\right)/B\text{?}\left(a_{\text{alt}},\text{?}{b}_{\text{alt}}\right)}{{\theta }_{\text{null}}^z\text{?}{\left(1-{\theta }_{\text{null}}\right)}^{(N-z)}}$

For a default alternative prior, the beta distribution is supposed to be uninformed, according to particular mathematical criteria. Intuition might suggest that a uniform distribution suits this requirement, that is, beta(θ|1, 1)Instead, some argue that the most appropriate uninformed beta distribution is beta(θ | ϵ, ϵ),where ϵ is a small number approaching zero (e.g., Lee & Webb, 2005; Zhu & Lu, 2004). This is called the Haldane prior (as was mentioned in Section 10.6). A Haldane prior is illustrated in the top-left plot of Figure 12.3, using ϵ = 0.01.

Let's compute the value of the Bayes’ factor in Equation 12.3, for the data used repeatedly in the previous chapter on NHST, namely z = 7 and N = 24. Here are the results for various values of a_alt and b_alt:

$\frac{p\text{?}\left(z,\text{?}N|M_{\text{alt}}\right)}{p\text{?}\left(z,\text{?}N|M_{\text{null}}\right)}=\{\begin{array}{l}3.7227\text{for}{a}_{\text{alt}}=2,b_{\text{alt}}=4\hfill \\ 1.9390\text{for}{a}_{\text{alt}}=\text{?}{b}_{\text{alt}}=1.000\hfill \\ 0.4211\text{for}{a}_{\text{alt}}=\text{?}{b}_{\text{alt}}=0.100\hfill \\ 0.0481\text{for}{a}_{\text{alt}}=\text{?}{b}_{\text{alt}}=0.010\hfill \\ 0.0049\text{for}{a}_{\text{alt}}=\text{?}{b}_{\text{alt}}=0.001\hfill \end{array}$

The first case, a_alt = 2, b_alt = 4, will be discussed later. For now, notice that when the alternative prior is uniform, with a_alt = b_alt = 1.000, the Bayes’ factor shows a (small) preference for the alternative hypothesis, but when the alternative prior approximates the Haldane, the Bayes’ factor shows a strong preference for the null hypothesis. As the alternative prior gets closer to the Haldane limit, the Bayes’ factor changes by orders of magnitude. Thus, as we have seen before (e.g. Section 10.6, p. 292), the Bayes’ factor is very sensitive to the choice of prior distribution.

You can see why this happens by considering Figure 12.3, which shows two priors in its two columns. The Haldane prior, in the left column, puts extremely small prior credibility on the parameter values that are most consistent with the data. Because the marginal likelihood used in the Bayes’ factor is the product of the prior and the likelihood, the marginal likelihood from the Haldane prior will be relatively small. The lower left plot of Figure 12.3 indicates that the marginal likelihood is p(D) = 2.87 × 10^-9, which is small compared to the probability of the data from the null hypothesis, $p\left(D\right)={\theta }_{\text{null}}^z{\left(1-{\theta }_{\text{null}}\right)}^{\left(N-z\right)}=5.96\times {10}^{-8}$. The right column of Figure 12.3 uses a mildly informed prior (discussed below) that puts modestly high probability on the parameter values that are most consistent with the data. Because of this, the marginal likelihood is relatively high, at p(D) = 2.22 × 10^-7.

If we consider the posterior distribution instead of the Bayes’ factor, we see that the posterior distribution on θ within the alternative model is only slightly affected by the prior. With z = 7 and N = 24, for the uniform a_alt = b_alt = 1.00, the 95% HDI is [0.1407,0.4828]. For the approximate Haldane a_alt = b_alt = 0.01, the 95% HDI is [0.1222,0.4710], as shown in the lower-left plot of Figure 12.3. And for the mildly informed prior a_alt = 2, b_alt = 4, the 95% HDI is [0.1449,0.4624],as shown in the lower-right plot of Figure 12.3. (These HDI limits were accurately determined by function optimization using HDI of ICDF in DBDA2E-utilities. R, not by MCMC.) The lower and upper limits vary by only about 2 percentage points. In all cases, the 95% HDI excludes the null value, although a wide $\mathrm{R}$OPE might overlap the HDI. Thus, the explicit estimation of the bias parameter robustly indicates that the null value should be rejected, but perhaps only marginally. This contrasts with the Bayes’ factor, model-comparison approach, which rejected the null or accepted the null depending on the alternative prior.

Of the Bayes’ factors in Equation 12.4, which is most appropriate? If your analysis is driven by the urge for a default, uninformed alternative prior, then the prior that best approximates the Haldane is most appropriate. Following from that, we should strongly prefer the null hypothesis to the Haldane alternative. While this is mathematically correct, it is meaningless for an applied setting because the Haldane alternative represents nothing remotely resembling a credible alternative hypothesis. The Haldane prior sets prior probabilities of virtually zero at all values of θ except θ = 0 and θ = 1. There are very few applied settings where such a U-shaped prior represents a genuinely meaningful theory.

I recommended back in Section 10.6.1, p. 294, that the priors of the models should be equivalently informed. In the present application, the null-hypothesis prior is, by definition, fixed. But the alternative prior should be whatever best represents a meaningful and credible hypothesis, not a meaningless default. Suppose, for example, that we have some mild prior information that the coin is tail-biased. We express this as fictional prior data containing 1 head in 4 flips. With a uniform “proto-prior,” that implies the alternative prior should be beta(θ|1+1, 3+1) = beta(θ|2, 4), as shown in the upper-right plot of Figure 12.3. The Bayes’ factor for this meaningful alternative prior is given as the first case in Equation 12.4, where is can be seen that the null is rejected. This agrees with the conclusion from explicit estimation of the parameter value in the posterior distribution. Exercise 12.1 has you generate these examples for yourself. More extensive discussion, and an example from extrasensory perception, can be found in Kruschke (2011a).

12.2.1.1 Bayes' factor can accept null with poor precision

Figure 12.4 shows two examples in which the Bayes’ factor favors the null hypothesis. In these cases, the data have a proportion of 50% heads, which is exactly consistent with the null value θ = 0.5. The left column of Figure 12.4 uses a Haldane prior (with ϵ = 0.01), and the data comprise only a single head in two flips. The Bayes' factor is 51.0 in favor of the null hypothesis! But should we really believe therefore that θ = 0.5? No, I don't think so, because the posterior distribution on θ has a 95% HDI from 0.026 to 0.974.

The right column of Figure 12.4 uses a uniform prior. The data show 7 heads in 14 flips. The resulting Bayes' factor is 3.14 in favor of the null. But should we really believe therefore that θ = 0.5? No, I don't think so, because the posterior distribution on θ has a 95% HDI from 0.266 to 0.734.

Thus, the Bayes' factor can favor the null hypothesis when the data are consistent with it, even for small quantities of data. For another example, involving continuous instead of dichotomous data, see Appendix D of Kruschke (2013a). The problem is that there is very little precision in the estimate of θ for small quantities of data. It seems inappropriate, even contradictory, to accept a particular value of θ when there is so much uncertainty about it.

When using the estimation approach instead of the model-comparison approach, accepting the null value demands that the HDI falls entirely inside a $\mathrm{R}$OPE, which typically demands high precision. Narrower $\mathrm{R}$OPEs require higher precision to accept the null value. For example, if we use the narrow $\mathrm{R}$OPE from θ = 0.48 to θ = 0.52 as in Figure 12.4, we would need N = 2400 and z = 1200 for the 95% HDI to fall inside the $\mathrm{R}$OPE. There is no way around this inconvenient statistical reality: high precision demands a large sample size (and a measurement device with minimal possible noise). But when we are trying to accept a specific value of θ, is seems logically appropriate that we should have a reasonably precise estimate indicating that specific value.

12.2.2 Are different groups equal or not?

In many research applications, data are collected from subjects in different conditions, groups, or categories. We saw an example with baseball batting abilities in which players were grouped by fielding position (e.g., Figure 9.14, p. 256), but there are many other situations. Experiments often have different subjects in different treatment conditions. Observational studies often measure subjects from different classifications, such as gender, location, etc. Researchers often want to ask the question, Are the groups different or not?

As a concrete example, suppose we conduct an experiment about the effect of background music on the ability to remember. As a simple test of memory, each person tries to memorize the same list of 20 words (such as “chair,” “shark” “radio,” etc.). They see each word for a specific time, and then, after a brief retention interval, recall as many words as they can. For simplicity, we assume that all words are equally memorable, and a person's ability to recall words is modeled as a Bernoulli distribution, with probability θ_ij for the ith person in the jth condition. The individual recall propensities θ_ij depends on condition-level parameters, ω_j and κ_j, that describe the overarching recall propensity in each condition, because θ_ij ~ dbeta(ω_j (κ_j − 2) +1, (1 − ω_j)(κ_j − 2) +1).

The only difference between the conditions is the type of music being played during learning and recall. For the four groups, the music comes from, respectively, the death-metal band “Das Kruschke”³, Mozart, Bach, and Beethoven. For the four conditions, the mean number of words recalled is 8.0, 10.0, 10.2, and 10.4.

The most straight-forward way to find out whether the different types of music produce different memory abilities is to estimate the condition-level parameters and then examine the posterior differences of the parameter estimates. The histograms in the upper part of Figure 12.5 show the distributions of differences between the ω_j parameters. It can be seen that ω₁ is quite different than ω₃ and ω₄, and possibly ω₂. A difference of zero falls well outside the 95% HDI intervals, even taking into account MCMC instability and a small $\mathrm{R}$OPE. From this we would conclude that Das Kruschke produces poorer memory than the classical composers.

A model-comparison approach addresses the issue a different way. It compares the full model, which has distinct ω_j parameters for the four conditions, against a restricted model, which has a shared ω₀ parameter to describe all the conditions simultaneously. The two models have equal (50/50) prior probabilities. The bottom panel of Figure 12.5 shows the results of a model comparison. The single ω₀ model is preferred over the distinct ω_j model, by about 85% to 15%. In other words, from the model comparison we might conclude that there is no difference in memory between the groups.

Which analysis should we believe? Is condition 1 different from some other conditions, as the parameter estimate implies, or are all the conditions the same, as the model comparison seems to imply? Consider carefully what the model comparison actually says: Given the choice between one shared mode and four different group modes, the one-mode model is less improbable. But that does not mean that the one-mode model is the best possible model. In fact, if a different model comparison is conducted, that compares the one-mode model against a different model that has one mode for group 1 and a second mode that is shared for groups 2 through 4, then the comparison favors the two-mode model. Exercise 12.2 has you carry out this alternative comparison.

In principle, we could consider all possible models formed by partitioning the four groups. For four groups, there are 15 distinct partitions. We could, in principle, put a prior probability on each of the 15 models, and then do a comparison of the 15 models (Gopalan & Berry, 1998). From the posterior probabilities of the models, we could ascertain which partition was most credible, and decide whether it is more credible than other nearly-as-credible partitions. (Other approaches have been described by D. A. Berry & Hochberg, 1999; Mueller, Parmigiani, & Rice, 2007; Scott & Berger, 2006). Suppose that we conducted such a large-scale model comparison, and found that the most credible model partitioned groups 2-4 together, separate from group 1. Does this mean that we should truly believe that there is zero difference between groups 2, 3, and 4? Not necessarily. If the group treatments are different, such as the four types of music in the present scenario, then there is almost certainly at least some small difference between their outcomes. (In fact, the simulated data do come from groups with all different means.) We may still want to estimate the magnitude of those small differences, even if they are small. An explicit posterior estimate will reveal the magnitude and uncertainty of those estimates. Thus, unless we have a viable reason to believe that different group parameters may be literally identical, an estimation of distinct group parameters will tell us what we want to know, without model comparison.

model {
 for(sin 1:nSubj) {
 nCorrOfSubj[s] ~ dbin(theta[s], nTrlOfSubj[s])
 theta[s] ~ dbeta( aBeta[CondOfSubj[s]], bBeta[CondOfSubj[s]])
 }

for (j in 1:nCond) {
 # Use omega[j] for model index 1, omega0 for model index 2:
 aBeta[j] <-  (equals(mdlIdx,1)*omega[j]
 + equals(mdlIdx,2)*omega0 ) * (kappa[j]-2)+1
 bBeta[j] <- (1 - (equals(mdlIdx,1)*omega[j]
 + equals(mdlIdx,2)*omega0)) * (kappa[j]-2)+1
 omega[j] ~ dbeta(a[j,mdlIdx], b[j,mdlIdx])
}
omega0 ~ dbeta(a0[mdlIdx], b0[mdlIdx])

for(j in 1:nCond) {
 kappa[j] <- kappaMinusTwo[j] + 2
 kappaMinusTwo[j] ~ dgamma(2.618, 0.0809) # mode 20 . sd 20
}

 # Constants for prior and pseudoprior:
 aP <- 1
 bP <- 1
 # a0[model] and b0[model]
 a0[1] <- 0.48*500  # pseudo
 b0[1] <- (1-0.48)*500  # pseudo
 a0[2] <- aP  # true
 b0[2] <- bP  # true
 # a[condition, model] and b[condition, model]
 a[1,1] <- aP  # true
 a[2,1] <- aP  # true
 a[3,1] <- aP  # true
 a[4,1] <- aP  # true
 b[1,1] <- bP  # true
 b[2,1] <- bP  # true
 b[3,1] <- bP  # true
 b[4,1] <- bP  # true
 a[1,2] <- 0.40*125  # pseudo
 a[2,2] <- 0.50*125  # pseudo
 a[3,2] <- 0.51*125  # pseudo
 a[4,2] <- 0.52*125  # pseudo
 b[1,2] <- (1-0.40)*125  # pseudo
 b[2,2] <- (1-0.50)*125  # pseudo
 b[3,2] <- (1-0.51)*125  # pseudo
 b[4,2] <- (1-0.52)*125  # pseudo
 # Prior on model index:
 mdlIdx ~ dcat(modelProb[])
 modelProb[1] <- 0.5
 modelProb[2] <- 0.5
}