<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/bayesian-blood/"></div></div><div><a data-url="https://www.r-bloggers.com/bayesian-blood/" data-counturl="https://www.r-bloggers.com/bayesian-blood/" data-count="horizontal" data-text="Bayesian Blood" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://fronkonstin.com/2016/11/17/bayesian-blood/"> R – Fronkonstin</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<div class="twitter-share"><a target="_blank" class="twitter-share-button" rel="nofollow" href="https://twitter.com/intent/tweet?url=https%3A%2F%2Ffronkonstin.com%2F2016%2F11%2F17%2Fbayesian-blood%2F&via=aschinchon">Tweet</a></div>
<blockquote><p>The fourth, the fifth, the minor fall and the major lift (Hallelujah, Leonard Cohen)</p>
</blockquote>
<p>Next problem is extracted from MacKay’s <em>Information Theory, Inference and Learning Algorithms</em>:</p>
<p><em>Two people have left traces of their own blood at the scene of a crime. A suspect, Oliver, is tested and found to have type ‘O’ blood. The blood groups of the two traces are found to be of type ‘O’ (a common type in the local population, having frequency 60%) and of type ‘AB’ (a rare type, with frequency 1%). Do these data give evidence in favor of the proposition that Oliver was one of the people who left blood at the scene?</em></p>
<p>To answer the question, let’s first remember the probability form of Baye’s theorem:</p>
<p><img title="p(H|D)=\dfrac{p(H)p(D|H)}{p(D)}" alt="p(H|D)=\dfrac{p(H)p(D|H)}{p(D)}" class="latex" src="https://s0.wp.com/latex.php?latex=p%28H%7CD%29%3D%5Cdfrac%7Bp%28H%29p%28D%7CH%29%7D%7Bp%28D%29%7D&bg=ffffff&%23038;fg=000&%23038;s=0"></p>
<p>where:</p>
<ul>
<li>p(H) is the probability of the hypothesis H before we see the data, called the <strong>prior</strong></li>
<li>p(H|D) is the probablity of the hyothesis after we see the data, called the <strong>posterior</strong></li>
<li>p(D|H) is the probability of the data under the hypothesis, called the <strong>likelihood</strong></li>
<li>p(D)is the probability of the data under any hypothesis, called the <strong>normalizing constant</strong></li>
</ul>
<p>If we have two hypothesis, A and B, we can write the ratio of posterior probabilities like this:</p>
<p><img title="\dfrac{p(A|D)}{p(B|D)}=\dfrac{p(A)p(D|A)}{p(B)p(D|B)}" alt="\dfrac{p(A|D)}{p(B|D)}=\dfrac{p(A)p(D|A)}{p(B)p(D|B)}" class="latex" src="https://s0.wp.com/latex.php?latex=%5Cdfrac%7Bp%28A%7CD%29%7D%7Bp%28B%7CD%29%7D%3D%5Cdfrac%7Bp%28A%29p%28D%7CA%29%7D%7Bp%28B%29p%28D%7CB%29%7D&bg=ffffff&%23038;fg=000&%23038;s=0"></p>
<p>If p(B)=1-p(B) (what means that A and B are <em>mutually exclusive</em> and <em>collective exhaustive</em>), then we can rewrite the ratio of the priors and the ratio of the posteriors as odds. Writing o(A) for odds in favor of A, we get the <strong>odds form</strong> of Baye’s theorem:</p>
<p><img title="o(A|D)=o(A)\dfrac{p(D|A)}{p(D|B)}" alt="o(A|D)=o(A)\dfrac{p(D|A)}{p(D|B)}" class="latex" src="https://s0.wp.com/latex.php?latex=o%28A%7CD%29%3Do%28A%29%5Cdfrac%7Bp%28D%7CA%29%7D%7Bp%28D%7CB%29%7D&bg=ffffff&%23038;fg=000&%23038;s=0"></p>
<p>Dividing through by o(A) we have:</p>
<p><img title="\dfrac{o(A|D)}{o(A)}=\dfrac{p(D|A)}{p(D|B)}" alt="\dfrac{o(A|D)}{o(A)}=\dfrac{p(D|A)}{p(D|B)}" class="latex" src="https://s0.wp.com/latex.php?latex=%5Cdfrac%7Bo%28A%7CD%29%7D%7Bo%28A%29%7D%3D%5Cdfrac%7Bp%28D%7CA%29%7D%7Bp%28D%7CB%29%7D&bg=ffffff&%23038;fg=000&%23038;s=0"></p>
<p>The term on the left is the ratio of the posteriors and prior odds. The term on the right is the likelihood ratio, also called the <strong>Bayes factor</strong>. If it is greater than 1, that means that the data were more likely under A than under B. And since the odds ratio is also greater than 1, that means that the odds are greater, in light of the data, than they were before. If the Bayes factor is less than 1, that means the data were less likely under A than under B, so th odds in favor of A go down.</p>
<p>Let’s go back to our initial problem. If Oliver left his blood at the crime scene, the probability of the data is just the probability that a random member of the population has type ‘AB’ blood, which is 1%. If Oliver did not leave blood at the scene, what is the the chance of finding two people, one with type ‘O’ and one with type ‘AB’? There are two ways it might happen: the first person we choose might have type ‘O’ and the second ‘AB’, or the other way around. So the probability in this case is 2(0.6)(0.01)=1.2%. Dividing probabilities of both scenarios we obtain a Bayes factor of 0.83, and we conclude that the blood data is evidence <strong>against Oliver’s guilt</strong>.</p>
<p>Once I read this example, I decided to replicate it using real data of blood type distribution by country from <a target="_blank" rel="nofollow" href="http://www.rhesusnegative.net/themission/bloodtypefrequencies/">here</a>. After cleaning data, I have this <em>nice</em> data set to work with:</p>
<p><br>
<iframe height="500" width="100%" src="https://fronkonstin.com/wp-content/uploads/2016/11/table.html"></iframe></p>
<p>For each country, I get <strong>the most common blood type</strong> (the one which the suspect has) <strong>and the least common</strong> and replicate the previous calculations. For example, in the case of Spain, the most common type is ‘O+’ with 36% and the least one is ‘AB-‘ with 0.5%. The Bayes factor is 0.005/(2(0.36)(0.005))=1.39 so data <em>support</em> the hypothesis of guilt in this case. Next chart shows Bayes factor accross countries:</p>
<p><br>
<iframe height="500" width="100%" src="https://fronkonstin.com/wp-content/uploads/2016/11/blood.html"></iframe></p>
<p>Just some comments:</p>
<ul>
<li>Sometimes data consistent with a hypothesis are not necessarily in favor of the hypothesis</li>
<li>How different is the distribution of blood types between countries!</li>
<li>If you are a estonian ‘A+’ murderer, choose carefully your accomplice</li>
</ul>
<p>This is the code of the experiment:</p>
<pre>
library(rvest)
library(dplyr)
library(stringr)
library(DT)
library(highcharter)

