Approximate Bayesian Computation

• In the previous section, we ran simulations under two demographic models.
• How can we say statistically which model fits the data better?

  

Basic idea of ABC

• Parameter:
  • Property of the population (not the samples)
  • Unknown in real life
  • Examples of parameters:
    • The probability of getting heads in a biased coin flip.
    • Proportion of Cherry flavor in Jelly Beans.
    • $\theta (=4 N \mu)$
    • Time of divergence between two populations.
• ABC approach to estimate the parameters

  

  1. Draw parameters randomly (say divergence time, T between 1000 and 10000 years, and $\theta$ between 0.001 and 1).
  2. Run a simulation with the parameters, and record the result
  3. Goto step 1, and repeat this to create 1 million sets of parameters and simulation results.
  4. Compare the results of simulations with observation (e.g., summary statistics vector).
  5. Select top 500 best matches.
  6. Create the histogram/density plots with the parameter values, which created the best matching simulations.
  7. Or if you are comparing the two hypotheses, count how many times each model is included in the top choices.

Bayesian inference

• Parameter: p = probability of getting a head from a biased coin.
• You flip the coin $n=10$ times, and observed $d=3$ heads (success).
• likelihood
  • If you happened to know the parameter (p), you know the probability of observing this data.
    Assuming that the actual parameter is p, the conditional probability of observing the data (d) is:

    $Pr[data \vert parameter] = {n\choose d} p^d (1-p)^{n-d}$

  • This conditional probability is called likelihood of parameters (or a hypothesis) given data ($L[p\vert data]$). A value of p which maximizes this conditional probability is called the maximum likelihood estimate (MLE) of the parameter p.

    $L[p\vert d=3] = Pr[d=3 \vert p] = {10 \choose 3} p^3 (1-p)^{10-3}$

    

• Won't it be nice if we can actually obtain the probability of parameter p under the observed data (i.e $Pr[p\vert d]$)?

• General Bayes' Theorem:

  $$Pr[A\vert B] = \frac{Pr[B\vert A] Pr[A]}{Pr[B]}$$

  • Simple probability example:
    • 40 % single
      • 50 % normal hearing
      • 50 % selective hearing
    • 60 % married
      • 100 % selective hearing

  • Question: obtain the probability that a guy with selective hearing is single.

  • Hint:
    $Pr[A]$: Pr[a random guy is single regardless of info about hearing]
    $Pr[B\vert A]$: Pr[selective hearing given single]
    $Pr[B]$: Pr[a random guy having selective hearing]
• Bayesian Inference:

  $$Pr[p \vert d=3] = \frac{Pr[d = 3 \vert p] \cdot Pr[p]}{Pr[d=3... ...= 3 \vert p] \cdot Pr[p]}{\sum_p Pr[d = 3 \vert p] \cdot Pr[p]}$$

  In Bayesian inference:
  $Pr[p\vert d]$: posterior probability.
  $Pr[p]$: prior probability.
  $Pr[d\vert p]$: likelihood.

• Interpretation:
  You may have certain prior data, information, or belief about the parameter. So this prior information is updated by the support from the current data, and this gives the updated posterior probability.

• For example, we have some belief or previous data showing that the coin is fair ($p=0.5$).

  

  • Our new belief about the true p is expressed in the posterior probability distribution.

  • If we want to find out the single most likely value of p, we can find the mode (or peak) of the posterior probability.
  • Or we can characterize our belief in terms of posterior mean (expectation of the posterior prob. distribution).
  • Posterior confidence in our estimate can be expressed by the variance of the posterior, 95 % quantiles, etc.
  • If you collect more data in the future, we can use the current posterior probability as the future prior probability distribution.
• When we have no information about the prior expectation about the parameters, people use flat prior probability (all parameter values are equally probable).

Connection to ABC

• When we talked about “Basic idea of ABC”, we drew random numbers for the parameters.
  • The shape of this random number distribution can be considered as the prior distribution (e.g., uniform distribution, gaussian distribution etc).
  • The histogram of accepted parameters become the posterior distribution.
• We did not have to calculate/derive the mathematical expression for the likelihood function!!!
  • This is a great advantage, because likelihood function can be only calculated for simple problems (e.g., the coin flip).
• The method described above is called Rejection method. There are more refinements in how to convert the simulation restults to posterior distribution.

Application

Let's go back to the two demography model.

1. Let's create a fake observed data:

   Let's say the divergence occured 0.7 coalescent unit ago (28000 generations), and model B was the real model.

   ms 8 1 -t 8 -I 2 3 5 -ma x 8.0 0.0 x  -n 1 0.5 -n 2 0.25  -en 0.7 1 1.0  -ej 0.7 2 1 | sample_stats | ./cleanSampleStats.pl > obs
   

   You can use the script which you created to clean up the sample_stats output. If you haven't made it, here is my version of cleanSampleStats.pl.

