<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/i-dont-know-fishers-exact-test-but-i-know-stan/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://tjmahr.github.io/bayesian-fisher-exact-test/"> Higher Order Functions</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>A few days ago, I watched a <a target="_blank" rel="nofollow" href="https://www.youtube.com/watch?v=qQFF4tPgeWI">terrific lecture by Bob Carpenter</a> on<br>
Bayesian models. He started with a Bayesian approach to Fisher’s exact test. I<br>
had never heard of this classical procedure, so I was curious to play with the<br>
example. In this post, I use the same data that he used in the lecture and<br>
in an earlier, <a target="_blank" rel="nofollow" href="https://lingpipe-blog.com/2009/10/13/bayesian-counterpart-to-fisher-exact-test-on-contingency-tables/">pre-Stan blog post</a>. I show how I would go about<br>
fitting the model in Stan and inspecting the results in R.</p>
<h2>Problem statement</h2>
<p>We observed the following data.</p>
<table>
<thead>
<tr>
<th>sex</th>
<th>n left handed</th>
<th>n right handed</th>
</tr>
</thead>
<tbody>
<tr>
<td>male</td>
<td>9</td>
<td>43</td>
</tr>
<tr>
<td>female</td>
<td>4</td>
<td>44</td>
</tr>
</tbody>
</table>
<p>Question: Is the rate of left-handedness different between the male and female<br>
groups? Specifically, is left-handedness more likely in the male group?</p>
<h3>Classical approach</h3>
<p>In frequentist statistics, we might run Fisher’s exact test. At least, that’s what<br>
the <a target="_blank" rel="nofollow" href="https://www.google.com/search?q=which+statistical+test+to+use&tbm=isch">flow charts</a> tell us. To do that, we first put the data in a<br>
matrix.</p>
<div>
<pre><code><span># Create a matrix representation of the data for the fisher test
</span><span>m</span><span> </span><span>&lt;-</span><span> </span><span>matrix</span><span>(</span><span>c</span><span>(</span><span>9</span><span>,</span><span> </span><span>43</span><span>,</span><span> </span><span>4</span><span>,</span><span> </span><span>44</span><span>),</span><span> </span><span>nrow</span><span> </span><span>=</span><span> </span><span>2</span><span>,</span><span> </span><span>byrow</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>,</span><span>
</span><span>dimnames</span><span> </span><span>=</span><span> </span><span>list</span><span>(</span><span>sex</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>&quot;male&quot;</span><span>,</span><span> </span><span>&quot;female&quot;</span><span>),</span><span>
</span><span>handedness</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>&quot;left&quot;</span><span>,</span><span> </span><span>&quot;right&quot;</span><span>)))</span><span>
</span><span>m</span><span>
</span><span>#> handedness
#> sex left right
#> male 9 43
#> female 4 44
</span></code></pre>
</div>
<p>We can run the two-tailed test: Are the two groups different?</p>
<div>
<pre><code><span>fisher.test</span><span>(</span><span>m</span><span>)</span><span>
</span><span>#> #> Fisher's Exact Test for Count Data
#> #> data: m
#> p-value = 0.2392
#> alternative hypothesis: true odds ratio is not equal to 1
#> 95 percent confidence interval:
#> 0.582996 10.927993
#> sample estimates:
#> odds ratio #> 2.283832
</span></code></pre>
</div>
<p>The output is a little verbose, but I do like how it spells out a sentence<br>
describing the alternative hypothesis.</p>
<p>We can also consider one-sided test: Is left-handedness greater in the male<br>
group?</p>
<div>
<pre><code><span>fisher.test</span><span>(</span><span>m</span><span>,</span><span> </span><span>alternative</span><span> </span><span>=</span><span> </span><span>&quot;greater&quot;</span><span>)</span><span>
</span><span>#> #> Fisher's Exact Test for Count Data
#> #> data: m
#> p-value = 0.1502
#> alternative hypothesis: true odds ratio is greater than 1
#> 95 percent confidence interval:
#> 0.7006563 Inf
#> sample estimates:
#> odds ratio #> 2.283832
</span></code></pre>
</div>
<p>In both tests, we cannot reject the null hypothesis because the <em>p</em>-value is<br>
greater than .05 :x:, so we would conclude that the two groups are not different.</p>
<p>But I don’t really know this test that well. We never covered it in any of my<br>
stats classes, and indeed, this post is the first time I ever used the function<br>
<code>fisher.test()</code>. If I had never heard of the test, I am not quite sure what I<br>
would have done. Maybe a logistic regression (<em>p</em> = .191 :x:).</p>
<h2>Creating a Stan model</h2>
<p>I’m very not fluent in the classical bag of tricks, but that’s okay. I know some<br>
Stan, I have an idea about how the data could have been generated, and that’s<br>
good enough :relaxed:. I can just write down my data-generating story in a model<br>
and let Stan compute a posterior distribution for the difference in handedness<br>
rates between the two groups.</p>
<p>For my model, I’m going to suppose that in each group, there is a probability of<br>
being left-handed called <em>θ</em> and that the counts we see result from a<br>
binomial process. The 9 left-handed males we observe are the number of successes<br>
from 52 observations of a process that “succeeds” with probability<br>
<em>θ</em><sub>male</sub>.</p>
<p>To fit the model, we need a prior distribution. The prior’s job is to generate<br>
possible <em>θ</em>’s, and we will use our data to update the prior. The Beta<br>
distribution generates values between 0 and 1, so it’s an obvious choice.</p>
<p>In math, the model for a group would be:</p>
<p></p>
<p>We can use a flat, uninformative prior by using <em>a</em> = 1, <em>b</em> = 1. This<br>
prior considers all probabilities from 0 to 1 as equally plausible.</p>
<div>
<pre><code><span>library</span><span>(</span><span>ggplot2</span><span>)</span><span>
</span><span>steps</span><span> </span><span>&lt;-</span><span> </span><span>seq</span><span>(</span><span>from</span><span> </span><span>=</span><span> </span><span>0</span><span>,</span><span> </span><span>to</span><span> </span><span>=</span><span> </span><span>1</span><span>,</span><span> </span><span>by</span><span> </span><span>=</span><span> </span><span>.01</span><span>)</span><span>
</span><span>ggplot</span><span>(</span><span>data.frame</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>steps</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>dbeta</span><span>(</span><span>steps</span><span>,</span><span> </span><span>shape1</span><span> </span><span>=</span><span> </span><span>1</span><span>,</span><span> </span><span>shape2</span><span> </span><span>=</span><span> </span><span>1</span><span>)))</span><span> </span><span>+</span><span> </span><span>geom_line</span><span>(</span><span>aes</span><span>(</span><span>x</span><span>,</span><span> </span><span>y</span><span>))</span><span> </span><span>+</span><span>
</span><span>scale_x_continuous</span><span>(</span><span>breaks</span><span> </span><span>=</span><span> </span><span>(</span><span>0</span><span>:</span><span>10</span><span>)</span><span> </span><span>/</span><span> </span><span>10</span><span>)</span><span> </span><span>+</span><span> </span><span>ylim</span><span>(</span><span>0</span><span>,</span><span> </span><span>2</span><span>)</span><span> </span><span>+</span><span>
</span><span>labs</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>&quot;p(left-handed)&quot;</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>&quot;density&quot;</span><span>,</span><span> </span><span>title</span><span> </span><span>=</span><span> </span><span>&quot;beta(1,1)&quot;</span><span>)</span><span>
</span></code></pre>
</div>
<p><img title="Density plot of the flat, uninformative prior." alt="Density plot of the flat, uninformative prior." src="https://i2.wp.com/tjmahr.github.io/figs//2017-05-16-bayesian-fisher-exact-test/flat-1.png?w=80%25&ssl=1" data-recalc-dims="1"></p>
<p>But I also think that 10-ish% of people are left handed. (I don’t know where I<br>
first heard this number, but it’ll serve as my prior information.) I toyed<br>
around with <code>shape1</code> and <code>shape2</code> parameters in <code>dbeta()</code> until I got the prior<br>
Beta(5, 40), which is peaked around .1-ish but wide enough to keep .5 and .15 as<br>
plausible values too.</p>
<div>
<pre><code><span>ggplot</span><span>(</span><span>data.frame</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>steps</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>dbeta</span><span>(</span><span>steps</span><span>,</span><span> </span><span>shape1</span><span> </span><span>=</span><span> </span><span>5</span><span>,</span><span> </span><span>shape2</span><span> </span><span>=</span><span> </span><span>40</span><span>)))</span><span> </span><span>+</span><span> </span><span>geom_line</span><span>(</span><span>aes</span><span>(</span><span>x</span><span>,</span><span> </span><span>y</span><span>))</span><span> </span><span>+</span><span>
</span><span>scale_x_continuous</span><span>(</span><span>breaks</span><span> </span><span>=</span><span> </span><span>(</span><span>0</span><span>:</span><span>10</span><span>)</span><span> </span><span>/</span><span> </span><span>10</span><span>)</span><span> </span><span>+</span><span>
</span><span>labs</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>&quot;p(left-handed)&quot;</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>&quot;density&quot;</span><span>,</span><span> </span><span>title</span><span> </span><span>=</span><span> </span><span>&quot;beta(5,40)&quot;</span><span>)</span><span>
</span></code></pre>
</div>
<p><img title="Density plot of the informative prior." alt="Density plot of the informative prior." src="https://i1.wp.com/tjmahr.github.io/figs//2017-05-16-bayesian-fisher-exact-test/informative-1.png?w=80%25&ssl=1" data-recalc-dims="1"></p>
<p>Let’s write out a really simple model in Stan. Okay, it <em>used to</em> be really<br>
simple. Then I made the parameters for the Beta prior data values, and then I<br>
created an option to just sample the prior distribution. But the core of it is<br>
simple. The most important lines are the one with <code>~</code> symbols. These<br>
correspond to the sampling statements in the mathematical description of the<br>
model.</p>
<div>
<pre><code><span>model_code</span><span> </span><span>&lt;-</span><span> </span><span>&quot;
data {
int&lt;lower=0> beta_a;
int&lt;lower=0> beta_b;
int&lt;lower=0> n_total_1;
int&lt;lower=0> n_total_2;
int&lt;lower=0> n_hits_1;
int&lt;lower=0> n_hits_2;
int&lt;lower=0, upper=1> sample_prior_only;
}
parameters { real&lt;lower=0, upper=1> theta_1;
real&lt;lower=0, upper=1> theta_2;
}
model {
theta_1 ~ beta(beta_a, beta_b);
theta_2 ~ beta(beta_a, beta_b);
if (sample_prior_only != 1) {
n_hits_1 ~ binomial(n_total_1, theta_1);
n_hits_2 ~ binomial(n_total_2, theta_2);
}
}
generated quantities {
real diff;
diff = theta_1 - theta_2;
}
&quot;</span><span>
</span></code></pre>
</div>
<p>The <code>generated quantities</code> block runs on every sample of the posterior<br>
distribution. Here, we compute the<br>
<em>θ</em><sub>male</sub> − <em>θ</em><sub>female</sub> on every<br>
draw. Computing the difference inside the model code means that Stan<br>
will treat the <code>diff</code> values like any other parameter of the model. It will show<br>
up in summary functions and in plots of model parameters. That saves us some<br>
work later on.</p>
<div>
<pre><code><span>library</span><span>(</span><span>rstan</span><span>)</span><span>
</span><span>#> Loading required package: StanHeaders
#> rstan (Version 2.15.1, packaged: 2017-04-19 05:03:57 UTC, GitRev: 2e1f913d3ca3)
#> For execution on a local, multicore CPU with excess RAM we recommend calling
#> rstan_options(auto_write = TRUE)
#> options(mc.cores = parallel::detectCores())
</span></code></pre>
</div>
<p>I begin by compiling the model. This step will create an executable<br>
program that can sample from the model. I do the compilation in its own step so<br>
that I can re-use the program for different versions of the model.</p>
<div>
<pre><code><span>model_program</span><span> </span><span>&lt;-</span><span> </span><span>stan_model</span><span>(</span><span>model_code</span><span> </span><span>=</span><span> </span><span>model_code</span><span>)</span><span>
</span></code></pre>
</div>
<p>For convenience, I wrote a function that fits different versions of this model.<br>
This step is not necessary, but I don’t like repeating myself. (Normally, you<br>
would use <code>sampling(model_program, stan_data)</code> to get samples from a<br>
model-program.)</p>
<div>
<pre><code><span>run_model</span><span> </span><span>&lt;-</span><span> </span><span>function</span><span>(</span><span>beta_a</span><span>,</span><span> </span><span>beta_b</span><span>,</span><span> </span><span>sample_prior_only</span><span>)</span><span> </span><span>{</span><span>
</span><span>stan_data</span><span> </span><span>&lt;-</span><span> </span><span>list</span><span>(</span><span>
</span><span>beta_a</span><span> </span><span>=</span><span> </span><span>beta_a</span><span>,</span><span>
</span><span>beta_b</span><span> </span><span>=</span><span> </span><span>beta_b</span><span>,</span><span>
</span><span>n_total_1</span><span> </span><span>=</span><span> </span><span>9</span><span> </span><span>+</span><span> </span><span>43</span><span>,</span><span>
</span><span>n_total_2</span><span> </span><span>=</span><span> </span><span>4</span><span> </span><span>+</span><span> </span><span>44</span><span>,</span><span>
</span><span>n_hits_1</span><span> </span><span>=</span><span> </span><span>9</span><span>,</span><span>
</span><span>n_hits_2</span><span> </span><span>=</span><span> </span><span>4</span><span>,</span><span>
</span><span>sample_prior_only</span><span> </span><span>=</span><span> </span><span>sample_prior_only</span><span>)</span><span>
</span><span># Use quietly() to hide the sampler's output text
</span><span> </span><span>model</span><span> </span><span>&lt;-</span><span> </span><span>purrr</span><span>::</span><span>quietly</span><span>(</span><span>sampling</span><span>)(</span><span>model_program</span><span>,</span><span> </span><span>stan_data</span><span>)</span><span>
</span><span># But print any warnings that would have appeared
</span><span> </span><span>invisible</span><span>(</span><span>lapply</span><span>(</span><span>model</span><span>$</span><span>warnings</span><span>,</span><span> </span><span>warning</span><span>,</span><span> </span><span>call.</span><span> </span><span>=</span><span> </span><span>FALSE</span><span>))</span><span>
</span><span>model</span><span>$</span><span>result</span><span>
</span><span>}</span><span>
</span></code></pre>
</div>
<p>The actual left-handed versus right-handed numbers are hard-coded, but I can<br>
adjust the parameters for the Beta prior and toggle between sampling from the<br>
prior and the posterior. Stan normally prints out verbose progress information<br>
but I suppress that by using <code>purrr::quietly()</code>. I still want warnings, so I<br>
print them if they arise.</p>
<h3>Checking our prior information</h3>
<p>Now, let’s draw samples from the priors of each model. This step lets us check<br>
that our program works as expected. We know how the values should be<br>
distributed—we made up the numbers!</p>
<div>
<pre><code><span>m_informative_pd</span><span> </span><span>&lt;-</span><span> </span><span>run_model</span><span>(</span><span>beta_a</span><span> </span><span>=</span><span> </span><span>5</span><span>,</span><span> </span><span>beta_b</span><span> </span><span>=</span><span> </span><span>40</span><span>,</span><span> </span><span>sample_prior_only</span><span> </span><span>=</span><span> </span><span>1</span><span>)</span><span>
</span><span>m_informative_pd</span><span>
</span><span>#> Inference for Stan model: 7e241131e72a0ec0b2cfd4f9a73290f1.
#> 4 chains, each with iter=2000; warmup=1000; thin=1; #> post-warmup draws per chain=1000, total post-warmup draws=4000.
#> #> mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
#> theta_1 0.11 0.00 0.05 0.04 0.08 0.11 0.14 0.22 3047 1
#> theta_2 0.11 0.00 0.05 0.04 0.08 0.11 0.14 0.21 2953 1
#> diff 0.00 0.00 0.06 -0.13 -0.04 0.00 0.04 0.13 3072 1
#> lp__ -32.40 0.02 0.99 -35.01 -32.81 -32.11 -31.68 -31.42 1776 1
#> #> Samples were drawn using NUTS(diag_e) at Tue May 16 14:05:44 2017.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split chains (at #> convergence, Rhat=1).
</span><span>
</span><span>m_flat_pd</span><span> </span><span>&lt;-</span><span> </span><span>run_model</span><span>(</span><span>beta_a</span><span> </span><span>=</span><span> </span><span>1</span><span>,</span><span> </span><span>beta_b</span><span> </span><span>=</span><span> </span><span>1</span><span>,</span><span> </span><span>sample_prior_only</span><span> </span><span>=</span><span> </span><span>1</span><span>)</span><span>
</span><span>m_flat_pd</span><span>
</span><span>#> Inference for Stan model: 7e241131e72a0ec0b2cfd4f9a73290f1.
#> 4 chains, each with iter=2000; warmup=1000; thin=1; #> post-warmup draws per chain=1000, total post-warmup draws=4000.
#> #> mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
#> theta_1 0.49 0.01 0.28 0.03 0.25 0.49 0.74 0.97 3134 1
#> theta_2 0.50 0.00 0.28 0.03 0.26 0.51 0.74 0.97 3247 1
#> diff -0.01 0.01 0.40 -0.77 -0.29 -0.01 0.28 0.75 3001 1
#> lp__ -3.90 0.03 1.09 -6.86 -4.34 -3.56 -3.12 -2.81 1362 1
#> #> Samples were drawn using NUTS(diag_e) at Tue May 16 14:05:45 2017.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split chains (at #> convergence, Rhat=1).
</span></code></pre>
</div>
<p>We can confirm (by inspecting the <code>diff</code> row) that the difference in<br>
left-handedness in both groups is 0 according to our priors. The<br>
informative prior says that the values between 0.04 and<br>
0.22 are plausible rates of left-handedness in each group.</p>
<p>It’s worth a moment to reflect on how obviously wrong the uninformative prior is.<br>
The central <em>θ</em> value in each group is .5. Therefore, half of the prior<br>
samples assert that there are more left-handed individuals than right-handed<br>
ones! If there is anything we know about handedness, it’s that left-handedness is<br>
less common than right-handedness. <em>Uninformative</em> sometimes connotes “unbiased”<br>
or “letting the data speak for itself”, but in this case, I would say<br>
“gullible”.</p>
<h2>Sampling the posterior</h2>
<p>Now, the fun part. We update our prior information with our data.</p>
<div>
<pre><code><span>m_informative</span><span> </span><span>&lt;-</span><span> </span><span>run_model</span><span>(</span><span>beta_a</span><span> </span><span>=</span><span> </span><span>5</span><span>,</span><span> </span><span>beta_b</span><span> </span><span>=</span><span> </span><span>40</span><span>,</span><span> </span><span>sample_prior_only</span><span> </span><span>=</span><span> </span><span>0</span><span>)</span><span>
</span><span>m_informative</span><span>
</span><span>#> Inference for Stan model: 7e241131e72a0ec0b2cfd4f9a73290f1.
#> 4 chains, each with iter=2000; warmup=1000; thin=1; #> post-warmup draws per chain=1000, total post-warmup draws=4000.
#> #> mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
#> theta_1 0.14 0.00 0.04 0.08 0.12 0.14 0.17 0.22 3288 1
#> theta_2 0.10 0.00 0.03 0.05 0.07 0.09 0.11 0.17 2837 1
#> diff 0.05 0.00 0.05 -0.05 0.02 0.05 0.08 0.14 3020 1
#> lp__ -70.65 0.02 1.05 -73.56 -71.09 -70.33 -69.89 -69.63 1834 1
#> #> Samples were drawn using NUTS(diag_e) at Tue May 16 14:05:45 2017.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split chains (at #> convergence, Rhat=1).
</span><span>
</span><span>m_flat</span><span> </span><span>&lt;-</span><span> </span><span>run_model</span><span>(</span><span>beta_a</span><span> </span><span>=</span><span> </span><span>1</span><span>,</span><span> </span><span>beta_b</span><span> </span><span>=</span><span> </span><span>1</span><span>,</span><span> </span><span>sample_prior_only</span><span> </span><span>=</span><span> </span><span>0</span><span>)</span><span>
</span><span>m_flat</span><span>
</span><span>#> Inference for Stan model: 7e241131e72a0ec0b2cfd4f9a73290f1.
#> 4 chains, each with iter=2000; warmup=1000; thin=1; #> post-warmup draws per chain=1000, total post-warmup draws=4000.
#> #> mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
#> theta_1 0.19 0.00 0.05 0.09 0.15 0.18 0.22 0.30 2987 1
#> theta_2 0.10 0.00 0.04 0.04 0.07 0.10 0.13 0.20 2511 1
#> diff 0.09 0.00 0.07 -0.04 0.04 0.09 0.13 0.21 3042 1
#> lp__ -43.14 0.02 1.00 -45.82 -43.54 -42.84 -42.43 -42.15 1768 1
#> #> Samples were drawn using NUTS(diag_e) at Tue May 16 14:05:46 2017.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split chains (at #> convergence, Rhat=1).
</span></code></pre>
</div>
<p>The flat model puts the difference at 0.09 and 90% of the plausible<br>
values fall in the interval [-0.02, 0.19]. The informative<br>
model is more skeptical of higher left-handedness rates, so it puts the difference<br>
at 0.05 with 90% of the values between [-0.03, 0.13].<br>
Both of these intervals contain 0 and negative values :x:, so there is not much<br>
evidence for higher left-handedness in the male group.</p>
<p>To compute a “Bayesian <em>p</em>-value”, we could ask what proportion of differences<br>
are 0 or negative. There are more proper ways to make this inference in a<br>
Bayesian framework, but this approach is the easiest and it works for a model<br>
this simple. If 10% of the plausible values for the group differences are<br>
negative, then we assign a 10% probability to a negative group difference.</p>
<div>
<pre><code><span>df_flat</span><span> </span><span>&lt;-</span><span> </span><span>as.data.frame</span><span>(</span><span>m_flat</span><span>)</span><span>
</span><span>mean</span><span>(</span><span>df_flat</span><span>$</span><span>diff</span><span> </span><span>&lt;=</span><span> </span><span>0</span><span>)</span><span>
</span><span>#> [1] 0.09575
</span><span>
</span><span>df_informative</span><span> </span><span>&lt;-</span><span> </span><span>as.data.frame</span><span>(</span><span>m_informative</span><span>)</span><span>
</span><span>mean</span><span>(</span><span>df_informative</span><span>$</span><span>diff</span><span> </span><span>&lt;=</span><span> </span><span>0</span><span>)</span><span>
</span><span>#> [1] 0.15525
</span></code></pre>
</div>
<p>It’s also worth comparing the two models. I’ve recently become a fan of the<br>
<a target="_blank" rel="nofollow" href="http://xavier-fim.net/packages/ggmcmc/">ggmcmc package</a> for quick visualization<br>
of Stan models. The package uses a function <code>ggs()</code> to create a long dataframe<br>
of MCMC samples. Then you plug those dataframes into various plotting functions<br>
that start with <code>ggs_</code>. I especially like how the package returns a plain<br>
ggplot2 plot that I can easily adjust with a few extra lines of code.</p>
<p>For example, not much effort is required—after some practice and<br>
trial-and-error, of course—to visualize the posterior samples in each model.</p>
<div>
<pre><code><span>library</span><span>(</span><span>ggmcmc</span><span>)</span><span>
</span><span>#> Loading required package: tidyr
#> #> Attaching package: 'tidyr'
#> The following object is masked from 'package:rstan':
#> #> extract
</span><span>
</span><span># A helper dataframe for relabeling parameters. I'm writing them in a way that
# works with ?plotmath conventions.
</span><span>labels</span><span> </span><span>&lt;-</span><span> </span><span>data.frame</span><span>(</span><span>
</span><span>Parameter</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>&quot;theta_1&quot;</span><span>,</span><span> </span><span>&quot;theta_2&quot;</span><span>,</span><span> </span><span>&quot;diff&quot;</span><span>),</span><span> </span><span>Label</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>&quot;theta[male]&quot;</span><span>,</span><span> </span><span>&quot;theta[female]&quot;</span><span>,</span><span> </span><span>&quot;theta[male] - theta[female]&quot;</span><span>))</span><span>
</span><span># Get ggmcmc's tidy dataframe of each model.
# ggs() doesn't like that labels I made have brackets so I am suppressing its
# warnings.
</span><span>ggs_flat</span><span> </span><span>&lt;-</span><span> </span><span>suppressWarnings</span><span>(</span><span>
</span><span>ggs</span><span>(</span><span>m_flat</span><span>,</span><span> </span><span>description</span><span> </span><span>=</span><span> </span><span>&quot;flat&quot;</span><span>,</span><span> </span><span>par_labels</span><span> </span><span>=</span><span> </span><span>labels</span><span>)</span><span>
</span><span>)</span><span>
</span><span>ggs_informative</span><span> </span><span>&lt;-</span><span> </span><span>suppressWarnings</span><span>(</span><span>
</span><span>ggs</span><span>(</span><span>m_informative</span><span>,</span><span> </span><span>description</span><span> </span><span>=</span><span> </span><span>&quot;informative&quot;</span><span>,</span><span> </span><span>par_labels</span><span> </span><span>=</span><span> </span><span>labels</span><span>)</span><span>
</span><span>)</span><span>
</span><span>ggs_density</span><span>(</span><span>ggs_flat</span><span>)</span><span> </span><span>+</span><span> </span><span>facet_grid</span><span>(</span><span>Parameter</span><span> </span><span>~</span><span> </span><span>.</span><span>,</span><span> </span><span>labeller</span><span> </span><span>=</span><span> </span><span>label_parsed</span><span>)</span><span> </span><span>+</span><span> </span><span>ggtitle</span><span>(</span><span>&quot;flat prior: beta(1, 1)&quot;</span><span>)</span><span> </span><span>+</span><span> </span><span>theme_grey</span><span>(</span><span>base_size</span><span> </span><span>=</span><span> </span><span>14</span><span>)</span><span> </span><span>+</span><span> </span><span>theme</span><span>(</span><span>legend.position</span><span> </span><span>=</span><span> </span><span>&quot;bottom&quot;</span><span>)</span><span> </span><span>+</span><span> </span><span># so the two models can be compared
</span><span> </span><span>xlim</span><span>(</span><span>-.2</span><span>,</span><span> </span><span>.5</span><span>)</span><span>
</span><span>ggs_density</span><span>(</span><span>ggs_informative</span><span>)</span><span> </span><span>+</span><span> </span><span>facet_grid</span><span>(</span><span>Parameter</span><span> </span><span>~</span><span> </span><span>.</span><span>,</span><span> </span><span>labeller</span><span> </span><span>=</span><span> </span><span>label_parsed</span><span>)</span><span> </span><span>+</span><span> </span><span>ggtitle</span><span>(</span><span>&quot;informative prior: beta(5, 40)&quot;</span><span>)</span><span> </span><span>+</span><span> </span><span>theme_grey</span><span>(</span><span>base_size</span><span> </span><span>=</span><span> </span><span>14</span><span>)</span><span> </span><span>+</span><span> </span><span>theme</span><span>(</span><span>legend.position</span><span> </span><span>=</span><span> </span><span>&quot;bottom&quot;</span><span>)</span><span> </span><span>+</span><span> </span><span># so the two models can be compared
</span><span> </span><span>xlim</span><span>(</span><span>-.2</span><span>,</span><span> </span><span>.5</span><span>)</span><span>
</span></code></pre>
</div>
<p><img title="Density plot of the MCMC samples for the parameters in each model." alt="Density plot of the MCMC samples for the parameters in each model." src="https://i0.wp.com/tjmahr.github.io/figs//2017-05-16-bayesian-fisher-exact-test/ggmcmc-density-1.png?w=50%25&ssl=1" data-recalc-dims="1"><img title="Density plot of the MCMC samples for the parameters in each model." alt="Density plot of the MCMC samples for the parameters in each model." src="https://i0.wp.com/tjmahr.github.io/figs//2017-05-16-bayesian-fisher-exact-test/ggmcmc-density-2.png?w=50%25&ssl=1" data-recalc-dims="1"></p>
<p>We can also compare the models together in a single plot by passing a list of<br>
model dataframes into <code>ggs_caterpillar()</code>.</p>
<div>
<pre><code><span>ggs_caterpillar</span><span>(</span><span>D</span><span> </span><span>=</span><span> </span><span>list</span><span>(</span><span>ggs_flat</span><span>,</span><span> </span><span>ggs_informative</span><span>),</span><span> </span><span>line</span><span> </span><span>=</span><span> </span><span>0</span><span>,</span><span>
</span><span>thick_ci</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>0.05</span><span>,</span><span> </span><span>0.95</span><span>),</span><span> </span><span>thin_ci</span><span> </span><span>=</span><span> </span><span>c</span><span>(</span><span>0.025</span><span>,</span><span> </span><span>0.975</span><span>))</span><span> </span><span>+</span><span>
</span><span># Parse the labels as formatted math
</span><span> </span><span>scale_y_discrete</span><span>(</span><span>
</span><span>breaks</span><span> </span><span>=</span><span> </span><span>as.character</span><span>(</span><span>unique</span><span>(</span><span>ggs_flat</span><span>$</span><span>Parameter</span><span>)),</span><span>
</span><span>labels</span><span> </span><span>=</span><span> </span><span>parse</span><span>(</span><span>text</span><span> </span><span>=</span><span> </span><span>as.character</span><span>(</span><span>unique</span><span>(</span><span>ggs_flat</span><span>$</span><span>Parameter</span><span>))))</span><span> </span><span>+</span><span>
</span><span>labs</span><span>(</span><span>caption</span><span> </span><span>=</span><span> </span><span>&quot;Intervals: thick 90%, thin 95%. Point: median.&quot;</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>NULL</span><span>,</span><span> </span><span>x</span><span> </span><span>=</span><span> </span><span>NULL</span><span>)</span><span> </span><span>+</span><span> </span><span>theme_grey</span><span>(</span><span>base_size</span><span> </span><span>=</span><span> </span><span>14</span><span>)</span><span>
</span></code></pre>
</div>
<p><img title="Caterpillar plot of the MCMC samples for the parameters in each model." alt="Caterpillar plot of the MCMC samples for the parameters in each model." src="https://i0.wp.com/tjmahr.github.io/figs//2017-05-16-bayesian-fisher-exact-test/ggmcmc-caterpillar-1.png?w=80%25&ssl=1" data-recalc-dims="1"></p>
<p>There’s a lot of useful information here. First, the intervals in the flat prior<br>
model are wider than the ones for the informative model. The two models largely<br>
agree on the values in the female group, although the flat prior model is wider.<br>
These wider intervals indicate greater uncertainty about the parameter values.<br>
The two models disagree on the male group, because the informative model assigns<br>
little prior probability to values greater than .25 but the flat model doesn’t<br>
discount those possibilities.</p>
<p>The models also demonstrate the <strong>regularizing effect of prior information</strong>.<br>
Regularization broadly refers to techniques to avoid overfitting a dataset.<br>
Priors can regularize a model by making it skeptical of certain parameter<br>
values—in this case, high values of left-handedness. The male probabilities in<br>
the flat model are basically pulled towards the values in the informative prior<br>
(.1-ish). We can see this effect in how the midpoint is shifted in the flat<br>
model versus the informative model.</p>
<hr>
<p>I would be remiss if I didn’t end with the following disclaimer/trivia. In a<br>
way, the question behind this post is ill-posed because handedness is not quite<br>
a binary measure. Some years ago, I had a class on stuttering and fluency<br>
disorders. (I used to be a speech pathologist.) There once was a lot of research<br>
on the association between handedness and stuttering, and at some point,<br>
researchers figured out that they could measure handedness as a continuous<br>
measure. They gave people a survey asking which hands they use for certain tasks<br>
and then computed a so-called <a target="_blank" rel="nofollow" href="https://scholar.google.com/scholar?hl=en&q=dextrality+quotient">dextrality quotient</a> from the responses. I<br>
had the concept driven home when I once saw my wife casually switch<br>
between her left and right hands while brushing her teeth. I would probably poke<br>
a hole through my cheek if I used my left hand! So: Handedness may be a matter<br>
of degree. (Or maybe not. I mostly wanted to mention the dextrality quotient.<br>
It’s fun to think about.)</p>
<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://tjmahr.github.io/bayesian-fisher-exact-test/"> Higher Order Functions</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/z5zaa5jyRe4"></a></div>