Phylogenetically Independent Contrasts

##############################################################
# ============================================================
# Lab 5: Phylogenetically Independent Contrasts using APE in R
# ============================================================
# by Nick Matzke and Nat Hallinan (and whoever else adds to this PhyloWiki page)
# Copyright 2011-infinity
# matzkeATberkeley.edu
# February 2011
#
# Please link/cite if you use this, email me if you have 
#   thoughts/improvements/corrections.
#
# Much of the text of this lab was derived from Nat
# Hallinan's 2009 lab for IB200b on PIC in R.
#
##############################################################
#
# Free to use/redistribute under:
# Attribution-NonCommercial-ShareAlike 3.0 Unported (CC BY-NC-SA 3.0) 
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the above license, linked here:
# 
# http://creativecommons.org/licenses/by-nc-sa/3.0/
# 
# Summary:
#
# You are free:
#
#   * to Share -- to copy, distribute and transmit the work
#   * to Remix -- to adapt the work
#
# Under the following conditions:
#
#   * Attribution -- You must attribute the work in the manner 
#     specified by the author or licensor (but not in any way that 
#     suggests that they endorse you or your use of the work).
#   * Noncommercial -- You may not use this work for commercial purposes. 
#
#   * Share Alike -- If you alter, transform, or build upon this work,
#     you may distribute the resulting work only under the same or
#     similar license to this one. 
#
# http://creativecommons.org/licenses/by-nc-sa/3.0/
# 
###################################################################
#
# Assignment for Thursday, Feb. 10
#
# See the questions under "Your Assignment" near the end of the lab
#
# Purpose: Learn about the general strategies one uses to look
#   at correlation between variables; applying these to the special
#   situation of when the data is character data from species that
#   are related by a phylogeny.
# 
###################################################################

# Stuff that is handy at the beginning of a script
######################################################
# R, close all devices/windows/ quartz windows
graphics.off()

# Turn off all BS stringsAsFactors silliness
# R tends to read 
options(stringsAsFactors = FALSE)
######################################################

# Working directories:
# One of the first things you want to do, usually, is 
# decide on your working directory

# On my Mac, this is the directory
wd = "/Users/nick/Desktop/_ib200a/ib200b_sp2011/lab04"

# On a PC, you might have to specify paths like this:
#wd = "c:\\Users\\nick\\Desktop\\_ib200a\\ib200b_sp2011\\lab03"

#desc: set working directory
setwd(wd)

#desc: get working directory
getwd()

#desc: list the files in the directory
list.files()

# Sourcing scripts
# Local sourcing:
# source a script (windows machines require \\ instead of /)
# sourcedir = '/_njm/'
#source3 = '_genericR_v1.R'
#source(paste(sourcedir, source3, sep=""))

# You need the APE package to source genericR_v1.R scripts
# install.packages("ape")
library(ape)

# Remote sourcing:
source("http://ib.berkeley.edu/courses/ib200b/scripts/_genericR_v1.R")
source("http://ib.berkeley.edu/courses/ib200b/scripts/_R_tree_functions_v1.R")

# 
#########################################################
# Lab 5: Independent Contrasts
######################################################### 
#
# Today we're going to use R to derive and analyze independent contrasts for continuous
# data.  First we're going to learn about some statistical models in R by analyzing the raw
# data from our last lab.  Then we'll generate some phylogenetically independent contrasts
# and see how they affect our analyses.  Finally, in the last section you will be given a
# second data set to analyze and explain.
#
# Several other programs can do PIC, including CAIC, PDAP and the PDAP:PDTree module in
# Mesquite.  We will stick with R for now, but you may want to explore some other options
# on your own.
# 
# A word to the wise:  I will use two abbreviations in this lab, "PC" (Principle Component)
# and "PIC" (Phylogenetically independent Contrast).  These are very different things and you
# should be sure to keep them straight.
# 
########################################################## 
# Statistical Analysis in R
##########################################################
#
# First open your workspace from the last lab (or re-run that script to re-generate
# everything).  Type "ls()" to see all your objects.
# 
# Statistical Distributions
#
# R has a number of statistical distributions automatically embedded in the software.
# These can be very useful for exploring and analyzing data.  To see a list of
# distributions look at the Intro to R documentation on line.  I will use an example of the
# normal distribution, which we are all familiar with.  Later we will use the bivariate
# normal distribution and the binomial distribution.
#
# For any distribution you can add a prefix before the name of the distribution to create
# different functions:
# 
# r --- (n,parameters) will return n random draw from a distribution.
# q --- (p,parameters) will return the value for the distribution which 
#       has the given p-value.
# p --- (q,parameters) will return a p-value for a given value from the distribution.
# d --- (x,parameters) will return the probability density for a given value.
# 
# For example, type:

