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.

Classes and R

I heard about this Coursera course through the Davis R Users' Group (DRUG) mailing list, which is something that has just started up.  I first heard about DRUG less than a week ago, and I immediately signed up.  The first meeting is this Friday.

It is apparently not a workshop but simply two hours of the week set aside to work with R in the company of other people working with R.  At least, that seems to be the plan so far.  Hopefully this will mean that I will get a lot of analysis done, and I might even come across some interesting problems I can talk about on this blog, both in my own work and in other people's.  And I believe it will be a good environment to work with any sort of phylogenetics software.

And of course, it can be time set aside for me to work on my online course material for "Computing for Data Analysis."  I enrolled in the class just today, so I haven't had a chance yet to check out the lectures, but I plan to do so soon.  They say it is designed for first-year biostatistics graduate students, so we will see how it fits into what I already know.  Regardless, I am sure it will be a blast!  I can always use more practice.

There are lectures, quizzes, and longer programming assignments associated with the class.  Coursera seems to have quite a variety of classes available.  I hope to take a more thorough look at what relevant courses are available later . . . I would like to take any classes involving R and maybe some other languages as well.

I am also interested in checking out some of the Udacity courses, which were previously recommended to me by a friend who works for Google.  I think none of their classes are about R (mostly Python, I believe), but I am sure many are relevant to those interested in programming languages and data analysis.  They offer classes for those with no programming background as well.

Both Coursera and Udacity are free.  While some Coursera courses don't give out certificates, this one does, as long as you get a passing grade (70/100).  I believe all of the Udacity courses have an option for certification, so you can have some record that you put in the effort to learn the material.

This will be my first time taking an online course, but from my experience taking a statistical computing class, I think the assignments will be very useful in retaining information and building good habits, even without the certification accompanying it.  Even still, it could be a nice addition to a C.V., and Udacity also offers to send resumes to companies.  That's not necessarily relevant to those wishing to stay in academia, but it's nice to know that we're developing transferable skills!

likelihood reconstruction of ancestral states

In a previous post, I talked about how summarizing the state at each internal node over many make.simmap mappings did not correspond exactly with the ace reconstructions.  I contacted Liam Revell about this, and he informed me that this may be because ace does not compute the scaled marginal likelihoods but the conditional likelihoods of the subtrees descending from each node.  He suggested I try rerooting the tree at each internal node and using the ace reconstruction at the root to find the marginal likelihoods.  Here is my attempt:

nodes = (exampleTree$Nnode+2):(exampleTree$Nnode*2+1) # This gives me the number associated with each internal node

reRootAnc = t(sapply(nodes, function(x) {
  tr = reroot(exampleTree, node = x, position = 0) # rerooting the tree at each internal node
  reconst = ace(x = discrete, phy = tr, type = "discrete", model = "SYM") # estimating the maximum likelihood ancestral state estimate
  reconst$lik.anc[1,] # taking only the value for the root
}))

This is the result I got for an example tree:

Seems a little strange to me...I wonder what might be going on.

I knew this shouldn't be so complicated...

I was trying to analyze my incomplete dataset, and I needed to remove data for species where I have measured fewer than four individuals.  Since my dataset was small enough, it was easy to just remove them by hand, doing something like this:

> df
   Species meas
1        a    9
2        a    1
3        a    7
4        a    5
5        b    7
6        c    0
7        c    5
8        c    2
9        c    9
10       c    1
11       d    3
12       d    2

> dfremoved = df[-which(df$Species=="b"),]
> dfremoved = dfremoved[-which(dfremoved$Species=="d"),]
> dfremoved
   Species meas
1        a    9
2        a    1
3        a    7
4        a    5
6        c    0
7        c    5
8        c    2
9        c    9
10       c    1

But in order to do this systematically, I used a couple of steps.

> toofew = names(which(table(df$Species) < 4))
> toofew
[1] "b" "d"

First, I found the species names for species with fewer than four individuals.  With this, I can remove all rows where df$Species match any of these names.

> dfremoved = df[!(df$Species %in% toofew),]
> dfremoved
   Species meas
1        a    9
2        a    1
3        a    7
4        a    5
6        c    0
7        c    5
8        c    2
9        c    9
10       c    1