# Webscapring of the table with the distribution of blood types
url &lt;- &quot;http://www.rhesusnegative.net/themission/bloodtypefrequencies/&quot;
blood &lt;- url %>%
   read_html() %>%
   html_node(xpath='/html/body/center/table') %>%
   html_table(fill=TRUE)

# Some data cleansing
blood %>% slice(-c(64:68)) -> blood

blood[,-c(1:2)] %>% 
  sapply(gsub, pattern=&quot;,&quot;, replacement=&quot;.&quot;) %>% 
  as.data.frame %>% 
  sapply(gsub, pattern=&quot;.79.2&quot;, replacement=&quot;.79&quot;) %>% 
  as.data.frame-> blood[,-c(1:2)]

blood %>% 
  sapply(gsub, pattern=&quot;%|,&quot;, replacement=&quot;&quot;) %>% 
  as.data.frame -> blood

blood[,-1] = apply(blood[,-1], 2, function(x) as.numeric(as.character(x)))


blood[,-c(1:2)] %>% mutate_all(funs( . / 100)) -> blood[,-c(1:2)]

# And finally, we have a nice data set
datatable(blood, 
          rownames = FALSE,
          options = list(
          searching = FALSE,
          pageLength = 10)) %>% 
  formatPercentage(3:10, 2)

# Calculate the Bayes factor
blood %>% 
  mutate(factor=apply(blood[,-c(1,2)], 1, function(x) {min(x)/(2*min(x)*max(x))})) %>% 
  arrange(factor)-> blood

# Data Visualization
highchart() %>% 
     hc_chart(type = &quot;column&quot;) %>% 
     hc_title(text = &quot;Bayesian Blood&quot;) %>%
     hc_subtitle(text = &quot;An experiment about the Bayes Factor&quot;) %>%  
     hc_xAxis(categories = blood$Country, 
             labels = list(rotation=-90, style = list(fontSize = &quot;12px&quot;)))  %>% 
     hc_yAxis(plotBands = list(list(from = 0, to = 1, color = &quot;rgba(255,215,0, 0.8)&quot;))) %>% 
     hc_add_series(data = blood$factor,
                   color = &quot;rgba(255, 0, 0, 0.5)&quot;,
                   name = &quot;Bayes Factor&quot;)%>% 
  hc_yAxis(min=0.5) %>% 
  hc_tooltip(pointFormat = &quot;{point.y:.2f}&quot;) %>% 
  hc_legend(enabled = FALSE) %>% 
  hc_exporting(enabled = TRUE) %>%
  hc_chart(zoomType = &quot;xy&quot;)
</pre>
<p><div><div><div><div data-width="55" data-show-faces="false" data-send="true" data-layout="box_count" data-height="62" data-href="https://www.r-bloggers.com/bayesian-blood/"></div></div><div><a data-url="https://www.r-bloggers.com/bayesian-blood/" data-counturl="https://www.r-bloggers.com/bayesian-blood/" data-count="vertical" data-text="Bayesian Blood" class="twitter-share-button" data-via="rbloggers" href="https://twitter.com/share"></a></div><div></div></div><div></div></div>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://fronkonstin.com/2016/11/17/bayesian-blood/"> R – Fronkonstin</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/?s=tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/?s=data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/?s=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/?s=ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/?s=map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/?s=animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/?s=RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/?s=sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/?s=LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/?s=SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/?s=eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/?s=git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/?s=hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/?s=Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/?s=regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/?s=PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/?s=time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/?s=trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/svgsfmsBYgc"></a></div>