## load packages, trees, and data. For this package only I parsed down the data file to only two columns (species and activity state) and created separate files for diurnal and nocturnal. This made referencing the data simpler in the corHMM function. library(ape) library(corHMM) ericson <- read.nexus("Ericson1914.tree") hackett <- read.nexus("Hackett1914.tree") dataD <- read.csv("appendixD.csv") dataN <- read.csv("appendixN.csv") ## run corHMM hidden rates model to determine fit, using default root state. Compare different numbers of hidden rates by running the model with all possible rate categories, from 1 to 5. The example listed below is for a 2 rate model. To run all additional corHMM models, simply change the rate.cat function. More rate catergories take more time to run. ## ericson nocturnal, 2 rates EN2.hmm<-corHMM(ericson, dataN[,c(1,2)], rate.cat = 2, node.states = c("marginal")) ## ericson diurnal, 2 rates ED2.hmm<-corHMM(ericson, dataD[,c(1,2)], rate.cat = 2, node.states = c("marginal")) ## hackett nocturnal, 2 rates HN2.hmm<-corHMM(hackett, dataN[,c(1,2)], rate.cat = 2, node.states = c("marginal")) ## hackett diurnal, 2 rates HD2.hmm<-corHMM(hackett, dataD[,c(1,2)], rate.cat = 2, node.states = c("marginal")) ## compare corHMM models with different rates using the AICc and loglik values within the same tree and coding scheme only EN1.hmm$loglik EN1.hmm$AICc EN2.hmm$loglik EN2.hmm$AICc ## test the affect of different root treatments on the fit of the model. The default in corHMM is root.p=NULL, assuming an equal weighting of all states at the root. The methods "yang" and "madfitz" can also be applied to weight the likelihoods of each state according to transition rates on the tree. This was only done on models using the best fitting number of hidden rate classes. ## testing "yang" root treatment ## ericson nocturnal EN2.yang <-corHMM(ericson, dataN[,c(1,2)], rate.cat = 2, node.states = c("marginal"), root.p="yang") ## ericson diurnal ED2.yang <-corHMM(ericson, dataD[,c(1,2)], rate.cat = 2, node.states = c("marginal"), root.p="yang") ## hackett nocturnal HN2.yang <-corHMM(hackett, dataN[,c(1,2)], rate.cat = 2, node.states = c("marginal"), root.p="yang") ## hackett diurnal HD3.yang <-corHMM(hackett, dataD[,c(1,2)], rate.cat = 3, node.states = c("marginal"), root.p="yang") ## testing "madfitz" root treatment ## ericson nocturnal EN2.madfitz <-corHMM(ericson, dataN[,c(1,2)], rate.cat = 2, node.states = c("marginal"), root.p="madfitz") ## ericson diurnal ED2.madfitz <-corHMM(ericson, dataD[,c(1,2)], rate.cat = 2, node.states = c("marginal"), root.p="madfitz") ## hackett nocturnal HN2.madfitz <-corHMM(hackett, dataN[,c(1,2)], rate.cat = 2, node.states = c("marginal"), root.p="madfitz") ## hackett diurnal HD3.madfitz <-corHMM(hackett, dataD[,c(1,2)], rate.cat = 3, node.states = c("marginal"), root.p="madfitz") ## for the best fitting number of rate classes and the best fitting root treatment, plot ancestral state reconstructions as pie charts at the nodes of the phylogeny use plotRECON ## colors are ordered (0_R1, 1_R1, 0_R2, 1_R2) with R1 (rate 1) being slower than R2 ## ericson nocturnal plotRECON(ericson, EN2.madfitz$states, piecolors = c("yellow", "cornflowerblue", "orange", "dark blue"), cex = 0.1, pie.cex = 0.15, file = "EricsonNoc2.corHMM.pdf", height = 50, width = 12, show.tip.label = TRUE, label.offset = 0.4, title = "Ericson Nocturnal Two Rates") ## ericson diurnal plotRECON(ericson, ED2.madfitz$states, piecolors = c("yellow", "cornflowerblue", "orange", "dark blue"), cex = 0.1, pie.cex = 0.15, file = "EricsonDiu2.corHMM.pdf", height = 50, width = 12, show.tip.label = TRUE, label.offset = 0.4, title = "Ericson Diurnal Two Rates") ## hackett nocturnal plotRECON(hackett, HN2.madfitz$states, piecolors = c("yellow", "cornflowerblue", "orange", "dark blue"), cex = 0.1, pie.cex = 0.15, file = "HackettNoc2.corHMM.pdf", height = 50, width = 12, show.tip.label = TRUE, label.offset = 0.4, title = "Hackett Nocturnal Two Rates") ## hackett diurnal plotRECON(hackett, HD3.madfitz$states, piecolors = c("yellow", "cornflowerblue", "goldenrod1", "blue2", "orange", "dark blue"), cex = 0.1, pie.cex = 0.15, file = "HackettDiu3.corHMM.pdf", height = 50, width = 12, show.tip.label = TRUE, label.offset = 0.4, title = "Hackett Diurnal Three Rates") ## end