wrightscape: wrightscape

# correct_labrid_data.R
require(wrightscape)
require(phytools)
require(geiger)

# This data has not been released
labrid_tree <- read.nexus("../data/labrid_tree.nex")
diet_data <- read.csv("../data/labriddata_parrotfish.csv")

corrected_data <- diet_data
# Use the simple names from Price et al 2010
traitnames <- c("Species", "group", "gape", "prot", "bodymass", "AM", "SH", "LP", "close", "open", "kt")
names(corrected_data) <- traitnames

# Lengths are log transformed 
corrected_data[["gape"]] <- log(corrected_data[["gape"]])
corrected_data[["prot"]] <- log(corrected_data[["prot"]])
#masses are log(cube-root) transformed
corrected_data[["bodymass"]] <- log(corrected_data[["bodymass"]])/3
corrected_data[["AM"]] <- log(corrected_data[["AM"]])/3
corrected_data[["SH"]] <- log(corrected_data[["SH"]])/3
corrected_data[["LP"]] <- log(corrected_data[["LP"]])/3
# ratios are fine as they are

# We'll get just the parrotfish data by taking only the traits for parrotfish.  
# treedata() will automatically drop the tips with wrasses from the phylogeny.  

parrotfish <- corrected_data[corrected_data[,2]=="parrotfish",]
ape <- treedata(labrid_tree, parrotfish[,3:11], parrotfish[,1])

ape$data["bodysize"]
out <- phyl.resid(ape$phy, ape$data[,"bodymass"], ape$data[,c("gape", "prot","AM", "SH", "LP")] )
## phyl.resid changes order of species listing. Merge for a set of uncorrected and corrected traits.  
traits <- merge(ape$data, out$resid, by="row.names")
# columns that are transformed now have gape.x for untransformed, gape.y for transformed.  

# format data gets regimes from column specified in "regimes" (e.g. this is a column id, not a # of regimes)
# This also converts the tree and data into ouch format
## Could just hand it all traits, but these are just the tranformed and size-corrected ones
labrid <- format_data(labrid_tree, traits[,2:length(traits)], species_names=traits[,1])  
#names(labrid$data)

## PAINTING REGIMES: Having taken care of traits, we paint on the 3 regime models.  
# Select common ancestor of a Chlorurus and a Hipposcarus as the changepoint
intra_ancestor <- mrcaOUCH(c("Chlorurus_sordidus", "Hipposcarus_longiceps"), labrid$tree)
intramandibular <- paintBranches(intra_ancestor, labrid$tree, c("other","intramandibular"))


# rename the ouch-formated data 
dat <- labrid$data
tree <- labrid$tree

# rename the ape-formatted data
ape.phy <- ape$phy
ape.dat <- traits


save(list=c("intramandibular", "tree", "dat", "ape.phy", "ape.dat"), file="parrotfish.rda")