#desc: generate 10 random numbers, uniformly distributed between 0 and 10
(random_numbers = runif(10, min=0, max=10))

#desc: show the relative density of each of these numbers under 
#desc: two different uniform distributions
dunif(0, min=0, max=1)
dunif(0, min=0, max=10)

#desc: show the PDF (probability distribution function, i.e. probability of
#desc: getting that number or lower) for each random number, under these
#desc: uniform distributions
punif(random_numbers, min=0, max=1)
punif(random_numbers, min=0, max=10)

#desc: generate 10 numbers between 0 and 1
(probabilities = seq(0, 1, 0.1))

#desc: the values corresponding to the input 
#desc: probabilities under uniform distributions
qunif(probabilities, min=0, max=1)
qunif(probabilities, min=0, max=10)

# Let's compare our snout-vent length data to a normal distribution.  One assumption of
# linear regression is that all the parameters are normally distributed.

hist(Anole.ordered[,1])
mn = mean(Anole.ordered[,1])
stand = sd(Anole.ordered[,1])

# That looks pretty crappy.  It looks like it has tons of positive skew. Let's see what the
# p-values look like for some of our biggest numbers:

quant  =  pnorm(sort(Anole.ordered[,1]),mn,stand)
quant[28:30]

# Considering that we only have 30 values here, the fact  that 3 of them are in the top
# 98.7% percentile is more than a little suspicious.  (What percentile would you expect
# them to be in?)  The p-values of our data could be considered the expected percentage of
# data points that should be less than that value.  Furthermore the cumulative probability
# of the sampled data points represent  the percentage of data points that are less than
# that value for the actual distribution. If our data is in fact normal, then our p-values
# should match up with the cumulative probability of our data points, which are distributed
# evenly between 0 and 1.  Let's plot our p-values against some evenly distributed
# fractions and draw a straight line with slope equals 1 to see if they are the same:

plot((1:30)/31, quant)
abline(0,1)

# No, definitely not.
# 
# On the flip side let's see what values we expect at the 90th percentile:

qnorm(0.90,mn,stand)

# While we actually got:
sort(Anole.ordered[,1], decreasing=TRUE)[3]

# For the normal distribution you can actually make plots of these comparisons
# automatically using qqnorm, but you could make these plots yourself anyway for any
# distribution.
#
# Finally let's look at what this distribution would look like if it had the same
# parameters, but was normal:

hist(rnorm(30, mn, stand))
hist(rnorm(3000, mn, stand))

########################################
# Parametric Correlation of Variables
########################################
#
# Now let's plot all our data against each other to see what it looks like:

pairs(Anole.ordered)

# As you can see all these variables seam to be positively correlated.  This is not
# surprising, since they are all measures of body size.  For our initial analysis we are
# going to look at mass as a function of snout-vent length, so plot them against each
# other:

plot(Anole.ordered[,1], Anole.ordered[,2])

# To test for this correlation we will use the cor.test function based on Pearson's product
# moment correlation.

cor.test(Anole.ordered[,1],Anole.ordered[,2],alternative= "g")

# alternative="g" makes the test one tailed for a positive correlation between the
# variables.  We can make this test one tailed, because we can a priori assume a positive
# relationship between measures of size.  If you could not make this assumption, then you
# should keep the test two-tailed. 
#
# A bunch of information will appear on the screen.  At the bottom you are given the
# calculated correlation (the covariance divided by the total variance) and 95% confidence
# intervals.  Above that you see statistics testing the hypothesis that there is no
# correlation between these data.  The t-statistic is used with 28 degrees of freedom and
# the null hypothesis is completely rejected.
# 
# Assumptions of Pearson's Correlation
#
# That looks great, right? Low p-values high R2?  No, it sucks ass.  The problem is that
# the data does not match the assumptions of the model, in particular the data are not from
# a bivariate normal distribution.  That means not only is each data set normally
# distributed, but for every value of one parameter the other is normally distributed.
# First let's test the assumption that each variable is normally distributed.

