<p>A few days ago, I watched a <a href="https://www.youtube.com/watch?v=qQFF4tPgeWI">terrific lecture by Bob Carpenter</a> on
Bayesian models. He started with a Bayesian approach to Fisher’s exact test. I
had never heard of this classical procedure, so I was curious to play with the
example. In this post, I use the same data that he used in the lecture and
in an earlier, <a 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
fitting the model in Stan and inspecting the results in R.</p>
<h2>Problem statement</h2>
<p>We observed the following data.</p>
<table><tr>sex
n left handed
n right handed
</tr><tr><td>male</td>
<td>9</td>
<td>43</td>
</tr><tr><td>female</td>
<td>4</td>
<td>44</td>
</tr></table><p>Question: Is the rate of left-handedness different between the male and female
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
the <a 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
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
describing the alternative hypothesis.</p>
<p>We can also consider one-sided test: Is left-handedness greater in the male
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 p-value is
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
stats classes, and indeed, this post is the first time I ever used the function
<code>fisher.test()</code>. If I had never heard of the test, I am not quite sure what I
would have done. Maybe a logistic regression (p = .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
Stan, I have an idea about how the data could have been generated, and that’s good enough :relaxed:. I can just write down my data-generating story in a model
and let Stan compute a posterior distribution for the difference in handedness
rates between the two groups.</p>
<p>For my model, I’m going to suppose that in each group, there is a probability of
being left-handed called θ and that the counts we see result from a binomial process. The 9 left-handed males we observe are the number of successes
from 52 observations of a process that “succeeds” with probability
θmale.</p>
<p>To fit the model, we need a prior distribution. The prior’s job is to generate
possible θ’s, and we will use our data to update the prior. The Beta
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>We can use a flat, uninformative prior by using a = 1, b = 1. This
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></p>
<p>But I also think that 10-ish% of people are left handed. (I don’t know where I first heard this number, but it’ll serve as my prior information.) I toyed around with <code>shape1</code> and <code>shape2</code> parameters in <code>dbeta()</code> until I got the prior Beta(5, 40), which is peaked around .1-ish but wide enough to keep .5 and .15 as
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></p>
<p>Let’s write out a really simple model in Stan. Okay, it used to be really simple. Then I made the parameters for the Beta prior data values, and then I
created an option to just sample the prior distribution. But the core of it is
simple. The most important lines are the one with <code>~</code> symbols. These
correspond to the sampling statements in the mathematical description of the
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
distribution. Here, we compute the
θmale − θfemale on every
draw. Computing the difference inside the model code means that Stan
will treat the <code>diff</code> values like any other parameter of the model. It will show
up in summary functions and in plots of model parameters. That saves us some 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 program that can sample from the model. I do the compilation in its own step so
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. This step is not necessary, but I don’t like repeating myself. (Normally, you would use <code>sampling(model_program, stan_data)</code> to get samples from a 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 adjust the parameters for the Beta prior and toggle between sampling from the prior and the posterior. Stan normally prints out verbose progress information
but I suppress that by using <code>purrr::quietly()</code>. I still want warnings, so I
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 that our program works as expected. We know how the values should be 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 left-handedness in both groups is 0 according to our priors. The
informative prior says that the values between 0.04 and
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.
The central θ value in each group is .5. Therefore, half of the prior samples assert that there are more left-handed individuals than right-handed ones! If there is anything we know about handedness, it’s that left-handedness is
less common than right-handedness. Uninformative sometimes connotes “unbiased”
or “letting the data speak for itself”, but in this case, I would say
“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
values fall in the interval [-0.02, 0.19]. The informative
model is more skeptical of higher left-handedness rates, so it puts the difference
at 0.05 with 90% of the values between [-0.03, 0.13].
Both of these intervals contain 0 and negative values :x:, so there is not much
evidence for higher left-handedness in the male group.</p>
<p>To compute a “Bayesian p-value”, we could ask what proportion of differences are 0 or negative. There are more proper ways to make this inference in a
Bayesian framework, but this approach is the easiest and it works for a model
this simple. If 10% of the plausible values for the group differences are
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
<a href="http://xavier-fim.net/packages/ggmcmc/">ggmcmc package</a> for quick visualization
of Stan models. The package uses a function <code>ggs()</code> to create a long dataframe
of MCMC samples. Then you plug those dataframes into various plotting functions
that start with <code>ggs_</code>. I especially like how the package returns a plain
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
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></p>
<p>We can also compare the models together in a single plot by passing a list of
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></p>
<p>There’s a lot of useful information here. First, the intervals in the flat prior
model are wider than the ones for the informative model. The two models largely
agree on the values in the female group, although the flat prior model is wider.
These wider intervals indicate greater uncertainty about the parameter values.
The two models disagree on the male group, because the informative model assigns
little prior probability to values greater than .25 but the flat model doesn’t
discount those possibilities.</p>
<p>The models also demonstrate the <strong>regularizing effect of prior information</strong>.
Regularization broadly refers to techniques to avoid overfitting a dataset.
Priors can regularize a model by making it skeptical of certain parameter
values—in this case, high values of left-handedness. The male probabilities in
the flat model are basically pulled towards the values in the informative prior
(.1-ish). We can see this effect in how the midpoint is shifted in the flat
model versus the informative model.</p>
<hr><p>I would be remiss if I didn’t end with the following disclaimer/trivia. In a way, the question behind this post is ill-posed because handedness is not quite a binary measure. Some years ago, I had a class on stuttering and fluency disorders. (I used to be a speech pathologist.) There once was a lot of research
on the association between handedness and stuttering, and at some point, researchers figured out that they could measure handedness as a continuous measure. They gave people a survey asking which hands they use for certain tasks
and then computed a so-called <a href="https://scholar.google.com/scholar?hl=en&q=dextrality+quotient">dextrality quotient</a> from the responses. I had the concept driven home when I once saw my wife casually switch between her left and right hands while brushing her teeth. I would probably poke
a hole through my cheek if I used my left hand! So: Handedness may be a matter of degree. (Or maybe not. I mostly wanted to mention the dextrality quotient.
It’s fun to think about.)</p>