Bayesian Testing of Constrained Hypotheses

Abstract

Statistical hypothesis testing plays a central role in applied research to determine whether theories or expectations are supported by the data or not. Such expectations are often formulated using order constraints. For example an executive board may expect that sales representatives who wear a smart watch will respond faster to their emails than sales representatives who don’t wear a smart watch. In addition it may be expected that this difference becomes more pronounced over time because representatives need to learn how to use the smart watch effectively. By translating these expectations into statistical hypotheses with equality and/or order constraints we can determine whether the expectations receive evidence from the data. In this chapter we show how a Bayesian statistical approach can effectively be used for this purpose. This Bayesian approach is more flexible than the traditional p-value test in the sense that multiple hypotheses with equality as well as order constraints can be tested against each other in a direct manner. The methodology can straightforwardly be used by practitioners using the freely downloadable software package BIEMS. An application in a human-computer interaction is used for illustration.

Appendices

Appendix 1: Gibbs sampler (theory)

We consider the general case of P repeated measurements. In the example discussed above, P was equal to 3. The following semi-conjugate prior is used for the model parameters,

$$\begin{aligned} p(\varvec{\mu }_A,\varvec{\mu }_B,\varvec{\varSigma })= & {} p(\varvec{\mu }_A)\times p(\varvec{\mu }_B)\times p(\varvec{\varSigma })\\ \nonumber\propto & {} N_{\varvec{\mu }_A}(\mathbf m _{A0},\mathbf S _{A0}) \times N_{\varvec{\mu }_A}(\mathbf m _{A0},\mathbf S _{A0})\times |\varvec{\varSigma }|^{-\frac{P+1}{2}}, \end{aligned}$$

(9.17)

where $\mathbf m _{A0}$ and $\mathbf m _{B0}$ are the prior mean of $\varvec{\mu }_A$ and $\varvec{\mu }_B$, respectively, and $\mathbf S _{A0}$ and $\mathbf S _{B0}$ the respective prior covariance matrices of $\varvec{\mu }_A$ and $\varvec{\mu }_B$.

The data are stored in the $(n_A+n_B)\times P$ data matrix $\mathbf Y =[\mathbf Y _A'~\mathbf Y _B']'$, where the ith row of $\mathbf Y$ contains the P measurements of ith sales executive and the first $n_A$ rows correspond to the responses of the executives in team A and the remaining $n_B$ rows contain the responses of the executives of team B. The likelihood of the data can be written as

$$\begin{aligned} p(\mathbf Y |\varvec{\mu }_A,\varvec{\mu }_B,\varvec{\varSigma })= & {} \prod _{i=1}^{n_A} p(\mathbf y _i|\varvec{\mu }_A,\varvec{\varSigma }) \times \prod _{i=n_A+1}^{n_A+n_B} p(\mathbf y _i|\varvec{\mu }_B,\varvec{\varSigma })\\\propto & {} N_{\varvec{\mu }_A|\varvec{\varSigma }}(\bar{\mathbf{y }}_A,\varvec{\varSigma }/n_A)\times N_{\varvec{\mu }_A|\varvec{\varSigma }}(\bar{\mathbf{y }}_B,\varvec{\varSigma }/n_B)\times IW_{\varvec{\varSigma }}(\mathbf S ,n_A+n_B-2), \end{aligned}$$

where $\bar{\mathbf{y }}_A$ and $\bar{\mathbf{y }}_B$ denote the sample means of team A and team B over the P measurements, the sums of squares equal

$$ \mathbf S =(\mathbf Y _A-\mathbf 1 _{n_A}\bar{\mathbf{y }}_A')'(\mathbf Y _A-\mathbf 1 _{n_A}\bar{\mathbf{y }}_A')+(\mathbf Y _B-\mathbf 1 _{n_B}\bar{\mathbf{y }}_B')'(\mathbf Y _B-\mathbf 1 _{n_B}\bar{\mathbf{y }}_B'), $$