qqnorm(Anole.ordered[,1])
qqline(Anole.ordered[,1])

# This is a plot of the quantiles of your data against the expected quantiles of a normal
# distribution.  qqline puts a line through the first and third quartiles.  If this
# distribution was close to normal, then all the dots would be close to the line.  As you
# can see they are not, the terminal values fall far from the line.  Repeat this for your
# mass as well.
#
# Since neither data set is normally distributed, you know that together they do not fit a
# bivariate normal distribution, but let's check anyway just for fun.  This is a little
# tricky, so I wrote a function that will plot the cumulative probability of our data
# against their p-values from the multivariate normal distribution.  The values from these
# distributions should match exactly, so we'll add a line of slope 1 through the origin for
# comparison.  Load the mnormt package...

#install.packages("mnormt")
library(mnmormt)

# ...then source the functions in "test.mvn.R" (written by previous GSI Nat Hallinan, 
# I believe)...

# Remote sourcing:
source("http://ib.berkeley.edu/courses/ib200b/scripts/test.mvn.R")

# ...and then:

qqmvn(Anole.ordered[,1:2])
abline(0,1)

# Wow, that really sucks.

################################# 
# Transforming Data
################################
#
# Luckily we do have some recourse, transforming the data.  Transforming data is a science
# and an art, we are not going to spend any time on it here.  We will just rely on the
# simple log-log transform that we already did.  Not only is this justified in that you may
# expect this type of data to be log normal, it also makes sense because these data should
# pass through the origin, so a log log transform will show a linear relationship between
# any two variables that have a relationship y=bxa, which is a more general form of the
# normal linear relationship passing through the origin. 
#
# Let's check all our assumptions for the log transformed data:

qqnorm(Anole.log[,1])
qqline(Anole.log[,1])

qqnorm(Anole.log[,2])
qqline(Anole.log[,2])

qqmvn(Anole.log[,1:2])
abline(0,1)

# Still not perfect, but over all pretty good, and a lot better than what we saw before.
# Let's test for this correlation:

cor.test(Anole.log[,1], Anole.log[,2], alternative= "g")

# That looks great, still highly correlated.
#
################################################## 
# Nonparametric Tests
#################################################
#
# So what do you do, if you can't get your data to fit the assumptions of your model?  You
# use a nonparametric test, which does not rely on any assumptions.  These tests are more
# robust than parametric tests, but less powerful.  There are a nearly infinite number of
# nonparametric tests that depend on what hypothesis you're testing.  In this case we want
# to see if one parameter is larger when the other parameter is larger.  We're going to try
# two test statistics, Kendall tau and Spearman's rho.  Both tests look to see if higher
# ranks of one variable are associated with higher or lower ranks of another.
#
# For the Kendall test:

cor.test(Anole.log[,1], Anole.log[,2], alternative= "g", method="k")

# This gives you the test statistic, z, an estimated p-value and Tau.  The p-value is
# one-tailed as before.  As you can see, the results are still highly significant, although
# the p-value is a little larger.  Check out this test for the data set that has not been
# log transformed:

cor.test(Anole.ordered[,1], Anole.ordered[,2], alternative= "g", method="k")

# The results are identical, because taking the log of a data set does not change the ranks
# of the data points.  If y > x, then ln(y) > ln(x).
#
#
# We can also do the Spearman test, which is very similar:

cor.test(Anole.log[,1], Anole.log[,2], alternative= "g", method="s")

# Obviously these two variables are highly correlated.
# 
# Testing Multiple Correlations at Once
#
# We are not limited to doing the tests for correlations between variables 2 at a time.
#
# The R package "picante" has a function cor.table, which can do every pairwise comparison
# at once.  Load picante...

install.packages("picante")
library(picante)

# ...then type:

cor.table(Anole.log)

# You will get two tables.  The first table shows the Pearson's correlation coefficient for
# each of the pairwise comparisons; the second shows the p-values for these correlations.
# Unlike cor.test, cor.table can only do two-tailed tests, since our test is one tailed, we
# can divide all these p-values by 2.  As you can see, all these correlations are highly
# significant.  However, do they match the assumptions of multivariate normality?

qqbvn(Anole.log)

