# Script to run PGLS analysis for fractal data in equids. Code was written
#by Nicholas A. Famoso on October 11th, 2015. This code requires .csv files
#for the raw data, a .nex file for the phylogeny.

#This code also has code for the t-test between the slopes of the PGLS for
#each tribe. Phylogenetic signal is also calculated. The code for the 
#phylogeny figure is also present.


#Clear the memory
rm(list=ls())

#Bring in data
Data<- read.csv("Fracal Data 2015.csv")
D<-Data$D
Tribe<-Data$Tribe
Area<-Data$Area
Genus<-Data$Genus
Species<-Data$Species
Age<-Data$Geologic_Age


#test for normality in D
shapiro.test(D)
#it's not significant so that means D is normally distributed 
#Test for homogeneity in D
bartlett.test(D~Tribe) 

#lets do the PGLS
horse.trees <- read.nex("horse.trees.nex", sep=";")

#load the libraries
library(caper)
library(paleotree)
#read tree from nexus file
tree<- read.nexus("horse.nex")
#read dates from csv.
PANCOVADATA<-read.csv("Dates_and_D.csv")
#make the names of the taxa into the rownames of the dataframe
rownames(PANCOVADATA)<-PANCOVADATA[,1]
#deleate the old column of taxon names
PANCOVADATA[,1]<-NULL
#make the data frame into a matrix
PANCOVADATA<-as.matrix(PANCOVADATA)
#make the paleo tree, zlba helps to correct for 0 length branchs and provides
#a 1 million year cushion
paleoTree<-timePaleoPhy(tree, PANCOVADATA, type = "zlba", vartime = 1, add.term = T)

#bring in data
stuff<-read.csv("D_Area.csv")

#creates dataframe
DF.Horse<-data.frame(stuff)
#makes compairative matrix
Horse <- comparative.data(phy = paleoTree, data = DF.Horse,
                          names.col = Taxon , vcv = TRUE,
                          na.omit = FALSE, warn.dropped = TRUE)
#Do the Maximum Likelyhood PGLS
MLmodel.pgls <- pgls(D~Area, data = Horse, lambda = "ML")
#print summary of PGLS
summary(MLmodel.pgls)
#plots PGLS result
plot(D~Area, data = Horse$data, pch=16)
abline(MLmodel.pgls)
title("Equidae")

#Now I need to do all of this for each tribe
#First Hipparionini
#load the libraries
library(caper)
library(paleotree)
#read tree from nexus file
hiptree<- read.nexus("Hipparion.nex")
#read dates from csv.
HippData<-read.csv("Hip_Dates.csv")
#make the names of the taxa into the rownames of the dataframe
rownames(HippData)<-HippData[,1]
#deleate the old column of taxon names
HippData[,1]<-NULL
#make the data frame into a matrix
HippData<-as.matrix(HippData)
#make the paleo tree, zlba helps to correct for 0 length branchs and provides
#a 1 million year cushion
HippaleoTree<-timePaleoPhy(hiptree, HippData, type = "zlba", vartime = 1, add.term = T)

#bring in data
hipstuff<-read.csv("Hipp_D_Area.csv")

#creates dataframe
DF.HipHorse<-data.frame(hipstuff)
#makes compairative matrix
HipHorse <- comparative.data(phy = HippaleoTree, data = DF.HipHorse,
                             names.col = Taxon , vcv = TRUE,
                             na.omit = FALSE, warn.dropped = TRUE)
#Do the Maximum Likelyhood PGLS
HipMLmodel.pgls <- pgls(D~Area, data = HipHorse, lambda = "ML")
#print summary of PGLS
summary(HipMLmodel.pgls)
#plots PGLS result
plot(D~Area, data = HipHorse$data, pch=16)
abline(HipMLmodel.pgls)
title("Hipparionini")

#Now the Equini
#load the libraries
library(caper)
library(paleotree)
#read tree from nexus file
Eqtree<- read.nexus("Equini.nex")
#read dates from csv.
EqData<-read.csv("Eq_Dates.csv")
#make the names of the taxa into the rownames of the dataframe
rownames(EqData)<-EqData[,1]
#deleate the old column of taxon names
EqData[,1]<-NULL
#make the data frame into a matrix
EqData<-as.matrix(EqData)
#make the paleo tree, zlba helps to correct for 0 length branchs and provides
#a 1 million year cushion
EqpaleoTree<-timePaleoPhy(Eqtree, EqData, type = "zlba", vartime = 1, add.term = T)

#bring in data
Eqstuff<-read.csv("Eq_D_Area.csv")

#creates dataframe
DF.EqHorse<-data.frame(Eqstuff)
#makes compairative matrix
EqHorse <- comparative.data(phy = EqpaleoTree, data = DF.EqHorse,
                            names.col = Taxon , vcv = TRUE,
                            na.omit = FALSE, warn.dropped = TRUE)
#Do the Maximum Likelyhood PGLS
EqMLmodel.pgls <- pgls(D~Area, data = EqHorse, lambda = "ML")
#print summary of PGLS
summary(EqMLmodel.pgls)
#plots PGLS result
plot(D~Area, data = EqHorse$data, pch=16)
abline(EqMLmodel.pgls)
title("Equini")

#t test of slopes
b1=0.013374
b2=-0.0158586
Sb1=0.018716
Sb2=0.0071118
n1=12
n2=22

sb1b2=sqrt((Sb1^2)+(Sb2^2))

t=(b1-b2)/(sb1b2)
df=((n1+n2)-4)
#critical t values two tailed at alpha=0.05 is 2.042

#is there phylo signal in the characters independent of one another
#just equini
Eqest.lambdaD <- pgls(D~1, data = EqHorse, lambda = "ML")
summary(Eqest.lambdaD)

Eqest.lambdaA <- pgls(Area~1, data = EqHorse, lambda = "ML")
summary(Eqest.lambdaA)

#just hipparions
Hipest.lambdaD <- pgls(D~1, data = HipHorse, lambda = "ML")
summary(Hipest.lambdaD)

Hipest.lambdaA <- pgls(Area~1, data = HipHorse, lambda = "ML")
summary(Hipest.lambdaA)

#entire system
est.lambdaD <- pgls(D~1, data = Horse, lambda = "ML")
summary(est.lambdaD)

est.lambdaA <- pgls(Area~1, data = Horse, lambda = "ML")
summary(est.lambdaA)

#map characters of tree
library(phytools)
#obj<-contMap(EqpaleoTree,EqD,lwd=6,res=200,legend=1.5)
D_Area<-read.csv("C:/Users/NickF/Google Drive/Horse Fractals/D_Area_1.csv")

d<-setNames(D_Area$D,D_Area$Taxon)
area<-setNames(D_Area$Area,D_Area$Taxon)

obj_D<-contMap(paleoTree,d,lwd=6,res=200,legend=20)
title("D")
obj_Area<-contMap(paleoTree,area,lwd=6,res=200,legend=20)
title("Area")