<div>
<p><a rel="nofollow" href="http://mc-stan.org/" target="_blank">The Stan project</a> for statistical computation has a great collection of <a rel="nofollow" href="http://mc-stan.org/users/documentation/case-studies" target="_blank">curated case studies</a> which anybody can contribute to, maybe even me, I was thinking. But I don’t have time to worry about that right now because I’m on vacation, being on the yearly visit to my old family home in the north of Sweden. </p>
<p>What I <em>do</em> worry about is that my son will be stung by a bumblebee. His name is Torsten, he’s almost two years old, and he loves running around barefoot on the big lawn. Which has its fair share of bumblebees. Maybe I should put shoes on him so he wont step on one, but what are the chances, really.</p>
<p>Well, what <em>are</em> the chances? I guess if I only had </p>
<ul>
<li>Data on the bumblebee density of the lawn.</li>
<li>Data on the size of Torsten’s feet and how many steps he takes when running around.</li>
<li>A reasonable Bayesian model, maybe implemented in Stan.</li>
</ul>
<p>I could figure that out. “How hard can it be?”, I thought. And so I made an attempt.</p>
<h2>Getting the data</h2>
<p>To get some data on bumblebee density I marked out a 1 m² square on a representative part of the lawn. During the course of the day, now and then, I counted up how many bumblebees sat in the square.</p>
<p><img alt="" src="https://i1.wp.com/www.sumsar.net/images/posts/2017-08-14-the-probability-my-son-will-be-stung-by-a-bumblebee/square_on_the_lawn.jpeg?w=456" data-recalc-dims="1"></p>
<p></p>
<p>Most of the time I saw zero bumblebees, but 1 m² is not that large of an area. Let’s put the data into R:</p>
<pre><code><div><pre>bumblebees <span>&lt;-</span> c(<span>1</span>, <span>0</span>,  <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>,
                <span>1</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>1</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>1</span>,
                <span>0</span>, <span>0</span>, <span>1</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>, <span>0</span>)