# This is another little function I wrote to test these assumptions.  The little plots
# along the diagonal are comparisons of the distribution of each variable to the normal
# distribution; the other little plots to the top right are comparisons of each pair of
# variables to their expectation under a bivariate normal distribution, just like the plot
# produced by qqmvn.  The big plot on the bottom left is a comparison of all the data
# together to their expectation under a multivariate normal distribution.  The numbers
# represent p-values for the Shapiro-Wilks test of normality; significant values indicate
# rejection of the normal distribution.  Simulating data showed that both these tests may
# be a little harsh, but these results clearly indicate that very little of this data
# really fits a normal distribution.
#
# Luckily cor.table can do the same nonparametric tests as cor.test can.

cor.table(Anole.log, cor.method="k")

# and

cor.table(Anole.log, cor.method="s")

# Obviously, all these results are significant, no matter what test you use.

###############################
# Principle Components Analysis
##############################
#
# When we originally looked at the data it was clear that all the characters were highly
# correlated with size.  What if we wanted to look at changes in theses quantitative
# characters independently of that one big factor?  Well, we can look at the principle
# components of this same set of characters.  Principle components are the same data points
# but on a different set of axes.  The originally set of axes that our characters were
# measured in are rotated, so that the total variance is minimized.  This new set of axes
# separates out the parts of the individual characters that are linearly correlated with
# each other into different components.  The idea is that each of the new axes explain how
# the data varies together, rather than leaving all the data separated as they were
# measured.  
#
# Why would you want to do this?  Let's take as an example a data set with only two columns
# hind leg length and fore leg length.  One thing you could learn from a PCA is how
# correlated these characters are.  If they are totally uncorrelated (that would be a crazy
# looking animal), then each of the PCs would explain 50% of the variation.  One PC would
# be fore leg length and the other hind leg.  On the other hand, if they are highly
# correlated, then the first PC would explain most of the variation in  both fore and hind
# limbs, and would be called something like limb length.  The second PC would explain only
# a little of the variance and would correspond to something like fore leg to hind leg
# ratio.  The idea is that measurements along the PCs are somehow more informative than
# those along the original axes.  You learn more talking about limb length and fore
# leg/hind leg ratio than you do by talking about fore leg and hind leg length
# independently.  For data sets with a number of variables PCA can be quite useful for
# picking out highly informative groupings of variables.

PCs = prcomp(Anole.log)

# First let's look at some basic information about our new axes:

summary(PCs)
screeplot(PCs)

# This will show three rows describing the amount of variance described by each PC.  Look
# at the proportion of variance rows.  The plot shows the same thing.  As you can see the
# first component explains 96.7% of the variance.  That's a lot, and we can assume that
# this component is essentially size.  The standard is to say that all the components that
# cumulatively explain 95% of the variance are significant, and the rest are just a bunch
# of noise.  In this case that is just the first component.  We could just as well set it
# at 99% or whatever other arbitrarily chosen number.
#
# Now let's look at the vectors these components define:

PCs$rotation

# The first component shows positive slopes for all of our values, this is what we would
# expect for size.  The slope is particularly biased to mass, as we might expect, since
# mass should scale to the cube of length in a species of uniform shape.  On the other hand
# the lamellae are relatively unaffected.  The second component is characterized mostly by
# an increase in tail length with a slight increase in the hind legs, and interestingly a
# decrease in mass and in the fore legs.  Thus this represents some shift in overall body
# shape.  To see how the different taxa line up on these new axes type:

biplot(PCs)

# The x-axis is PC1, so taxa to your left are smaller.  The y-axis is PC2, so divergence
# along this axis shows divergence in this first shape parameter; taxa to the top have
# relatively long tails and hind legs, and relatively low mass and short fore legs for
# their size.  The vectors in the middle show how much each of our original variables
# contributes to each component. As you can see, these two components are completely
# uncorrelated, that is a consequence of PCA, and tells you nothing.
#
# To compare other PCs (like say the 2nd and 3rd) use the choices command:

biplot(PCs,choices=c(2,3))

# If for some reason you want values of each of your taxa for each of the PCs, you can find
# them in PCs$x.

################################################ 
# Phylogentically Independent Contrasts
################################################ 
# 
# Finally, we're going to get into the meat of this lab.  We've already discussed the
# theory behind PICs in lecture, so let's just get right down to it.  To calculate PICs for
# the logs of our traits on our tree:

pic(Anole.log[,1],Anole.ultra,var.contrasts=TRUE)

# This will return a matrix of 29 rows and 2 columns.  Each row refers to one internal
# node, and the name of the row gives the name of that node.  Column 1 shows the contrast
# across that node and column 2 shows it's expected variance.  Let's make an array
# containing the contrasts for all our characters and all their variances:

# create array 29x2x6
Anole.PIC = array(dim=c(29, 2, 6))

# fill it in
for (i in 1:6)
    {
    Anole.PIC[ , , i] = pic(Anole.log[,i], Anole.ultra, var.contrasts=TRUE)
    }

# "array" will create an array with the dimensions given.  for(i in 1:6) will repeat the
# functions given in the command for every element in the vector 1:6 being the value of i.
# I bet you can figure out the rest.
# 
# Confirming Assumptions of PIC
#
# So, we have our contrasts, but just like with our linear correlations, we need to make
# sure that our data fit the assumptions of the model.  Luckily David Ackerly has provided
# us with a function that makes this all very easy.  
#
# First load the picante  package.  Source the R script diagnostics_v3.R (originally 
# written by David Ackerly; modified/improved slight by Nick Matzke) and run it; if it
# doesn't work try diagnostics2.R.  

source("http://ib.berkeley.edu/courses/ib200b/scripts/diagnostics_v3.R")

# ...then type:

diagnostics(Anole.log[,1],Anole.ultra)

# Six plots will appear.  The two on the left will test whether the positivized PICs fit a
# half-normal distribution centered at 0.  The top left is just a histogram of the PICs and
# it should look like just half of a normal distribution.  The plot on the bottom right is
# a QQ-plot of the PICs against a half normal.  It should fit approximately a line.  Both
# those plots look good.
#
# You want the next three plots to show no significant relationship.  The top middle is the
# most important; it plots the standardized contrasts against their branch lengths.  If
# there is a significant negative relationship here, then that implies your branch lengths
# are bad.  Because contrasts are standardized by dividing by branch lengths, long
# arbitrary branch lengths can make the contrasts too small.  This problem can be
# alleviated by replacing all your branch lengths with their logs.  There is a negative
# slope here; it's not significant, so we wouldn't normally worry about it, but just for
# fun let's log convert the branches:

Anole.tree.log = Anole.ultra
Anole.tree.log$edge.length = log(Anole.ultra$edge.length + 1)
diagnostics(Anole.log[,1], Anole.tree.log)

# The p-value is now even higher for that relationship.  The plot on the top right shows
# the relationship between the contrasts and the age of the node; 0 is the time of the root
# and the present can be found to the right.  A positive relationship indicates extremely
# low signal and a negative extremely high.  The plot in the bottom middle shows the
# relationship between the reconstructed value at a node and the contrasts; a positive
# relationship implies that the data should be log transformed, before the contrast is
# taken.  To see an example of this look at the untransformed mass contrasts to the
# transformed ones:

diagnostics(Anole.ordered[,2],Anole.ultra)
diagnostics(Anole.log[,2],Anole.ultra) 

# See, more justification for a log transformation.  The plot on the bottom right is of
# course the plots of k that we learned about during lecture on Thursday.  Values of k
# close to 1 show the expected amount of signal under Brownian Motion.  If you have a lot
# of taxa, you may want to suppress the calculation of k with "calcK=FALSE", as it can be
# cumbersome.
#
# Investigate the contrasts for our other traits.  Are you satisfied that they fit the
# assumptions?  Great then let's move on.
# 
# 
# Parametric Significance of PICs
#
# The first thing that we want to do is test for a significant relationship between the
# PICs.  It is important to recognize the difference between the relationship between the
# raw data and the relationship between the PICs.  In the first case you are asking whether
# taxa with greater length have larger mass; in the second your question is do evolutionary
# increases in length correspond with increases in mass.  
# As always let's start by plotting all our factors against each other:

pairs(Anole.PIC[,1,])

# They all look very correlated.  First we'll make sure that our contrasts match the
# assumptions of Pearson's test.  We've already concluded that our data is normal from the
# half-normal plots in the tests of PIC assumptions.  So all we need to test is that our
# data fits a  bivariate normal distribution.  Remember that PICs must have a mean of 0, so
# we're going to have to have to force the fit through the origins by setting the means
# ourselves:

