

require(ape)
require(phytools)
require(geiger)
require (mvMORPH)
require(devtools)
require(nlme)
require(qpcR)
require(reshape)
require(ouch)
require(calibrate)
require(surface)
require(PhylogeneticEM)

tree<- read.nexus("Chapter 1_trees.nexus")
dat<- read.csv("Chapter 1_data_all PCs.csv", row.names=1)
dd<-read.csv("Chapter 1_PhylogeneticEM_data.csv", check.names=FALSE)

TreeOnly <- setdiff(tree$tip.label,rownames(dat)) #checks to make sure the tips of the phylogeny matches the morpho data

DataOnly <- setdiff(rownames(dat), tree$tip.label) #checks to make sure the tips of the phylogeny matches the morpho data

#########################################
#Prune Tree
########################################

phyloTime <- drop.tip(tree,TreeOnly)

phy<-ladderize(phyloTime)


#########SURFACE

res<- PhyloEM(phylo = phy,
				Y_data= dd,
				process= "scOU",
				#independent=TRUE,
				#random.root=TRUE,
				#method.selection = c("LINselect","DDSE", "Djump"),
				K_Max=17)

#method.init.alpha = "estimation",
				#method.init.alpha.estimation = c("regression", "regression.MM", "median"),

plot(res)