</pre>
</div>
<p></p>
<p>During the same day I kept an eye on Torsten, and sometimes when he was engaged in <em>active play</em> I started a 60 s. timer and counted how many steps he took while running around. Here are the counts:</p>
<pre><code><div><pre>toddler_steps <span>&lt;-</span> c(<span>26</span>, <span>16</span>, <span>37</span>, <span>101</span>, <span>12</span>, <span>122</span>, <span>90</span>, <span>55</span>, <span>56</span>, <span>39</span>, <span>55</span>, <span>15</span>, <span>45</span>, <span>8</span>)
</pre>
</div>
<p></p>
<p>Finally I needed to know the area of Torsten’s feet. At a moment when he was <em>not</em> running around I put his foot on a piece of graph paper and traced its edge, which made him giggle. I then calculated how many of the 0.5 cm² squares were fully covered and how many were partially covered.</p>
<p><img alt="" src="https://i2.wp.com/www.sumsar.net/images/posts/2017-08-14-the-probability-my-son-will-be-stung-by-a-bumblebee/torstens_foot_area.jpeg?w=456" data-recalc-dims="1"></p>
<p>To estimate of the area of his foot I took the average of the number of full squares and the number partial squares, and then converted to m²:</p>
<pre><code><div><pre>full_squares <span>&lt;-</span> <span>174</span>
partial_squares <span>&lt;-</span> <span>251</span>
squares <span>&lt;-</span> (full_squares <span>+</span> partial_squares) <span>/</span> <span>2</span>
foot_cm2 <span>&lt;-</span> squares <span>/</span> <span>4</span>
foot_m2 <span>&lt;-</span> foot_cm2 <span>/</span> <span>100^2</span>
foot_m2
</pre>
</div>
<p></p>
<pre><code><div><pre><span>## [1] 0.0053125</span>
</pre>
</div>
<p></p>
<p>Turns out my son’s feet are about 0.0053 m² each.</p>
<h2>Specifying the model</h2>
<p>The idea I had here was that <em>if</em> you know the average number of steps Torsten takes per minute while playing, and <em>if</em> you know the average number of bumblebees per m², <em>then</em> you could calculate the average area covered by Torsten’s tiny feet while running around, and then <em>finally</em> you could calculate the probability of a bumblebee being in that area. This would give you a probability for how likely Torsten is to be stung by a Bumblebee. However, while we have data on all of the above, we don’t know any of the above for certain. </p>
<p>So, to capture some of this uncertainty I thought I would whip together a small Bayesian model and, even though it’s a bit overkill in this case, why not fit it using Stan? (If you are unfamiliar with Bayesian data analysis and Stan I’ve made <a rel="nofollow" href="http://www.sumsar.net/blog/2017/02/introduction-to-bayesian-data-analysis-part-one/" target="_blank">a short video tutorial you can find here</a>.)</p>
<p>Let’s start our Stan program by declaring what data we have:</p>
<pre><code><div><pre>data {
  int n_bumblebee_observations;
  int bumblebees[n_bumblebee_observations];
  int n_toddler_steps_observations;
  int toddler_steps[n_toddler_steps_observations];
  real foot_m2;
}
</pre>
</div>
<p></p>
<p>Now we need to decide on a model for the number of <code>bumblebees</code> in a m² and a model for the number of <code>toddler_steps</code> in 60 s. For the bees I’m going with a simple <em>Poisson</em> model, that is, the model assumes there is an average number of bees per m² ($\mu_\text{bees}$) and that this is all there is to know about bees. </p>
<p>For the number of steps I’m going to step it up a notch with a <em>negative binomial</em> model. The negative binomial distribution can be viewed in a number of ways, but one particularly useful way is as an extension of the Poisson distribution where the mean number of steps for each outcome is not fixed, but instead is varying, where this is controlled by the precision parameter $\phi_\text{steps}$. The <em>larger</em> $\phi_\text{steps}$ is the <em>less</em> the mean number of steps for each outcome is going to vary around the overall mean $\mu_\text{steps}$. That is, when $\phi_\text{steps} \rightarrow \infty$ then the negative binomial becomes a Poisson distribution. The point with using the negative binomial is to capture that Torsten’s activity level when playing on the lawn is definitely <em>not</em> constant. Sometimes he can run full speed for minutes, but sometimes he spends a long time being fascinated by the same little piece of rubbish and then he’s not taking many steps. </p>
<p>So we almost have a Bayesian model for our data, we just need to put priors on the parameters $\mu_\text{bees}$, $\mu_\text{steps}$, and $\phi_\text{steps}$. All three parameters are positive and a quick n’ dirty prior when you have positive parameters is a <em>half-Cauchy</em> distribution centered at zero: It’s just like half a normal distribution centered at zero, but with a much fatter tail. It has one free parameter, the scale, which also happens to be the median of the half-Cauchy. Set the scale/median to a good guess for the parameter in question and you are ready to go! If your guess <em>is</em> good then this prior will give the parameter the slightest nudge in the right direction, if your guess is bad then, as long as you have enough data, the fat tail of the half-Cauchy distribution will help you save face.</p>
<p>My guess for $\mu_\text{bees}$ is (4.0 + 4.0) / (10 × 10) = 0.08. I base this on that while crossing the lawn I usually scan an area of about 10×10 m², I then usually see around four bees, and I assume I’m missing half of the bees. So eight bees on a 100 m² lawn gives 0.08 bees per m². My guess for the number of steps per minute is going to be 60. For the precision parameter $\phi_\text{steps}$ I have no clue, but it shouldn’t be too large nor to small, so my guess is going to be 1.0. Here are the priors:</p>
<p><img alt="" src="https://i0.wp.com/www.sumsar.net/figures/2017-08-14-the-probability-my-son-will-be-stung-by-a-bumblebee/unnamed-chunk-4-1.png?w=456" data-recalc-dims="1"></p>
<p>With the priors specified we now have all the required parts and here is the resulting Bayesian model:</p>
<p>$$% <![CDATA[<br />
 \begin{align}<br />
   &\text{bees} \sim \text{Poisson}(\mu_\text{bees}) \<br />
   &\text{steps} \sim \text{Neg-Binomial}(\mu_\text{steps}, \phi_\text{steps}) \<br />
   &\mu_\text{bees} \sim \text{Half-Cauchy}((4 + 4) / (10 \times 10)) \<br />
   &\mu_\text{steps} \sim \text{Half-Cauchy}(60) \<br />
   &\phi_\text{steps} \sim \text{Half-Cauchy}(1)<br />
   \end{align} %]]>$$</p>
