# Phylogenetics: R Primer

 EEB 349: Phylogenetics This lab represents an introduction to the R computing language, which will be useful not only for today's lab, which explores the most common probability distributions used in Bayesian phylogenetics, but also for Phylogenetics software written as R extensions (i.e. APE)

## What is R?

R is one of two general-purpose programming languages that we will use this semester (the other being Python). While R and Python share many similarities, R excels in tasks that involve statistics and creating plots. While a spreadsheet program like Excel can create charts, it is often easier to create the plot you need in R. The purpose of this lab is to introduce you to R by showing you how to make various plots related to probability distributions. This will also enable you to learn about several important probability distributions that serve as prior distributions in Bayesian phylogenetics. This introduction to R will serve us well later in the course when we make use of phylogenetic packages (e.g. APE) that are written in the R language.

## Some R basics

To use R, type commands into the console window at the R prompt (>). Try the following commands for starters. Note that I've started each line with the > character. That > character is the R command prompt: you should type in only what follows the > character. Also, you need not type in the # character or any characters that follow # on the same line. These represent comments in R, and will be ignored completely by the R interpreter.

> help()
> help(plot)            # get help on a particular topic
> help.search("gamma")  # search for all help topics that mention the word "gamma"


That last one will generate a confusing list of topics, some of which do not appear to have anything to do with the word "gamma"! Nevertheless, this search feature is good to know about if you are desperate for information about some topic for which the basic help command turns up nothing. (Google is also an excellent help search tool: start your Google search with a capital letter R followed by a keyword.)

> y = rnorm(1000, 0, 1) # generate 1000 values from a standard normal distribution (mean 0, std. dev. 1) and store in the variable y
> y                     # view the stored standard normal deviates
> hist(y)               # create a histogram of the values stored in y


Important time-saving hint: you can save a lot of typing by using your up-arrow to recall previously-issued commands. You can then modify the previously-issue command and hit the Enter key to submit the modified version.

### Creating a plot

Let's begin by making a simple scatterplot using these data:

 x y 1 1 2 1 1 3 4 5 2 4 3 5 4 7 5 8

In the R console window, enter the data and plot it as follows (note that you should not type the initial > on each line; that's the R prompt):

> x <- c(1,2,1,4,2,3,4,5)
> y <- c(1,1,3,5,4,5,7,8)
> plot(x,y)


The left-pointing arrow is the assignment operator: x <- y means assign y to x. The "c" in c(...) means "combine", so the c function's purpose is to combine its arguments. Thus, both x and y are not single values, but instead each is a collection of 8 numbers. Collections of numbers like this are known as vectors in mathematics, or arrays in many computer programming languages (although in Python it would be called a list).

### Modifying a plot

R does a fairly nice job of selecting defaults for plots, but I find myself wanting to tweak things. This section explains how to modify various aspects of plots.

#### Changing the plotting symbol

Use the pch parameter to change the symbol R uses to plot points:

> plot(x, y, pch=19)


Common values for pch include 19 (solid dot), 20 (bullet), 21 (circle), 22 (square), 23 (diamond), 24 (up-pointing triangle), and 25 (down-pointing triangle). The stroke and fill color for symbols 21-25 can be changed using the col and bg parameters. For example, to use red squares with yellow fill, try this:

> plot(x, y, pch=22, col="red", bg="yellow")


> help(points)


#### Plot type

The default plot type is "p" for points, but other plot types can be used (use the command help(plot) for details). Let's use the "l" (that's a lower-case L) type to connect our points with lines:

> plot(x, y, type="l")


Clearly, this is not something that you would want to do with these data, but we will use type="l" extensively later in this tutorial. To use both points and lines, set the type to both ("b"):

> plot(x, y, type="b")


#### Line width

To make the line connecting the points thicker, use the lwd option (lwd=1 is the default line width):

> plot(x,y,type="l",lwd=3)


#### Plot and axis labels

The xlab, ylab and main parameters can be used to change the x-axis label, the y-axis label and the main plot title, respectively:

> plot(x, y, type="l", xlab="width", ylab="height", main="Leaf width vs. height")


#### Losing the box

The bty parameter can be set to "n" to remove the box around the plot:

> plot(x, y, type="l", bty="n")


#### Losing an axis