and $IW_{\varvec{\varSigma }}(\mathbf S ,n)$ denotes an inverse Wishart probability density for $\varvec{\varSigma }$. Note that the likelihood function of $\varvec{\varSigma }$ given $\varvec{\mu }_A$ and $\varvec{\mu }_B$ is proportional to an inverse Wishart density $IW(\mathbf S _{\varvec{\mu }},n_A+n_B)$, where $\mathbf S _{\varvec{\mu }}=(\mathbf Y _A-\mathbf 1 _{n_A}\varvec{\mu }_A')'(\mathbf Y _A-\mathbf 1 _{n_A}\varvec{\mu }_A')+(\mathbf Y _B-\mathbf 1 _{n_B}\varvec{\mu }_B')'(\mathbf Y _B-\mathbf 1 _{n_B}\varvec{\mu }_B')$. These results can be found in most classic Bayesian text books, such as Gelman et al. (2004), for example.

Because the prior in (9.17) is semi-conjugate, the conditional posterior distributions of each model paramater given the other parameters have known distributions from which we can easily sample,

$$\begin{aligned} p(\varvec{\mu }_A|\mathbf Y ,\varvec{\varSigma })= & {} N\left( \left( \mathbf S _{A0}^{-1}+n_A\varvec{\varSigma }^{-1}\right) ^{-1}\left( \mathbf S _{A0}^{-1}{} \mathbf m _{A0}+n_A\varvec{\varSigma }^{-1}\bar{\mathbf{y }}_A \right) ,\left( \mathbf S _{A0}^{-1}+n_A\varvec{\varSigma }^{-1}\right) ^{-1}\right) \\ p(\varvec{\mu }_B|\mathbf Y ,\varvec{\varSigma })= & {} N\left( \left( \mathbf S _{B0}^{-1}+n_B\varvec{\varSigma }^{-1}\right) ^{-1}\left( \mathbf S _{B0}^{-1}{} \mathbf m _{B0}+n_A\varvec{\varSigma }^{-1}\bar{\mathbf{y }}_B \right) ,\left( \mathbf S _{B0}^{-1}+n_B\varvec{\varSigma }^{-1}\right) ^{-1}\right) \\ p(\varvec{\varSigma }|\mathbf Y ,\varvec{\mu }_A,\varvec{\mu }_B)= & {} IW(\mathbf S _{\varvec{\mu }},n_A+n_B). \end{aligned}$$

We can use a Gibbs sampler to get a sample from the joint posterior of $(\varvec{\mu }_A,\varvec{\mu }_B,\varvec{\varSigma })$. In a Gibbs sampler we sequentially draw each model parameter from its conditional posterior given the remaining parameters. The Gibbs sampler algorithm can be written as

• 1. Set initial values for the model parameters: $\varvec{\mu }_A^{(0)}$, $\varvec{\mu }_B^{(0)}$, and $\varvec{\varSigma }^{(0)}$.

• 2. Draw $\varvec{\mu }_A^{(s)}$ from its conditional posterior $p(\varvec{\mu }_A|\mathbf Y ,\varvec{\varSigma }^{(s-1)})$.

• 3. Draw $\varvec{\mu }_B^{(s)}$ from its conditional posterior $p(\varvec{\mu }_B|\mathbf Y ,\varvec{\varSigma }^{(s-1)})$.

• 4. Draw $\varvec{\varSigma }^{(s)}$ from its conditional posterior $p(\varvec{\varSigma }|\mathbf Y ,\varvec{\mu }_A^{(s)},\varvec{\mu }_A^{(s)})$.

• 5. Repeat steps 2–4 for $s=1,\ldots ,S$.

In the software program R, drawing from a multivariate normal distribution can be done using the function ‘rmvnorm’ in the ‘mvtnorm’-package and drawing from an inverse Wishart distribution can be done using the function ‘riwish’ in the ‘MCMCpack’-package.

It may be that the initial values, $\varvec{\mu }_A^{(0)}$, $\varvec{\mu }_B^{(0)}$, and $\varvec{\varSigma }^{(0)}$, are chosen far away from the subspace where the posterior is concentrated. If this is the case, a burn-in period of, say, 100 draws is needed. After the burn-in period convergence is reached and the remaining draws come from the actual posterior of the model parameters.

Appendix 2: Gibbs sampler (R code)

Conditional posteriors for $\varvec{\mu }_A$, $\varvec{\mu }_B$, and $\varvec{\varSigma }$.

Gibbs sampler

Generate data matrix Y

Compute classical estimates

Set priors for Gibbs sampler

Run Gibbs sampler

Compute descriptive statistics from Gibbs output

Create data matrix for BIEMS

Appendix 3: Derivation of the Bayes factor