<p>And here is the Stan code implementing the model above:</p>
<pre><code><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">parameters {
  real<span style="color: #666666"><</span>lower<span style="color: #666666">=0></span> mu_bees;
  real<span style="color: #666666"><</span>lower<span style="color: #666666">=0></span> mu_steps;
  real<span style="color: #666666"><</span>lower<span style="color: #666666">=0></span> precision_steps;
}

model {
  <span style="color: #408080; font-style: italic"># Since we have contrained the parameters to be positive we get implicit</span>
  <span style="color: #408080; font-style: italic"># half-cauchy distributions even if we declare them to be 'full'-cauchy.</span>
  mu_bees <span style="color: #666666">~</span> cauchy(<span style="color: #666666">0</span>, (<span style="color: #666666">4.0</span> <span style="color: #666666">+</span> <span style="color: #666666">4.0</span>) <span style="color: #666666">/</span> (<span style="color: #666666">10.0</span> <span style="color: #666666">*</span> <span style="color: #666666">10.0</span>) );
  mu_steps <span style="color: #666666">~</span> cauchy(<span style="color: #666666">0</span>, <span style="color: #666666">60</span>);
  precision_steps <span style="color: #666666">~</span> cauchy(<span style="color: #666666">0</span>, <span style="color: #666666">1</span>);
  bumblebees <span style="color: #666666">~</span> poisson(mu_bees);
  toddler_steps <span style="color: #666666">~</span> neg_binomial_2(mu_steps, precision_steps);
}
</pre>
</div>
<p></code></p>
<p>The final step is to predict how many bumblebees Torsten will step on during, say, one hour of active play. We do this in the <code>generated quantities</code> code block in Stan. The code below will step through 60 minutes (a.k.a. one hour) and for each minute: (1) Sample <code>pred_steps</code> from the negative binomial, (2) calculate the area (m²) covered by these steps, and (3) sample the number of bees in this area from the Poisson giving a predicted <code>stings_by_minute</code>. Finally we sum these 60 minutes worth of strings into <code>stings_by_hour</code>.</p>
<pre><code><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">generated quantities {
  int stings_by_minute[<span style="color: #666666">60</span>];
  int stings_by_hour;
  int pred_steps;
  real stepped_area;

  <span style="color: #008000; font-weight: bold">for</span>(minute <span style="color: #008000; font-weight: bold">in</span> <span style="color: #666666">1:60</span>) {
    pred_steps <span style="color: #666666">=</span>  neg_binomial_2_rng(mu_steps, precision_steps);
    stepped_area <span style="color: #666666">=</span> pred_steps <span style="color: #666666">*</span> foot_m2;
    stings_by_minute[minute] <span style="color: #666666">=</span> poisson_rng(mu_bees <span style="color: #666666">*</span> stepped_area);
  }
  stings_by_hour <span style="color: #666666">=</span> sum(stings_by_minute);
}
</pre>
</div>
<p></code></p>
<h2 id="looking-at-the-result">Looking at the result</h2>
<p>We have data, we have a model, and now all we need to do is fit it. After we’ve put the whole Stan model into <code>bee_model_code</code> as a string we do it like this using the <code>rstan</code> interface:</p>
<pre><code><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span style="color: #408080; font-style: italic"># First we put all the data into a list</span>
data_list <span style="color: #666666"><-</span> list(
  n_bumblebee_observations <span style="color: #666666">=</span> length(bumblebees),
  bumblebees <span style="color: #666666">=</span> bumblebees,
  n_toddler_steps_observations <span style="color: #666666">=</span> length(toddler_steps),
  toddler_steps <span style="color: #666666">=</span> toddler_steps,
  foot_m2 <span style="color: #666666">=</span> foot_m2
)
<span style="color: #408080; font-style: italic"># Then we fit the model which description is in bee_model_code .</span>
library(rstan)
fit <span style="color: #666666"><-</span> stan(model_code <span style="color: #666666">=</span> bee_model_code, data <span style="color: #666666">=</span> data_list)
</pre>
</div>
<p></code></p>
<p>With the model fitted, let’s take a look at the posterior probability distribution of our parameters.</p>
<pre><code><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">stan_hist(fit, c(<span style="color: #BA2121">"mu_bees"</span>, <span style="color: #BA2121">"mu_steps"</span>, <span style="color: #BA2121">"precision_steps"</span>))
print(fit, c(<span style="color: #BA2121">"mu_bees"</span>, <span style="color: #BA2121">"mu_steps"</span>, <span style="color: #BA2121">"precision_steps"</span>))
</pre>
</div>
<p></code></p>
<p><img src="https://i2.wp.com/www.sumsar.net/figures/2017-08-14-the-probability-my-son-will-be-stung-by-a-bumblebee/unnamed-chunk-8-1.png?w=456" alt="" data-recalc-dims="1"></p>
<pre><code><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span style="color: #408080; font-style: italic">##                  mean se_mean    sd  2.5%   25%   50%   75% 97.5% n_eff Rhat</span>
<span style="color: #408080; font-style: italic">## mu_bees          0.12    0.00  0.05  0.04  0.08  0.11  0.15  0.23  3551    1</span>
<span style="color: #408080; font-style: italic">## mu_steps        50.88    0.22 11.65 32.56 43.09 49.34 57.13 78.00  2801    1</span>
<span style="color: #408080; font-style: italic">## precision_steps  1.80    0.01  0.68  0.78  1.32  1.70  2.17  3.41  2706    1</span>
</pre>
</div>
<p></code></p>
<p>This seems reasonable, <em>but</em> note that the distributions are pretty wide, which means there is a lot of uncertainty! For example, looking at <code>mu_bees</code> it’s credible that there is everything from 0.04 bees/m² to 0.2 bees/m². </p>
<p>Anyway, let’s take a look at what we are really interested in, the predictive distribution over how many stings per hour Torsten will suffer during active play:</p>
<pre><code><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span style="color: #408080; font-style: italic"># Getting the sample representing the prob. distribution over stings/hour .</span>
stings_by_hour <span style="color: #666666"><-</span> as.data.frame(fit)<span style="color: #666666">$</span>stings_by_hour
<span style="color: #408080; font-style: italic"># Now actually calculating the prob of 0 stings, 1 stings, 2 stiongs, etc.</span>
strings_probability <span style="color: #666666"><-</span> prop.table( table(stings_by_hour) )
<span style="color: #408080; font-style: italic"># Plotin' n' printin'</span>
barplot(strings_probability, col <span style="color: #666666">=</span> <span style="color: #BA2121">"yellow"</span>,
        main <span style="color: #666666">=</span> <span style="color: #BA2121">"Posterior predictive stings per hour"</span>,
        xlab <span style="color: #666666">=</span> <span style="color: #BA2121">"Number of stings"</span>, ylab <span style="color: #666666">=</span> <span style="color: #BA2121">"Probability"</span>)