qqbvn(Anole.PIC[,1,],means=rep(0,6))

# That looks pretty good.  Why would you expect the PICs to fit a multivariate normal
# distribution better than the raw data?  The p-values for the Shapiro-Wilks test are not
# displayed, because its assumptions are not met by PICs, so we're going to have to go off
# of fit alone.  I might be concerned about the fits of all the comparisons with the
# lamellae or the tail lengths, but the rest look great.
#
# Now we can look for a linear correlation.  We also have to force our fit through the
# origin here.  cor.test can't do that, but cor.table can, if we use the cor.type="c" (for
# contrast) command:

cor.table(Anole.PIC[,1,],cor.type="c")

# As you can see all these results are very significant. Keep in mind that we made 15
# separate comparisons; you would expect to find one comparison with a p-value less than
# 0.05 more than 1 in 20 times, when you look at 15 comparisons.  As a general rule you
# should multiply your p-values by the number of comparisons you make.  The number of
# comparisons can easily be calculated as n(n-1)/2, where n is the number of variables that
# we have.  Let's do that to see the effects on our results.

results.cor = cor.table(Anole.PIC[,1,], cor.type="c")
results.cor$P * 15/2

# We divided by two, because these p-values are for a two tailed-test and our test is one
# tailed.  Of course our results are still all significant.  They were so significant
# before, that we wouldn't expect multiplying by 7.5 to have much of an effect.  For a
# quick view of which comparisons are significant try:

(results.cor$P * 7.5) < 0.05

# Nonparametric Significance of PICs
#
# As with raw data, PICs may not always fit the assumptions of linear regression.  In that
# case you need to do a nonparametric test.  The first and most basic test is to ask
# whether one character increases, when the other increases.  This is a very straight
# forward test.  The first thing to do is create a logical vector that is true when both
# PICs are positive or negative, false when one is negative and the other positive, and
# excludes contrasts for which either character shows no change.

ref = (Anole.PIC[,1,1] != 0) && (Anole.PIC[,1,2] != 0)
both.pos = (Anole.PIC[ref,1,1] * Anole.PIC[ref,1,2]) > 0
length(both.pos)
sum(both.pos)

# For 26 out of 29 comparisons between contrasts either both show an increase or a
# decrease.  We can test the significance of this result by comparing it to the binomial
# distribution.  The binomial distribution is a distribution that takes one of two values.
# The probability of getting one result or the other is one of the parameters in the model.
# We will assume that there is a 50% chance of getting each result.  In that way we will
# try to reject the null hypothesis that it is equally likely for the contrasts to change
# in the same direction as it is for them to change in the opposite.

pbinom(3,29,0.5) 

# We used 29-26=3, instead of 26, because this formula evaluates the probability of getting
# that many results or fewer, so for a one tailed test of positive correlation, we want to
# know the probability of getting as many negative results or fewer.  As you can see this
# relationship is highly significant.
#
# cor.table also has been programmed for you to do non-parametric rank correlation tests
# for independent contrasts.  As with parametric tests, you must use the cor.type="c"
# command:

cor.table(Anole.PIC[,1,], cor.method="k", cor.type="c")

# and

cor.table(Anole.PIC[,1,], cor.method="s", cor.type="c")

# Are all those results significant?  Don't forget to multiply by 7.5.  (Would you like to
# make a chart of all the possible comparisons between PICs using the binomial
# distribution? If so see the appendix.)
# 
# PCA on PICs
#
# You can also generate PCs from PICs.  Why would you want to do that?  Why is that any
# better than doing it on the raw data?  The idea is that the PCs of the PICs represent the
# principle axes of evolutionary change.  In other words, the ways that things actually
# change when they evolve.  This can obviously be a very interesting bit of information.
# We have to do one little trick to do PCA on PICs to force them through the origin.  We
# are going to create a mirror image of the PICs and stick that together with the other
# PICs.  This will double the amount of data, but since PCA does not do any statistical
# tests, degrees of freedom don't matter, so we'll be fine.  

PIC.mirror = rbind(Anole.PIC[,1,], (Anole.PIC[,1,] * (-1)) )
PCs.PIC = prcomp(PIC.mirror)

# Now let's check them out:

summary(PCs.PIC)

# As you can see, the first PC still describes the vast majority of the variance, and it is
# probably change in size.  Let's look at the vectors:

PCs.PIC$rotation

