Fun with vectorizing

I found myself staying up way too late and having way too much fun helping a friend with her R homework. One of the things we did was to vectorize a function! Her assignment was to create plots using the Ricker model, with varying values for r:

So the first thing we did was to write a function that would calculate the population size over a certain time period for a single value of r:

Ricker = function(r, K = 100, n0 = 50, tf = 100) {
    n = numeric(length = tf)
    n[1] = n0
    for(t in 1:(tf-1)) {
          # We have tf - 1 because when t = 1, we are actually
          # calculating n[2], so if we want to end on n[tf], t 
          # should end on tf - 1
         n[t+1] = n[t]*exp(r*(1-n[t]/K))
    }
    n
}

Changing this to work with a longer vector of r values was fairly simple:

vectorizedRicker = function(r, K = 100, n0 = 50, tf = 100) {
     NumofRs = length(r)
      # First, we're changing our output to be
      # a matrix instead of a vector
     n = matrix(0, NumbofRs, tf)
      # We could vectorize over n0 as well, or any of the 
      # other parameters, if we set the number of rows for n
      # to be equal to the length of the longest vector
      # of parameter inputs
     n[,1] = n0
     for(t in 1:(tf-1)) {
           # Now, we just change the for loop to reflect that it
           # is working on entire columns of a matrix rather
           # than just a single element of a vector
          n[,t+1] = n[,t]*exp(r*(1 - n[,t]/K))
     }
     n
}

rfishbase

I was given a task to use the data in Fishbase to find percentages of different fish in different habitats. Carl Boettiger created an R package (rfishbase) to make the data very accessible in R, but unfortunately, some of the functions in the version on CRAN didn't work for me. Instead, I used the code he has on github. So if anyone else is having trouble getting updateCache() to work or wants to modify the code to extract a piece of information that isn't being extracted in the original fishbase() function, this is the place to go! Another nice thing is that you can also use getData() to download smaller chunks of data if you don't want to overwhelm the server.

Moving on from base

So I gave myself a project in my last post to remove hard-coded numbers in my code to plot a stacked histogram. I succeeded, but it may have been more work than it's worth. Perhaps it's time for me to move on to ggplot2 or lattice?

cats = levels(dat$Type) # same as before

xaxis = round(range(scores$x)*5)/5 # the range, rounded to every .2

breaks = seq(xaxis[1]-.2, xaxis[2]+.2, .2) # the breaks to use for histograms

histData = sapply(cats, function(x) {
     hist(scores$x[which(dat$Type == x)],
          breaks = breaks)$counts
})
allhistData = do.call(rbind, list(histData))
barplot(t(allhistData), space = 0, ylab = "number of sequences", xlab = "LD")
zero = which(breaks==0)-2 # find the 0 on my x-axis. The plot starts at -1

marks = seq(floor(xaxis[1]), ceiling(xaxis[2])) # the tick-marks I want to use, every integer

ticks = numeric() # this vector will hold where I want the tick labels to appear

while(zero >= -1) {
  ticks[length(ticks)+1] = zero
  zero = zero - 5
} # goes down by five, which corresponds to one LD unit

while(length(ticks) < length(marks)) {
  ticks[length(ticks)+1] = max(ticks) + 5
} # go up by five, until I have as many locations as I have labels

axis(side = 1, at = ticks[order(ticks)], labels = marks) #I need to order the tick marks in sequential order

legend("topleft", legend = cats, text.col = c("black", "gray"))

This produces the same plot from the previous post.

Linear Discriminant Analysis

This follows fairly naturally from the PCA I did on the data in my previous post. I have a dataset of several quantitative variables that can be grouped by a categorical variable. This time, I am going to maximize the separation between the two groups to see what traits are important in determining group inclusion.

While looking for information about linear discriminant analyses (LDA), I came across a very informative site by Dr. Avril Coghlan. I used a couple of the functions available from that website, including groupStandardise and calcWithinGroupsVariance, to obtain meaningful coefficients to determine what traits are more informative about group membership.

The lda() function is in the MASS package. Since I only have two groups, I only have one discriminant axis, as the number of discriminant axes is equal to the number of groups minus one. Thus I won't get a nice scatterplot the way I did with my PCA. Instead, if I try to plot the output to lda(), I get this:

I wanted a stacked histogram instead, so I had to do a little more work. 

cats = levels(dat$Type) # the categories I'm using

histData = sapply(cats, function(x) {

  hist(scores$x[which(dat$Type == x)],

       breaks = seq(-3, 2.2, .2))$counts

# bad me, I hard-coded the breaks. My next project can be to use the same breaks as the above plot

})

# this gives me a matrix of the counts for each interval for each type

