## load packages
library(ape) 
library(geiger)
library(caper) 
## load tree and data 
ericson <- read.nexus("Ericson1914.tree") 
hackett <- read.nexus("Hackett1914.tree") 
labels <- read.csv("appendix1914.csv")
Nstates <- labels$ActNoc
Dstates <- labels$ActDiu
names(Nstates) <- labels$Species
names(Dstates) <- labels$Species 

## Run fitDiscrete function using estimated lambda transformation 
## ericson nocturnal 
EN.lambda<- fitDiscrete(ericson, Nstates, model = c("ARD"), transform = c("lambda")) 
## ericson diurnal 
ED.lambda<- fitDiscrete(ericson, Dstates, model = c("ARD"), transform = c("lambda"))
## hackett nocturnal 
HN.lambda<- fitDiscrete(hackett, Nstates, model = c("ARD"), transform = c("lambda"))
## hackett diurnal 
HD.lambda<- fitDiscrete(hackett, Dstates, model = c("ARD"), transform = c("lambda"))

## Run fitDiscrete function using the white-noise transformation 
EN.wn<- fitDiscrete(ericson, Nstates, model = c("ARD"), transform = c("white")) 
## ericson diurnal 
ED.wn<- fitDiscrete(ericson, Dstates, model = c("ARD"), transform = c("white"))
## hackett nocturnal 
HN.wn<- fitDiscrete(hackett, Nstates, model = c("ARD"), transform = c("white"))
## hackett diurnal 
HD.wn<- fitDiscrete(hackett, Dstates, model = c("ARD"), transform = c("white"))

## Run phylo.d to calculate D-statistic with caper
## create data frame
dfN <- data.frame(species, Nstates)
dfD <- data.frame(species, Dstates)
## ericson nocturnal D-statistic
EN.Dstat <- phylo.d(Nstates, ericson, names.col = species, binvar = Nstates) 
## ericson diurnal 
ED.Dstat <- phylo.d(Dstates, ericson, names.col = species, binvar = Dstates)
## hackett nocturnal 
HN.Dstat <- phylo.d(Nstates, hackett, names.col = species, binvar = Nstates)
## hackett diurnal 
HD.Dstat <- phylo.d(Dstates, hackett, names.col = species, binvar = Dstates)

## end