---
title: "Computer lab: Bayesian analysis"
subtitle: "Methods of Scientific Working (3502-450)"
author: "Prof. Karl Schmid, University of Hohenheim"
date: "`r Sys.Date()`" 
format:
  html:
    embed-resources: true
bibliography: ../scimeth.bib
---

```{r}
library(ggplot2)
```


# Background {#sec-1}

The purpose of this computer lab is to simulate an urn experiment to demonstrate the key principles of a Bayesian analysis.

Then, as an application of the theory we will apply Bayesian reasoning to estimate a parameter (allele frequency) in a simple population genetic model.

The third examples is a Bayesian test of a hypothesis given the data observed in a sample. We will test the hypothesis that that genotype frequencies of a genetic locus are in Hardy Weinberg equilibrium.

# Marble experiment

We first start with a simple marble experiment in which we compare two hypotheses about the number of blue and white marbles in an urn.

## Design of the experiment {#sec-2-1}

The experiment is set up by flipping a fair coin to select one of two possible hypothesis as true hypothesis.

Depending on the outcome, we put the marbles into the urn:

 - If head, place in an urn 1 white and 3 blue marbles
 - If tail, place in an urn 3 white and 1 blue marble

From this, we generate two hypotheses, we will test:

-   H$_B$: 1 white and 3 blue marbles (WBBB)
-   H$_W$: 3 white and 1 blue marble (WWWB)

The goal of the Bayesian analysis is to test which of the hypotheses is probably true (or more likely true) given the data.

To find out, we conduct an experiment: 

1. Mix the marbles in the urn
2. Draw a marble and record its color
3. Put it back into the urn
4. Repeat the procedure as necessary

Such an experiment can be easily set up with R

The goal of the setup is to find out which of the two hypotheses have a higher probability given the data. We flip the coin first to randomly choose between one of the hypotheses. Because we do not know the outcome of the coin toss, we assume that both hypotheses are equally likely and for this reason the prior probability is 1:1 for both hypotheses.

```{r}
coin <- c("head", "tail")
draw <- sample(coin, 1)
if (draw == "head") {
  urn <- c('W', 'B', 'B', 'B')
} else {
  urn <- c('W', 'W', 'W', 'B')
}
```

We have now randomly constructed an urn, from which we can sample. To run our experiment, we need to calculate the posterior probability.  Since we are comparing two hypotheses, we can use the ratio form of Bayes's rule,

$$\frac{P(H_B | D)}{P(H_W | D)} = \frac{P(D | H_B)}{P(D | H_W)} \times \frac{P(H_B)}{P(H_W)}$$

The first term describes the ratio of posterior probabilities, the second the likelihood odds ratio, and the third the ratio of the prior probabilities.

It is now possible to construct the likelihoods. If hypothesis $H_B$ is true, the probability of obtaining a blue draw is 0.75 (three out of four marbles are blue). Under hypothesis $H_W$ it is 0.25. Therefore the likelihood odds ratio is

$$\frac{P(D|H_B)}{P(D|H_W)} = \frac{0.75}{0.25} = \frac{3}{1}.$$ 

If we draw a white marble, the likelihood odds ratio is

$$\frac{P(D|H_B)}{P(D|H_W)} = \frac{0.25}{0.75} = \frac{1}{3}.$$ 

For the first draw, the prior probability of both hypotheses is $P(H_B) = P(H_W) = 0.5$, and we can calculate the posterior probabilities if we drew a blue marble as

$$\frac{P(H_B | D)}{P(H_W | D)} = \frac{0.75}{0.25} \times \frac{0.5}{0.5} = \frac{3}{1}$$

The posterior probability of the first draw is the prior probability of the second draw. Since we are only interested in the ratio of the new prior probabilities, we can just set 

$$\frac{P(H_B)}{P(H_W)}=\frac{3}{1}$$ 

as the new prior ratio and start from new.

If the second draw is a white marble, the posterior ratio is then

$$\frac{P(H_B | D)}{P(H_W | D)} = \frac{0.25}{0.75} \times \frac{3}{1} = \frac{0.75}{0.75} = \frac{1}{1}$$

and so on. We can implement this experiment now in R:

```{r}
    urnB <- c('W','B','B','B')
    urnW <- c('W','W','W','B')

    priorB <- 0.5 # prior probability of hypothesis B
    priorW <- 0.5 # prior probability of hypothesis W
    
    # randomly select the urn (i.e., the true hypothesis)
    hypothesis <- sample(c('H_B', 'H_W'), 1)
    if (hypothesis == 'H_B') { 
      urn <- c('W','B','B','B') 
    } else if (hypothesis == 'H_W')  {
      urn <- c('W','W','W','B') 
    }
    print(urn)
    stoppingrule <- 0.999999

    # set up experiment
    i <- 0
    draws <- c(i)
    results <- c("(Prior)")
    ratio <- c("1:1")
    posteriorB <- c(priorB)
    posteriorW <- c(priorW)
    
    while (priorB < stoppingrule && priorW < stoppingrule) {
      draw <- sample(urn, 1) # Randomly choose the true hypothesis
      if (draw == 'B') {
        likelihoodB <- 0.75
        likelihoodW <- 0.25
      } else {
        likelihoodB <- 0.25
        likelihoodW <- 0.75
      }
      
      # Bayes' theorem for odds ratios
      posteriorRatio <- (likelihoodB / likelihoodW) * (priorB / priorW)
      posteriorInteger <- paste(as.integer(posteriorRatio), ": 1")
      priorB <- posteriorRatio / (posteriorRatio + 1) #see footnote
      priorW <- 1 - priorB
      
      i = i + 1 # increase the counter for the next draw
      
      draws <- c(draws, i)
      results <- c(results, draw)
      ratio <- c(ratio, posteriorInteger)
      posteriorB <- c(posteriorB, priorB)
      posteriorW <- c(posteriorW, priorW)
    }

# put the data into a dataframe and plot the outcome
df = data.frame(draws, H_B = posteriorB, H_W = posteriorW)  

# Melting the dataframe to make it suitable for ggplot
df_long <- reshape2::melt(df, id.vars = 'draws', measure.vars = c('H_W', 'H_B'))

# Plotting
ggplot(df_long, aes(x = draws, y = value, color = variable)) +
  geom_line() +
  geom_point() +
  labs(title = paste("True hypothesis:", hypothesis, "; ", nrow(df)-1, "draws until posterior probability >", stoppingrule), x = "Draw number", y = "Posterior probability", color = "Hypothesis") +
  ylim(0, 1) +
  scale_color_manual(values = c("H_B" = "blue", "H_W" = "black")) 

```