</pre>
</div>
<p></code></p>
<p><img src="https://i2.wp.com/www.sumsar.net/figures/2017-08-14-the-probability-my-son-will-be-stung-by-a-bumblebee/unnamed-chunk-9-1.png?w=456" alt="" data-recalc-dims="1"></p>
<pre><code><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">round(strings_probability, <span style="color: #666666">2</span>)
</pre>
</div>
<p></code></p>
<pre><code><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span style="color: #408080; font-style: italic">## stings_by_hour</span>
<span style="color: #408080; font-style: italic">##    0    1    2    3    4    5    6    7    8    9   10 </span>
<span style="color: #408080; font-style: italic">## 0.27 0.30 0.22 0.12 0.05 0.02 0.01 0.00 0.00 0.00 0.00</span>
</pre>
</div>
<p></code></p>
<p>Ok, it seems like it is most probable that Torsten will receive one sting per hour, but we should not be surprised if it’s two or even three stings. I’d better put some shoes on him! The <em>problem</em> is that after a couple of days full of active barefoot play, my son Torsten’s feet look like this:</p>
<p><img src="https://i1.wp.com/www.sumsar.net/images/posts/2017-08-14-the-probability-my-son-will-be-stung-by-a-bumblebee/torstens_feet.jpeg?w=456" alt="" data-recalc-dims="1"></p>
<p>As you can see his feet are <em>not</em> swollen from all the bee stings they should receive according to the model. Actually, even after a week, he has not gotten a single bee sting! Which is good, I suppose, in a way, but, well, it means that my model is likely pretty crap. How can this be? Well,</p>
<ul>
<li>I should maybe have gotten more and better data. Perhaps the square I monitored for bumblebees didn’t yield data that was really representative of the bumblebee density in general.</li>
<li>The assumption that bumblebees always sting when stepped upon might be wrong. Or maybe Torsten is so smart so that he actively avoids them…</li>
<li>Maybe the model was too simplistic. I really should have made a hierarchical model incorporating data from multiple squares and several days of running around. To factor in the effect of the weather, the flower density, and the diet of Torsent also couldn’t hurt.</li>
</ul>
<p>I guess I could spend some time improving the model. And I guess there is a lesson to be learned here about that it is hard to predict the predictive performance of a model. But all that has to wait. I am on vacation after all.</p>
<p><em>The full data and code behind this blog post can be found <a href="https://gist.github.com/rasmusab/850e5c76950c3e755ee0293b8552000a" rel="nofollow" target="_blank">here</a>.</em></p>
<p></body></html></p>

		<script type='text/javascript'>
		  var vglnk = { key: '949efb41171ac6ec1bf7f206d57e90b8' };

		  (function(d, t) {
			var s = d.createElement(t); s.type = 'text/javascript'; s.async = true;
			s.src = '//cdn.viglink.com/api/vglnk.js';
			var r = d.getElementsByTagName(t)[0]; r.parentNode.insertBefore(s, r);
		  }(document, 'script'));
		</script>		
		