2. run the simulations 1 million times, and capture the output to a file (sim.out).

   nohup ./msDrive.pl > sim.out &
   

   • This could take a long time. By specifying nohup (no hung up) and running in the background, even if you log out, the simulation will keep running.
   • In LSI cluster with many users, you should use batch queing to submit your job.
   • Basically, the batch queing system will send your request to a dedicated compute node, and run it there.
   • Example:
     • Create a submit shell script (e.g. "submit.sh") with the following contents:

       #$ -N Naoki_Sim
       #$ -S /bin/bash
       #$ -cwd
       PATH=$SGE_O_PATH
       
       ./msDrive.pl > sim.out
       

       • You can change the name of your batch job in the first line (I used Naoki_sim). Don't use spaces as the part of the name. The job name is useful when you check if your job is still running.

       • For most batch jobs, you don't have to change the 2nd to 4th line.

       • Then you can type in a unix command (or multiple commands) which you want to execute.

         Notes:
         In this case, I used ./msDrive.pl. Make sure to submit the job in a directory where msDrive.pl is stored. In other words, the current working directory (cwd) where you submit the job will be the cwd of the batch job.

         Don't use ampersand (&), which indicates running a process in the background, within submit.sh.

     • Submit the job.

       $ qsub submit.sh
       Your job 18363 ("Naoki_Sim") has been submitted
       

     • Check if your job is running:

       $ qstat
       job-ID  prior   name       user         state submit/start at     queue                          slots ja-task-ID 
       -----------------------------------------------------------------------------------------------------------------
         18363 0.55500 Naoki_Sim  naoki        r     11/07/2011 17:58:06 all.q@tux01.cluster                1
       

     • Submitting a batch job creates two files (e.g. Naoki_Sim.o18363 Naoki_Sim.e18363). If there are STDOUT or STDERR outputs from your batch script, these files capture them.
     • If you noticed that you screwed up something, and stop the running job:

       $ qdel 18363
       naoki has registered the job 18363 for deletion
       

       You'll need to specify the job-ID (find it with qstat).

3. We'll use Kevin Thornton's ABCreg to create the posterior distribution. Thornton, K. R. 2009, BMC Genetics 10:35 Here is the documentation.

   reg -p sim.out -d obs -P 2 -S 5 -t 0.001 -b output -L
   

   • -P: number of parameter columns
   • -S: number of summary statistics columns
   • -t: tolerance, closest 0.1% of the simulation are used.
   • -b: base name for the output file
   • Methods;
     • -T: tangent transformation
     • -L: natural log transformation
     • no transformation if (both -T and -L are not specified)
   • For this analysis, -T and -L may not work, so you may need to run reg without these flags. For a discrete parameter (Model A vs B), linear regression used in ABCreg is not appropriate, but let's ignore the minor details for now.
4. Look at the output file (output.1.xxxx.post), xxxx can be "tangent", "log", or nothing.

   The two columns corresponds to the two parameters (model and T).

5. Exercises:
   • Write a small Perl script to check how many times model B was selected. Also, calculate the mean for the second column, which will be a estimate (posterior mean) of divergence time, T.
   • Check how well this ABC inference performs.

     To do this, you can try many different parameters to create the fake data, and check how well the estimated parameters matches with the real parameter values for each fake data.

     • You can create 100 fake observed data and put them in a file (100 lines). To do this, you can run the msDrive.pl with 100 iterations (or take 100 lines from the sim.out). Let's say this file is obsWPara. Then you create an input file (let's say this file is obsSSOnly) for reg by removing the two parameter columns (use cut).
     • If you provide 100 lines (obsSSOnly) for the -d option of reg, it will run the analyses 100 times (each line is treated as the observed data), and creates 100 output files (output.1.xxxx.post, output.2.xxxx.post, etc).
     • Now you need to run the script you created in Exercies 1 for each file to summarize the results. Then you can compare the results with the real parameter (the first two columns of obsWPara file). You can write a small script to automate this.
     • By running this kind of analyses will tell you how well your estimation (hypothesis testing) works.

6. Revised summary stats.
   • For distinguishing these two models, additional summary statistics may be useful.
   • msSS.pl

     msSS.pl -h
     

     • I wrote this script to add two additional summary statistics to sample_stats: $\pi$ (average pairwise differences) within each subpopulation, and $\pi$ between subpopulations.

     • First 5 columns are same as the cleaned output of sample_stats: pi, ss, D, thetaH, H.

     • Then the additional columns are $pi_{within,1}, pi_{within,2}, ..., pi_{within,n}$.
       After that, $pi_{bet,1-2}, pi_{bet,1-3}, ..., pi_{bet,1-n}, pi_{bet,2-3}, pi_{bet,2-4}...$
     • So in our case with only 2 subpopulations, we have 3 additional summary statistics $pi_{within,1},pi_{within,2}, pi_{bet,1-2}$ (total of 8 summary statistics)

   • Exercises:
     
     • This script is a little slow, but see if this will increase the ABC Estimation (i.e. repeat the Exerecises described in Step 5).
     • Evaluate the performances of the ABCreg analyses with -L, -T, or without -L or -T. Does one of the method work better?
     • Does using a subset of summary statistics improve the estimation?