## Excercises

1. Run a single simulation run using this code and inspect the output
2. Repeat the simulation run. How does the output change?
3. Change the prior to check out how much the experiment depends on the prior used in the early stage of the experiment
4. The ratio of the posterior probabilities give the *relative* believe in each of the hypotheses. Check how the ratio is effected when the true hypothesis differs from the one which are tested.

## Study questions {#sec-2-6}

1.  Why can we use the posterior probability of the last experiment as the prior of the next experiment?
2.  How would the calculations be changed if we have three different marble colors instead of two different marbles?

:::{.content-hidden .comment}

# (HOMEWORK) A slightly expanded Bayesian analysis of the urn experiment

The above example served to demonstrate how a basic Bayesian analysis works.

In the following example, we investigate all possible hypotheses for these urn experiment and then ask, which hypothesis is most likely with a given number of draws.

We first decribe the hypotheses:

```{r}
h1 <- c('W', 'W', 'W', 'W')
h2 <- c('B', 'W', 'W', 'W')
h3 <- c('B', 'B', 'W', 'W')
h4 <- c('B', 'B', 'B', 'W')
h5 <- c('B', 'B', 'B', 'B')
hypotheses <- list(h1, h2, h3, h4, h5)
priors <- c(0.2, 0.2, 0.2, 0.2, 0.2) # set prior probability for each hypothesis
likelihoodsBlue <- c(0, 0.25, 0.5, 0.75, 1) # likelihood to draw a blue ball under
                                           # each hypothesis
```

Then we randomly select one hypothesis

```{r}
urn <- unlist(sample(hypotheses, 1))
urn
```

Then make a sample of draws with replacement from the urn

```{r}
drawNumber <- 3 # number of draws in experiment
draws <- sample(urn, replace= TRUE, drawNumber)
```

To calculate the likelihood of a series of draws given a hypothesis, one has to multiply the individual probabilities of the observed result of a draw.

```{r}
calculateTotalLikelihood <- function(draws, likelihoodBlue) {
  totalLikelihood <- 1
  for (draw in draws) {
    if (draw == "B") {
      totalLikelihood <- totalLikelihood * likelihoodBlue
    } else {
      totalLikelihood <- totalLikelihood * (1 - likelihoodBlue)
    } 
  }
  return(totalLikelihood)
}

# Iterate through all hypotheses
totalLikelihoods <- rep(0, 5)
for (hypothesis in seq(1:5)) {
  totalLikelihoods[hypothesis] <- calculateTotalLikelihood(draws, likelihoodsBlue[hypothesis])
}

# Calculate product of likelihood with prior
priorXlikelihoods <- totalLikelihoods * priors
priorXlikelihoods

sumLikelihood <- sum(priorXlikelihoods)
sumLikelihood
posteriorProbabilities <- priorXlikelihoods / sumLikelihood
round(posteriorProbabilities, 3)
```

Finally, plot the posterior probabilities for visualization

```{r}
library(ggplot2)
hypothesisNames <- c('WWWW', 'BWWW', 'BBWW', 'BBBW', 'BBBB')
df <- data.frame(hypothesisNames, posteriorProbabilities)
p<-ggplot(data = df, aes(x = hypothesisNames, y = posteriorProbabilities)) +
  geom_bar(stat = "identity") +
  labs(title = paste("True hypothesis:", paste(urn, collapse = ""), ", Draw number", drawNumber ))
p
```

## Exercises

1. Increase the number of draws to 100, 1000, 10000. What do you observe in the distribution of posterior probabilities? Which computational problem arises with the calculation of likelihoods with very large draw sizes? Do you have an idea how to solve this issue?
2. Choose a certain hypothesis and change the prior probability to different values. How does it effect the distribution of the posterior probabilities?

## Discussion and thinking questions

1. If you reflect the two exercises above, can you name advantages and disadvantages of a Bayesian versus a frequentist analysis of an urn experiment? How would you describe the key differences between the frequentist and the Bayesian analysis?
2. How can one model the Bayes theorem with an infinite number of possible hypotheses? Under which experimental circumstances do you expect an infinite number of hypotheses?
3. Why is the order of the marbles in the draws not important?
4. Bayesian analysis allows to estimate model parameters. What is the single model parameter of the urn experiment, and why does the Bayesian analysis estimate the model parameter.
5. Assume that in the design of a real experiment you have a limited budget. Which parameter in the urn experiment could be influenced by a limited budget?


# Further Reading

The book by @kruschke2014doing is an excellent introduction into Bayesian analysis for non-mathematicians. Check out also his webpage ([www.indiana.edu/\~kruschke](www.indiana.edu/~kruschke)) which contains further introductory papers and a blog on Bayesian statistics.

In genetics, Bayesian methods are widely used. These two reviews are very good introductions into Bayesian thinking and methods for genetics:

- @beaumont_bayesian_2004
- @shoemaker_bayesian_1999

# References
:::