The xaxt and yaxt parameters can be set to "n" to remove the x- and y-axis, respectively, from the plot. Here for example is a plot with no y-axis (note that I've also set the y-axis label to an empty string with ylab="" and removed the box around the plot with bty):

> plot(x, y, type="l", yaxt="n",ylab="",bty="n")


#### Changing axis labeling

It is often desirable to change the values used on the x- or y-axis. You can use the ylim parameter to make the y-axis extend from 0 to 10 (instead of the default, which is 1 to 8 in this case because that is the observed range of y values):

> plot(x, y, type="l", ylim=c(0,10))


Note that there is a little overhang at both ends on the y-axis (the 0 tick mark is slightly above the bottom of the plot and the tick mark for 10 is slightly below the top of the plot). This is because, by default, the y-axis style parameter yaxs is set to "r", which causes R to extend the range you defined (using ylim) by 4% at each end. You can set yaxs="i" to make R strictly honor the limits you set:

> plot(x, y, type="l", ylim=c(0,10), yaxs="i")


You can change the x-axis in the same way using xlim and xaxs instead of ylim and yaxs.

### Data frames

The c(...) notation used above is fine for small examples, but if you have a lot of data, you will prefer to read in the data from a file. R uses the concept of a data frame to store a collection of data. For example, save the following data in a text file on your hard drive named myfile.txt:

    one two three four
A    0   1    2     3
B    4   5    6     7
C    8   9   10    11
D   12  13   14    15


Note that the first row has only 4 elements whereas all others have 5. If the first row contains one fewer element than the other rows, then R assumes that the first row is a list of variable names.

Use this command to read these data into a data frame:

> d <- read.table("myfile.txt")


If you get an error, it is probably because the current working directory used by R is not the directory in which you stored your file. You can ask R to tell you the current working directory:

> getwd()


and you can set the current working directory as follows (here I'm assuming that you are using Windows and myfile.txt is stored in C:\phylogenetics\rlab:

> setwd("C:\\phylogenetics\\rlab")


Note that (for rather esoteric reasons) the single backslash characters in the actual path must be replaced by double backslashes in R commands.

If instead you were using a Mac and myfile.txt was stored on your desktop, then you would do this:

> setwd("~/Desktop")


To ensure that your data were read correctly, simply type d to see the contents of the variable you created to hold the data:

> d
one two three four
A   0   1     2    3
B   4   5     6    7
C   8   9    10   11
D  12  13    14   15


To extract one column of your stored data frame, use a $followed by the column name: > col1 <- d$one
> col1
[1]  0  4  8 12


You can also extract single values, although this is not commonly needed:

> x <- d[2,3]
> x


[1] 6

## Exploring probability distributions using R

R's plotting abilities and its built-in knowledge of many probability distributions used in statistics make it useful for exploring these distributions. The distributions described here are used extensively as prior distributions for parameters in Bayesian phylogenetic analyses, so it is worthwhile to familiarize yourself with their basic properties. With R we can examine them graphically.

### Gamma distribution

The gamma distribution has two parameters that govern its shape (parameter $\alpha$) and scale (parameter $\beta$). Gamma distributions range from 0.0 to infinity, so they provide excellent prior distributions for parameters such as relative rates (e.g. kappa, GTR relative rates, omega) and branch lengths that can be arbitrarily large but not negative.

Here are some basic facts about gamma distribuitons:

 Shape parameter $\alpha$ Scale parameter $\beta$ Mean $\alpha \; \beta$ Variance $\alpha \; \beta^2$ Density function $f(x) = \frac{x^{\alpha - 1} e^{-x/\beta}}{\beta^{\alpha} \Gamma(\alpha)}$

Note that that $\Gamma(\alpha)$ in the density function refers to the gamma function (which is distinct from the gamma distribution). The gamma function is easy to calculate when its argument is a positive integer:

$\Gamma(x) = (x-1)!$

When, however, the argument is an arbitrary positive real number (which is the usual situation for $\alpha$), then computing it gets more complicated. Fortunately, R handles details like this for us!

When the gamma distribution is used to model rate heterogeneity in phylogenetics, its mean is set to 1 because it is being used to model relative rates. In this application, only the shape parameter $\alpha$ is ever specified, and that is because the scale parameter $\beta$ is determined by the constraint that the mean is 1.0.

What is the value of $\beta$ if $\alpha$ is known and the mean must be 1.0?

#### Creating a Histogram

Let's begin by drawing 1000 random numbers (deviates) from a gamma distribution (shape=0.5, scale=2.0) and creating a histogram. Type the following command into the R console to generate the random deviates and store them in the variable x:

> x <- rgamma(1000, shape=0.5, scale=2.0)


The "r" in "rgamma" stands for "random". There is a corresponding "r" command for all the common probability distributions (e.g. rbeta, rbinom, rchisq, rexp, rnorm, runif) and in each case it allows you to draw a sample of values from that distribution (e.g. the beta, binomial, chi-square, exponential, normal and uniform distributions, respectively).

It may appear that nothing has happened, but you can view the 1000 values stored in the variable x by typing x at the prompt.

Creating a basic histogram from these 1000 values is easy:

> hist(x)


Now refine the histogram by asking R to give you 50 bars:

> hist(x, breaks=50)


Note that breaks represents a suggestion only. R will often not give you exactly as many bars as you want.

Create a histogram with approximately 40 bars from 10000 deviates from a Gamma(shape=10, scale=0.1) distribution
Create a histogram with 30 bars from 10000 deviates from a gamma distribution in which the mean is 2 and the variance is 1
Assuming you have stored the 10000 deviates in the variable x, type summary(x) to get the sample mean, median, etc. How close is the sample mean from the expected mean? Draw another sample with 1 million deviates and summarize. How close is the sample mean to the expected mean now?

#### Plotting the density function

The "d" commands (e.g. dgamma, dbeta, dexp, dnorm, dunif) can be used to compute the probability density function for a series of values. In this part of the tutorial, you will use the dgamma function to compute the gamma density function for all 501 values in the series 0.0, 0.01, 0.02, ..., 5.0. First, generate the values in the series and store in a variable named x:

> x <- seq(0.0, 5.0, 0.01)


The seq generates a series of values starting at the first supplied value (0.0 in this case), stopping at the second supplied value (5.0 in this case), with spacing between values in the sequence determined by the third supplied value (0.01 in this case). To see if 501 values were generated, type

> length(x)


You can list all the values in x by simply typing x at the console.

R makes it easy to compute a quantity for every number in a list. Here is how to compute the gamma density function value at each value stored in the variable x:

> y <- dgamma(x, shape=10, scale=0.1)


Before plotting the density function, you should convince yourself that R is computing the density correctly by calculating one value yourself by hand. First print out the 100th x value as follows:

> x[99]


(Note that indexes begin at 0 in R (as in most programming languages), so the first value is x[0], the second is x[1], and so forth. That's why the 100th value has index 99.) The value displayed should be 0.98.

Calculate (using the calculator application on your computer) the gamma density for the value 0.98 using $\alpha$=10 and $\beta$ = 0.1

Compare your value to the value computed by R:

> y[99]


Now plot the density function as follows:

> plot(x,y)


This command says to plot all possible pairs of points, where one value is taken from x and the other from y. There are 501 values in both x and y, so 501 points will be plotted. This plot looks kindof ugly with lots of overlapping points. More reasonable would be to connect the points with a line but not to print all 501 points themselves. You can do this by adding type="l" to your plot command:

> plot(x,y,type="l")


The type is specified to be a lower case L character here (L for line). For other plot types, see the online help by typing:

> help(plot)


#### Using the distribution function

We've seen the "d" and "r" commands (for density function and random number generation), but there are two additional standard commands defined for probability distributions. The "p" command computes the cumulative distribution function for a distribution. For the gamma distribution, you can use the pgamma command to compute the integral of the density function up to a specified value. First, plot a gamma density with shape 2 and scale 0.5:

> x <- seq(0.0, 5.0, 0.01)
> y <- dgamma(x, shape=2.0, scale=0.5)
> plot(x, y, type="l")

The total area under this curve is 1.0. What is your guess as to the area of the curve from 0 to 1?

Compute the area under the density from 0 to 1 as follows using the pgamma command:

> pgamma(1, shape=2, scale=0.5)
[1] 0.5939942


Thus, nearly 60% of random draws from a Gamma(2,0.5) distribution would be less than 1.0. You can compute several cumulative probabilities by supplying a vector to the pgamma command. Suppose you wanted to know the cumulative probability for each of the values on the x-axis (i.e. 1, 2, 3, 4 and 5):

> pgamma(c(1,2,3,4,5), shape=2, scale=0.5)
[1] 0.5939942 0.9084218 0.9826487 0.9969808 0.9995006


Note that 5 values are produced by this command, corresponding to the 5 values supplied in the vector c(1,2,3,4,5). The last one (0.9995006) means that 99.95% of random draws from this distribution would be less than 5.0.

Use the pgamma command to compute the probability that a value drawn from a Gamma(2,0.5) distribution would fall between 0.5 and 4.5? Hint: you will need to either invoke the command twice or provide it with a vector containing two values, then do some subtraction with the results.

#### Quantiles

A quantile is a point along the x-axis of a density plot that corresponds to a particular cumulative probability. For example, the 60% quantile corresponds to a cumulative probability of 0.6. If you computed the area under the density curve up to some point x and the cumulative probability was 0.6, then x is the 60% quantile. For example, in the previous section, the value 1 is very close to the 60% quantile. Suppose we wanted to divide up a Gamma(2, 0.5) distribution into four equal-area pieces, and needed to find the points along the x-axis corresponding to the boundaries between these pieces (does this problem sound familiar?). We could do this in R as follows:

> qgamma(c(0.25, 0.5, 0.75), shape=2, scale=0.5)
[1] 0.4806394 0.8391735 1.3463173


Note that we have supplied three values to the qgamma function via the vector c(0.25, 0.5, 0.7). We did not need to supply the fourth value (1.0) because we know that the 100% quantile is equal to infinity. In case you haven't guessed, these would be the boundaries used if we applied discrete gamma rate heterogeneity to a model, specifying 4 categories and shape=2.0. You can check this in PAUP if you like using PAUP's gammaplot command (note that PAUP requires you to first execute a data file):

gammaplot shape=2 ncat=4;


To save you the trouble of logging into the cluster, here is the table PAUP produces:

Cut-points and category rates for discrete gamma approximation

          ------ cut-points ------
category      lower        upper      rate (mean)
-------------------------------------------------
1      0.00000000   0.48063938     0.29327472
2      0.48063938   0.83917350     0.65501368
3      0.83917350   1.34631726     1.06998966
4      1.34631726     infinity     1.98172195


So the pgamma command can be used to ask questions like "What are the chances that a Gamma(2,0.5) deviate would be less than 1.0?" whereas the qgamma command is used to ask the related question "What interval accounts for 95% of the values drawn from a Gamma(2,0.5) distribution?"

What value x produces an interval (0,x) that accounts for 95% of the values drawn from a Gamma(2,0.5) distribution?
What values x1,x2 produces an interval (x1,x2) that accounts for 95% of the values drawn from a Gamma(2,0.5) distribution (2.5% of the distribution should be less than x1 and 2.5% greater than x2)?

### Exponential distribution

The exponential distribution is a special case of the gamma distribution. Exponential distributions are gamma distributions in which the shape parameter is 1.0. The single parameter of an exponential distribution is called the rate parameter and is the inverse of the shape parameter of the corresponding gamma distribution. Exponential distributions range from 0.0 to infinity, so they provide excellent prior distributions for all the parameters for which gamma distributions are appropriate (various relative rates and branch lengths). Exponential distributions are popular as priors because they force you to only come up with one number (the rate parameter) rather than the two (shape and scale) needed for gamma distributions, but this does not in any way justify using Exponential distributions over Gamma distributions!

Here are some basic facts about gamma distributions:

 rate parameter $\lambda$ Mean $\textstyle \frac{1}{\lambda}$ Variance $\textstyle \left(\frac{1}{\lambda}\right)^2$ Density function $f(x) = \lambda e^{-\lambda x}$

#### Exponential densities

Let's plot a couple of exponential densities and compare them. First, plot an Exponential(10) density function as follows:

> x <- seq(0.0,5.0,0.01)
> y <- dexp(x, rate=10)
> plot(x, y, type="l", col="blue", lwd=3, ylim=c(0, 10), xaxs="i")


Note that I've colored the line blue (col="blue"), made it thick (lwd=3), set the y-axis limits to 0.0 and 10.0 (ylim=c(0,10)) and specified that the x-axis limits should be adhered to srictly (no added 4% to the left and right ends).

Now let's add an Exponential(0.1) density (in red) to the plot using the lines function:

> y2 <- dexp(x, rate=0.1)
> lines(x, y2, col="red", lwd=3, ylim=c(0, 10), xaxs="i")


The lines command is like the plot command except that it adds a line to an existing plot rather than creating a new plot (also, notice that type="l" is not necessary). (As you might imagine, there is also a points command that adds points to an existing plot.)

Is the height of the density at x=0 related to the rate parameter?
Which of these two densities has the larger mean? (Try to answer this question first by looking at the plot, then compute the numbers to check your intuition.)
Which of these two densities has the larger variance? (Again, answer the question first by looking at the plot.)
What fraction of the Exponential(10) distribution is less than 1? (Hint: use the pexp command.)
Fill in the blank: 95% of the Exponential(0.1) distribution is less than _____. (Hint: use the qexp command.)

### Beta distribution

The Beta distribution differs from the Gamma distribution in that Beta random variables have an upper bound (1.0) as well as a lower bound (0.0). Two parameters (a and b, but often referred to as $\alpha$ and $\beta$) govern the shape of the Beta distribution. Beta distributions are natural priors for parameters such as pinvar (and the proportion of heads parameter in coin flipping examples), which are restricted to the [0.0, 1.0] interval.

Here are some basic facts about beta distributions:

 shape parameter 1 $a$ shape parameter 2 $b$ Mean $\textstyle \frac{a}{a + b}$ Variance $\textstyle \frac{a b}{(a + b)^2 (a + b + 1)}$ Density function $f(x) = \frac{\Gamma(a + b)}{\Gamma(a) \Gamma(b)} x^{a - 1} (1-x)^{b-1}$

The $\Gamma$ function above is the same one used in the Gamma distribution.

#### Symmetrical Beta distributions

Beta densities are symmetrical about 0.5 if $a=b$. Try plotting a Beta(3,3) distribution in R as follows:

> x <- seq(0.0, 1.0, 0.01)
> y <- dbeta(x, 3, 3)
> plot(x, y, type="l")


As $a$ and $b$ get larger (still constraining a to equal b), the density becomes more and more sharply peaked at 0.5. Regenerate y based on a Beta(1000,1000) density and plot it to see this:

> y <- dbeta(x, 1000,1000)
> plot(x, y, type="l")


Note that if $a=b$, the variance formula given above simplifies to just $1/(8a + 4)$, so the variance of a Beta(3,3) distribution is $1/28 = 0.0357$, whereas the variance of a Beta(1000,1000) distribution is $1/8004 = 0.000125$, which is 286 times smaller.

If a=b=1, the Beta distribution is identical to a Uniform(0,1) distribution, and the density function is simply a straight line at height 1.0:

> y <- dbeta(x, 1,1)
> plot(x, y, type="l")

Try plotting a Beta density when $a=b=0.1$

This distribution might be useful for a coin-flipping analysis if you believed strongly that the world was full of trick coins (two-headed or two-tailed coins)!

#### Asymmetrical Beta distributions

If $a \neq b$, Beta densities are skewed to the right or left and are not symmetrical. Take a look at the formula for the mean: $a/(a+b)$. If $a < b$, the mean is less than 0.5 and the Beta distribution leans to the left. If $a > b$, the Beta distribution leans to the right.

Try plotting these Beta densities: Beta(1,100), Beta(100,1), Beta(2,5), Beta(5,2)
Knowing what Beta(2,5) is like, try to predict what a Beta(20,50) density would look like, then plot it to check your intuition

Calculate the mean and variance of a Beta(2,5) distribution using the formulas provided at the beginning of this section. Now sample 10000 values from a Beta(2,5) distribution and compute the sample mean and variance as follows:

> x <- rbeta(10000, 2, 5)
> mean(x)
[1] 0.2850546
> var(x)
[1] 0.02490558


Create a histogram of the values in x:

> hist(x)


Save the histogram as a PDF file named "beta25a.pdf" (we'll be using this file again in a minute). On a Mac, you can create the PDF file by clicking on the window containing the histogram, then choosing File, Print... to bring up the "Print" dialog box, and finally choosing "Save as PDF..." from the PDF drop-down list. In Windows, click on the window containing the histogram, then choose File, Save As, PDF... to bring up the "PDF file" dialog box, where you can specify the file name.

#### Relationship to the Gamma distribution

Beta distributions are related to Gamma distributions. If $x$ is a Gamma($a$,1) random variable, and $y$ is a Gamma($b$,1) random variable, then $x/(x+y)$ is distributed as Beta($a$,$b$). Now that you have just generated 10000 Beta(2,5) random variables, let's compare that to 10000 random variables generated as described above:

> x <- rgamma(10000, 2, 1)
> y <- rgamma(10000, 5, 1)
> z <- x/(x+y)
> mean(z)
> var(z)


Note in line 3 that R automatically applies formulas repeatedly to each value of a list. Both x and y are lists of 10000 values. Line 3 produces 10000 values stored in z where z[i] = x[i]/(x[i] + y[i]) for i = 1, 2, ..., 10000.

Create a histogram of the values in z, and save as a PDF file named "beta25b.pdf" (the procedure is described in the previous section). Compare "beta25a.pdf" with "beta25b.pdf" (they should be very similar because we have simply generated 10000 Beta(2,5) random variables in two different ways).

### Dirichlet distribution

Gamma distributions can be used as priors for branch lengths, relative rates, the among-site rate heterogeneity shape parameter and other parameters whose range is 0 to infinity. Beta distributions are good for proportions, such as pinvar, and parameters whose range is 0 to 1. What about base frequencies? The four base frequencies all lie between 0 and 1, but any prior distribution used for base frequencies must take account of the fact that the four base frequencies must add up to 1. The Dirichlet distribution is perfect for constrained proportions like this.

The Dirichlet behaves in many ways like a Beta distribution. In fact, the Beta distribution is a special case of the Dirichlet distribution or, said another way, the Dirichlet distribution is a multivariate version of the Beta distribution. Think of the Beta distribution as a distribution governing two proportions, p and 1-p, that are also subject to the same constraint as base frequencies (that they must add to 1). A Beta(a,b) distribution is thus equivalent to a Dirichlet(a,b) distribution. There are four base frequency parameters, so four Dirichlet parameters are needed (i.e., a Dirichlet(a,b,c,d) distribution) to create a prior appropriate for a vector of four base frequencies: $\pi_A, \pi_C, \pi_G, \pi_T$.

Here are some basic facts about the Dirichlet(a,b,c,d) distribution, specifically when used as a prior distribution for base frequencies:

 shape parameters 1, 2, 3, and 4, respectively $a$, $b$, $c$, $d$ sum of shape parameters (used in formulas below) $s = a + b + c + d$ Mean of $\pi_A$, $\pi_C$, $\pi_G$, $\pi_T$, respectively $\textstyle \frac{a}{s}, \frac{b}{s}, \frac{c}{s}, \frac{d}{s}$ Variance of $\pi_A$, $\pi_C$, $\pi_G$, $\pi_T$, respectively $\textstyle \frac{a (s - a)}{s^2 (s+1)}, \frac{b (s - b)}{s^2 (s+1)},\frac{c (s - c)}{s^2 (s+1)},\frac{d (s - d)}{s^2 (s+1)}$ Density function $f(\pi_A, \pi_C, \pi_G, \pi_T) = \frac{\Gamma(s)}{\Gamma(a) \Gamma(b) \Gamma(c) \Gamma(d)} \pi_A^{a - 1} \pi_C^{b - 1} \pi_G^{c - 1} \pi_T^{d - 1}$

#### Relationship to the Gamma distribution

Just like Beta distributions, Dirichlet distributions are related to the Gamma distribution. If

$x_1 \sim \mbox{Gamma}(a,1)$

$x_2 \sim \mbox{Gamma}(b,1)$

$x_3 \sim \mbox{Gamma}(c,1)$

$x_4 \sim \mbox{Gamma}(d,1)$

(where $x_1 \sim \mbox{Gamma}(a, 1)$ is shorthand for "$x_1$ is a random number drawn from a Gamma distribution with shape $a$ and scale 1"), then a Dirichlet(a,b,c,d) random variable (a vector of 4 values) can be created by dividing each $x_i$ by the sum of all four $x_i$ values. Create 10 random Dirichlet(1,1,1,1) variables using this method and list them:

> x1 <- rgamma(10, 1, 1)
> x2 <- rgamma(10, 1, 1)
> x3 <- rgamma(10, 1, 1)
> x4 <- rgamma(10, 1, 1)
> s <- x1 + x2 + x3 + x4
> d <- c(x1/s, x2/s, x3/s, x4/s)


Note the use of the c(...) function to construct a vector of four values, each of which is one of the $x_i$ values divided by $s$, which is the sum of all four $x_i$ values.

Now display the 10 Dirichlet vectors:

> d
[1] 0.09234290 0.37789542 0.08640752 0.28049539 0.14485492 0.05511074 0.14355295
[8] 0.05358819 0.03350763 0.19362171 0.39036151 0.25134264 0.29663099 0.43909681
[15] 0.13981216 0.29531499 0.34553287 0.76487943 0.13311580 0.23832980 0.11153418
[22] 0.35743685 0.03021077 0.12139093 0.63609473 0.17125022 0.27907977 0.09122259
[29] 0.44937431 0.19609418 0.40576141 0.01332508 0.58675072 0.15901687 0.07923820
[36] 0.47832406 0.23183441 0.09030979 0.38400226 0.37195431


Unfortunately, R dumped all 10 vectors together, creating a single vector of 40 values! We can coerce R into displaying these 4 at a time by using the dim command to redimension d into a matrix with 10 rows and 4 columns:

> dim(d) <- c(10,4)
> d
[,1]      [,2]       [,3]       [,4]
[1,] 0.09234290 0.3903615 0.11153418 0.40576141
[2,] 0.37789542 0.2513426 0.35743685 0.01332508
[3,] 0.08640752 0.2966310 0.03021077 0.58675072
[4,] 0.28049539 0.4390968 0.12139093 0.15901687
[5,] 0.14485492 0.1398122 0.63609473 0.07923820
[6,] 0.05511074 0.2953150 0.17125022 0.47832406
[7,] 0.14355295 0.3455329 0.27907977 0.23183441
[8,] 0.05358819 0.7648794 0.09122259 0.09030979
[9,] 0.03350763 0.1331158 0.44937431 0.38400226
[10,] 0.19362171 0.2383298 0.19609418 0.37195431


Each row now represents a random draw from a Dirichlet(1,1,1,1) distribution. If we were using this Dirichlet distribution as a prior for base frequencies, each row would represent one representative vector of base frequencies from our prior. Given the impossibility of plotting a 4-dimensional density function, it would thus be useful to use the method above to produce samples from various Dirichlet distributions to see what they imply about the distribution of base frequencies.

To make it easier to explore the properties of Dirichlet distributions, let's put the commands you just entered into a function that can be called with different combinations of parameters. Copy the following into your clipboard and paste it into the R console window.

Dirichlet <- function(a,b,c,d) {
x1 <- rgamma(10, a, 1)
x2 <- rgamma(10, b, 1)
x3 <- rgamma(10, c, 1)
x4 <- rgamma(10, d, 1)
s <- x1 + x2 + x3 + x4
d <- c(x1/s, x2/s, x3/s, x4/s)
dim(d) <- c(10,4)
d
}


#### Symmetrical Dirichlet(a,b,c,d) distributions

Like Beta distributions, Dirichlet distributions are symmetrical if all parameters are equal. Use your new R function, which I've arbitrarily (but appropriately!) named "Dirichlet", to explore a series of symmetrical Dirichlet distributions:

> Dirichlet(1,1,1,1)
> Dirichlet(100,100,100,100)
> Dirichlet(1000000,1000000,1000000,1000000) # yes, that's 1 million!

Can you see a trend? What is the effect of increasing the magnitude of all four parameters from 1 to 100 to 1 million?
Which of these three Dirichlet distributions comes closest to fitting your pre-conceived notions about base frequencies?

In the case of a Dirichlet(a,b,c,d) distribution, as the parameters become larger (assuming $a=b=c=d$), the density becomes more and more sharply peaked at (0.25, 0.25, 0.25, 0.25). The magnitude of the Dirichlet parameters thus determines the informativeness of the prior: if all four parameters equal 1, the Dirichlet is flat or uninformative (any possible combination of base frequencies has the same probability density as any other combination), whereas if all four parameters equal 1 million, then approximately equal base frequencies (e.g. 0.250241, 0.250133, 0.249963, 0.249663) get much more weight than very unequal base frequencies (e.g., 0.000404, 0.531065, 0.000002, 0.468529) and the distribution would be considered informative.

We can visualize these samples better by plotting them. Below is an R function that plots the four base frequencies in different colors (note that it draws 1000 samples rather than just 10):

PlotDirichlet <- function(a,b,c,d) {
x1 <- rgamma(1000, a, 1)
x2 <- rgamma(1000, b, 1)
x3 <- rgamma(1000, c, 1)
x4 <- rgamma(1000, d, 1)
s <- x1 + x2 + x3 + x4
d <- c(x1/s, x2/s, x3/s, x4/s)
dim(d) <- c(1000,4)
plot(1:1000, d[,1], type="l", col="red", ylim=c(0,1), ylab="frequencies")
lines(1:1000, d[,2],col="blue")
lines(1:1000, d[,3],col="green")
lines(1:1000, d[,4],col="black")
}


The 1:1000 in the plot and lines commands is shorthand for seq(1, 1000, 1) (i.e., a sequence of values starting at 1, ending at 1000, with a step size of 1. Also, d[,1] means "the first column of d", d[,2] means "the second column of d, etc.

Now plot a few Dirichlet distributions:

> PlotDirichlet(1,1,1,1)
> PlotDirichlet(10,10,10,10)
> PlotDirichlet(100,100,100,100)


The x-axis in these plots doesn't have any meaning (it is just the row number of the matrix d), but looking at the plot it is easy to see which Dirichlet distributions have higher variance (and hence are less informative) and which have lower variance (and thus are more informative).

#### Asymmetrical Dirichlet distributions

If some of the four parameters differ from the others, Dirichlet(a,b,c,d) densities are not symmetrical. It might be useful to consider using asymmetrical Dirichlet distributions if you know your sequences are GC rich, for example. Here are three GC-rich Dirichlet distributions varying in informativeness but all having mean $(\pi_A, \pi_C, \pi_G, \pi_T) = (0.2, 0.3, 0.3, 0.2)$:

> PlotDirichlet(2, 3, 3, 2)
> PlotDirichlet(20, 30, 30, 20)
> PlotDirichlet(200, 300, 300, 200)


#### Dirichlet prior for GTR relative rates

As you know, the GTR model specifies six relative rates corresponding to the six possible types of substitutions in the rate matrix (actually there are 12, but the GTR model is symmetrical in that the relative rate for A to G substitutions is identical to that for G to A substitutions). Bayesian programs such as MrBayes typically use a 6-parameter Dirichlet(a,b,c,d,e,f) distribution as the prior distribution of these GTR relative rates.

Below is a revised Dirichlet function that takes 6 parameters rather than 4. Copy this into your R console:

Dirichlet <- function(a,b,c,d,e,f) {
x1 <- rgamma(10, a, 1)
x2 <- rgamma(10, b, 1)
x3 <- rgamma(10, c, 1)
x4 <- rgamma(10, d, 1)
x5 <- rgamma(10, e, 1)
x6 <- rgamma(10, f, 1)
s <- x1 + x2 + x3 + x4 + x5 + x6
d <- c(x1/s, x2/s, x3/s, x4/s, x5/s, x6/s)
dim(d) <- c(10,6)
d
}


Now apply the new function:

> Dirichlet(1,1,1,1,1,1)
> Dirichlet(1000,4000,1000,1000,4000,1000)


Now each row represents a set of GTR relative rates. The first case is a flat prior: it says that we have no idea which of six substitution classes we expect to be evolving at a relatively fast rate compared to the other substitution classes. The second case is an informative prior that says we expect the transition rate to be 4 times that of the transversion rate (assuming the ordering is $r_{AC}, r_{AG}, r_{AT}, r_{CG}, r_{CT}, r_{GT}$).

If you are not completely exhausted, and want a challenge, try creating a function to plot 6-parameter Dirichlet distributions in the same way we plotted 4-parameter Dirichlet distributions above

## Conclusion

This concludes the introduction to the R programming language and its use in exploring various probability distributions used as priors in Bayesian phylogenetics. Next week you will see how this knowledge is used in running analyses with the widely-used Bayesian phylogenetics software program MrBayes.