# These look very much like the original PCs that we generated. For some data sets there
# can be large differences.  When would you expect to see such a difference?
# Finally, let's plot them:

biplot(PCs.PIC)

# or

biplot(PCs.PIC,choices=c(2,3))

# You will see your contrasts spread out over the PCs.  This tells you how much change of
# each type happened along each contrast.  Remember you mirrored your contrasts, so each
# one will appear twice on opposite sides of the plot.
# 
# Your Assignment
# 
# Here is your assignment for this week.  I'm really going to grade this one, not just mark
# off that you did it, so do a good job.  Turn it in to me on paper in class next Tuesday,
# the 24th of February.  Please try to be concise.
# 
# You should already have downloaded the nexus file PIC.assignment from the web site.  This
# file has a tree with 49 taxa and a data matrix with 3 continuous characters.  This is a
# modified version of a file from the Mesquite examples.
# 
# 1.  Are these characters correlated without considering the phylogeny?  What are the
# p-value of the correlations?  Is the correlation positive or negative?  What statistical
# test did you use?  Why did you choose that test?  Was the test two-tailed or one-tailed?
# Why?
# 
# 2.  Are these characters correlated, when you do consider the phylogeny?  Does the data
# fit the assumptions of independent contrasts?  What are the p-value of the correlations?
# Is the correlation positive or negative?  What statistical test did you use?  Why did you
# choose that test?  Was the test two-tailed or one-tailed?
# 
# 3  Why is there a difference in result between these two tests for significance? Be
# specific; you should both look over the contrast plots and plot reconstructions of all
# these characters on a tree (in R or Mesquite).  Does this say anything interesting about
# the evolution of these taxa?
# 
# You should always ask yourself these questions for any analysis involving PICs.

#######################################################################
# Appendix: A Brief Digression on Writing Functions
#######################################################################
#
# Up until now I've avoided the subject of writing functions.  Let's write a function
# that will produce a table of results from a two-tailed sign test for every possible
# comparison between traits.  I just need to explain a few things first.
# In R functions are treated as objects, just like variables.  When you write a function,
# you write "function" followed by parentheses, which have the possible input parameters.
# You can define defaults for those parameters at this time.  You also assign that function
# to a named object like any other object in R.  The function will then appear in your list
# of objects.
#
# "for (num in vector)" will run the commands that follow it for every element in the
# vector.  Replacing num with that element.
#
# "if" is followed by some conditional statement in parentheses.  The subsequent commands
# are only carried out if the conditional statement is true.  If it is not true the
# commands following "else" are executed.
#
# All of these execute a series of commands that appear after them; by after I mean either
# on the same line, or surrounded by "{}".

# Create a function object named sign.test 
#   with one input paramater, x
sign.test = function(x)
    {
    # How many columns in x
    n = dim(x)[2]                     

    # Create two n by n matrices
    counts = ps = array(0,dim=c(n,n))

    # logical vector indicating which changes are 0
    not0 = x!=0

    # start a loop that will run variable i from 1 to n-1
    for (i in 1:(n-1))
        {
        for (j in (i+1):n)
            {
            # Neither should be 0
            ref = not0[,i] && not0[,j]
            both.pos = (x[ref,i] * x[ref,j]) > 0

            # The number of positive results
            pos = sum(both.pos)

            # Add those to above the diagonal
            counts[i,j] = pos

            # The number of negative results
            neg = length(both.pos)-pos

            # Add those to below the diagonal
            counts[j,i] = neg

            # For 2-tailed want pos or neg correlation
            times = min(pos,neg)

            # we do the test
            # times 2 cause two-tailed
            ps[i,j] = 2*pbinom(times,pos+neg,0.5)

            # Below diahonal matches the top
            ps[j,i] = ps[i,j]
        }
    }
    # return your results as a list
    list_to_return = list(counts=counts,ps=ps)
    return(list_to_return)
}

# Now we can run the test on all our data:

sign.test(Anole.PIC[,1,])

# There, don't you feel accomplished.  The first table shows the number of positive results
# above the diagonal and the number of negative results below.  The second shows the
# p-values.  As you can see many of our correlations are significant, but not all.  Don't
# forget to multiply by 7.5.  To see exactly which are significant.

(sign.test(Anole.PIC[,1,])$ps * 7.5) < 0.05

# How does that compare to the other non-parametric tests?