<p class="syndicated-attribution">
<div style="border: 1px solid; background: none repeat scroll 0 0 #EDEDED; margin: 1px; font-size: 13px;">
<div style="text-align: center;">To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://www.sumsar.net/blog/2017/08/the-probability-my-son-will-be-stung-by-a-bumblebee/"> Publishable Stuff</a></strong>.</div>
<hr />
<a href="https://www.r-bloggers.com/" rel="nofollow">R-bloggers.com</a> offers <strong><a href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers" rel="nofollow">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" href="https://www.r-project.org/" rel="nofollow">R</a> news and <a title="R tutorials" href="https://www.r-bloggers.com/search/tutorial" rel="nofollow">tutorials</a> on topics such as: <a title="Data science" href="https://www.r-bloggers.com/search/data%20science" rel="nofollow">Data science</a>, <a title="Big Data" href="https://www.r-bloggers.com/search/Big%20Data" rel="nofollow">Big Data, <a title="R jobs" href="https://www.r-users.com/" rel="nofollow">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" href="https://www.r-bloggers.com/search/ggplot2" rel="nofollow">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" href="https://www.r-bloggers.com/search/boxplot" rel="nofollow">Boxplots</a>, <a title="Maps and gis" href="https://www.r-bloggers.com/search/map" rel="nofollow">maps</a>, <a title="Animation in R" href="https://www.r-bloggers.com/search/animation" rel="nofollow">animation</a>), programming (<a title="RStudio IDE for R" href="https://www.r-bloggers.com/search/RStudio" rel="nofollow">RStudio</a>, <a title="Sweave and literate programming" href="https://www.r-bloggers.com/search/sweave" rel="nofollow">Sweave</a>, <a title="LaTeX in R" href="https://www.r-bloggers.com/search/LaTeX" rel="nofollow">LaTeX</a>, <a title="SQL and databases" href="https://www.r-bloggers.com/search/SQL" rel="nofollow">SQL</a>, <a title="Eclipse IDE for R" href="https://www.r-bloggers.com/search/eclipse" rel="nofollow">Eclipse</a>, <a title="git and github, Version Control System" href="https://www.r-bloggers.com/search/git" rel="nofollow">git</a>, <a title="Large data in R using Hadoop" href="https://www.r-bloggers.com/search/hadoop" rel="nofollow">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" href="https://www.r-bloggers.com/search/Web+Scraping" rel="nofollow">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" href="https://www.r-bloggers.com/search/regression" rel="nofollow">regression</a>, <a title="principal component analysis tutorial" href="https://www.r-bloggers.com/search/PCA" rel="nofollow">PCA</a>, <a title="Time series" href="https://www.r-bloggers.com/search/time+series" rel="nofollow">time series</a>, <a title="finance trading" href="https://www.r-bloggers.com/search/trading" rel="nofollow">trading</a>) and more...
</div></p><img src="http://feeds.feedburner.com/~r/RBloggers/~4/4aVFWV6lgJw" height="1" width="1" alt=""/>]]></code></pre></code></pre></code></pre></code></pre></code></pre></div>