The Bayes factor is derived for a one-sided hypothesis $H_1:\delta <0$ versus the unconstrained hypothesis $H_u:\delta \in \mathbb {R}$. In the encompassing prior approach, the prior under $H_1$, $p_1(\delta ,\sigma ^2)$, is a truncation of the unconstrained (or encompassing) prior under $H_u$, $p_u(\delta ,\sigma ^2)$, in the region where $\delta <0$, i.e., $p_1(\delta ,\sigma ^2)=p_u(\delta ,\sigma ^2)I(\delta <0)/\text{ Pr }(\delta <0|H_u)$, where the prior probability $\text{ Pr }(\delta <0|H_u)=\int _{\delta <0}p_u(\delta )d\delta$, where $I(\cdot )$ is the indicator function. Note that $\text{ Pr }(\delta <0|H_u)=\frac{1}{2}$ if the unconstrained prior is centered at 0, such as $p_u(\delta )=N(0,\sigma _0^2)$. Also note that the likelihood under $H_1$ is a truncation of the likelihood under $H_u$, i.e., $p_1(\mathbf y |\delta ,\sigma ^2)=p_u(\mathbf y |\delta ,\sigma ^2)I(\delta <0)$. For this reason we can omit the hypothesis index u in the likelihood functions in the derivation below. The Bayes factor of $H_1$ versus $H_u$ can then be derived as follows

$$\begin{aligned} \nonumber B_{1u}= & {} \frac{\iint _{\delta <0}p(\mathbf y |\delta ,\sigma ^2)p_1(\delta ,\sigma ^2) d \delta d\sigma ^2}{\iint p(\mathbf y |\delta ,\sigma ^2)p_u(\delta ,\sigma ^2)d \delta d \sigma ^2}\\ \nonumber= & {} \frac{1}{\text{ Pr }(\delta <0|H_u)}\iint _{\delta <0} \frac{p(\mathbf y |\delta ,\sigma ^2)p_u(\delta ,\sigma ^2)}{\iint p(\mathbf y |\delta ,\sigma ^2)p_u(\delta ,\sigma ^2)d \delta d \sigma ^2}d \delta d \sigma ^2 \\ \nonumber= & {} \frac{1}{\text{ Pr }(\delta <0|H_u)}\iint _{\delta <0} p_u(\delta ,\sigma ^2|\mathbf y )d \delta d \sigma ^2\\= & {} \frac{\text{ Pr }(\delta <0|\mathbf y ,H_u)}{\text{ Pr }(\delta <0|H_u)}, \end{aligned}$$

which corresponds to (9.15), where $\text{ Pr }(\delta <0|\mathbf y ,H_u)=\iint _{\delta <0}p_{u}(\delta |\mathbf y ) d\delta d\sigma ^2$ is the posterior probability that the constraints hold under $H_u$. For $H_2:\delta >0$ versus the unconstrained hypothesis $H_u:\delta \in \mathbb {R}$ we can follow the same steps to obtain (9.16).

For $H_0:\delta =0$ versus the unconstrained hypothesis $H_u:\delta \in \mathbb {R}$, the encompassing prior approach implies that $p_0(\sigma ^2)=p_u(\sigma ^2|\delta =0)$. Consequently,

$$\begin{aligned} \nonumber B_{0u}= & {} \frac{\int p(\mathbf y |\delta =0,\sigma ^2)p_0(\sigma ^2) d\sigma ^2}{\iint p(\mathbf y |\delta ,\sigma ^2)p_u(\delta ,\sigma ^2)d\delta d \sigma ^2}\\= & {} \frac{\int p(\mathbf y |\delta =0,\sigma ^2)p_u(\sigma ^2|\delta =0) d\sigma ^2}{\iint p(\mathbf y |\delta ,\sigma ^2)p_u(\delta ,\sigma ^2)d \delta d \sigma ^2} \\= & {} \frac{1}{p_u(\delta =0)}\int \frac{p(\mathbf y |\delta =0,\sigma ^2)p_u(\delta =0,\sigma ^2)}{\iint p(\mathbf y |\delta ,\sigma ^2)p_u(\delta ,\sigma ^2)d \delta d \sigma ^2} d\sigma ^2\\= & {} \frac{1}{p_u(\delta =0)}\int p_u(\delta =0,\sigma ^2|\mathbf y ) d\sigma ^2\\= & {} \frac{p_u(\delta =0|\mathbf y )}{p_u(\delta =0)}, \end{aligned}$$

which is equal to the Savage-Dickey density ratio in (9.14).