allhistData = do.call(rbind, list(histData))

barplot(t(allhistData),

        space = 0,

        ylab = "number of sequences",

        xlab = "LD")

axis(side = 1,

     at = c(-1, 4, 9, 14, 19, 24),

     labels = c(-3, -2, -1, 0, 1, 2)) # more bad hard-coding

legend("topleft", legend = cats, text.col = c("black", "gray"))

And here is what I get:

It doesn't look great, but when did stacked histograms ever look good? Probably better if I have a lot more observations.

PCAs and Plotting

Principal Components Analysis, or PCA, is fairly straightforward using the princomp() function in R. But the data I have is divided by two factors: the type of response and the individual. I wanted to plot using different colors for the type of response and different plotting symbols for the individuals. Simple enough, but I also didn't want to need to modify the code if the number of types of responses or the number of individuals changed (the latter is probably more likely, but I still want my code to be as general as possible). Here is an example generated from random data:

The points themselves are fairly simple. I simply referred to the two variables I wanted in my col and pch arguments of plot().

col = as.numeric(dat[,x])+2

and

pch = dat$Individual+14

The reason for the +2 and +14 are to get the colors and plotting symbols I wanted. I could also assign specific colors I want by doing this:

colors = c("green", "blue")
ptype = c(15, 16, 17)

Then add these as arguments to plot().

col = colors[as.numeric(dat[,x])
pch = ptype[dat$Individual]

In this case, I would need to make sure that I have enough colors and plotting symbols that I don't run out.

The fun happens with the legend. For the text, I used this argument:

legend = c(levels(dat[,x]), unique(dat$Individual))

This remains flexible for any number of types or individuals. For the colors, I used a combination of seq() and rep() to get the numbers I wanted. If I didn't want my code to be general, I could simply use this:

text.col = c(3, 4, 1, 1, 1)

Instead, I used this:

text.col = c(seq(3, length(levels(dat[,x]))+2), rep(1, length(unique(dat$Individual))))

seq(3, length(levels(dat[,x]))+2) gives me a sequence of integers from three all the way to the number of types I have plus two (because I started with three instead of one). rep(1, length(unique(dat$Individual))) gives me 1 repeated for every unique individual.

Finally, the plotting symbols:

pch = c(rep(NA, length(levels(dat[,x]))), unique(dat$Individual)+14)

This is essentially the opposite of what I just did for the text color, except that I don't want any symbol next to the two types. Even though this looks like (and is) a lot more typing than simply hard-coding the appropriate numbers, this lets me use exactly the same code to make the figure even after I have doubled or tripled the amount of data I have.

jModelTest2

I had been having difficulties with jModelTest, and while searching for information about one of the problems I was having**, I discovered jModelTest2. I am quite surprised I hadn't come across it before, since it apparently has been around in some form or another for over a year, but I somehow failed to hear about it until now.

The first thing of note is that jModelTest2 can test all submodels of the GTR model. I am attempting to do this with my first trial to see how long it takes. Unlike jModelTest, jModelTest2 can use a hill-climbing algorithm to systematically find the best-fit models. Because the software I use (Phycas, MrBayes) doesn't accomodate all submodels, this may not be too relevant. However, BEAST or RevBayes can accomodate any submodel, so it may be useful in conjunction with those.

The second thing of note is high performance computing, supporting up to eight threads. Unfortunately, I don't really know much about this. My computer has a single-core processor, but it still seems to benefit from the threading. I may do a test run to see if the time it takes is dramatically reduced. Also not quite sure if "iddle" is a typo.

All-in-all, it seems like jModelTest2 is having no difficulties with the files that jModelTest threw errors for, so I take that as a good sign.

** In the old version of jModelTest, I was getting this error:

Notice "BIONJ-JC tree: null." This was the error I was trying to find more information on when I found out about jModelTest2.

RColorBrewer

Another great package that I was first introduced to in my Statistical Computing class, RColorBrewer.

This is a package that has built-in palettes that allows you to choose colors that have enough contrast for making plots.

The palettes come in three basic types: sequential, diverging, and qualitative. The sequential and diverging are great for plots where you want the colors to show an order. The difference between sequential and diverging seems to be a little subtle in terms of need: sequential shows more of a gradient, while diverging emphasizes both high and low extremes.

The qualitative palettes are best for categorical data with no ordering among categories. There are many sets, but they differ in the number of maximum colors, from 8 to 12. One interesting palette is the Paired palette, which consists of 6 hues, each with a light and dark color. I used this recently in a plot of different species, with males and females of varying lightness.

The same effect can be produced, perhaps to better effect, with different plotting symbols, but the Paired palette does a pretty good job.