#-----------------------------------------# # Script to calculate PIC for in class # example and conduct a regression analysis # Tree and data = in Ginger's lecture notes # # Ginger Jui # Feb. 10, 2011 #-----------------------------------------# rm(list=ls()) # clear the workspace in R library(ape) # load the ape package #-----------------------------------------# # First, input the data into R #-----------------------------------------# # Tree with all branchlengths = 1 newick = "(((t1:1,t2:1):1,(t3:1,t4:1):1):1,((t5:1,t6:1):1,(t7:1,t8:1):1):1):1;" # Read in tree as phylo object phy <- read.tree(text=newick) # Plot the tree and label the internal nodes plot(phy) nodelabels() # Create vectors with the trait values X <- c(2,6,1,11,8.5,11.5,19,13) Y <- c(8,16,6,14,9,11,18,12) # Names of vectors should be same as tips names(X) <- names(Y) <- c("t1", "t2", "t3", "t4", "t5", "t6", "t7","t8") #-----------------------------------------# # Calculate phylogenetic independent contrasts #-----------------------------------------# pic.X <- pic(X, phy, var.contrasts=T) pic.Y <- pic(Y,phy, var.contrasts=T) var.con.X <- pic.X[,2] var.con.Y <- pic.X[,2] std.con.X <- pic.X[,1] std.con.Y <- pic.Y[,1] #-----------------------------------------# # Plot the contrasts, and a regression # line through the origin #-----------------------------------------# plot(std.con.Y~std.con.X) abline(lm(std.con.Y~std.con.X-1)) #-----------------------------------------# # Statistics #-----------------------------------------# # Pearson's Correlation Coefficient = .84; # not calculated through origin! cor.test(std.con.X, std.con.Y) # Correlation Coefficient through the origin = .87 sum(std.con.X*std.con.Y)/sqrt((sum(std.con.X2)*sum(std.con.Y2)))