Challenge problem 4: Discrete character evolution solution

Here is a solution to the challenge problem. First we need to get the data prepped - you probably already have this from the exercise!

Files here:

1. brandley_table.csv
2. squamate.tre

library(ape)
library(phytools)

## Loading required package: maps

## 
##  # maps v3.1: updated 'world': all lakes moved to separate new #
##  # 'lakes' database. Type '?world' or 'news(package="maps")'.  #

library(geiger)

sqData<-read.csv("brandley_table.csv")
sqTree<-read.nexus("squamate.tre")
hindLimbs<-sqData[,"HLL"]!=0
foreLimbs<-sqData[,"FLL"]!=0

speciesNames<-sqData[,1]
species<-sub(" ", "_", speciesNames)

limbState<-paste(foreLimbs, hindLimbs)
names(limbState)<-species

Run the previous models:

limbState<-paste(foreLimbs, hindLimbs)
names(limbState)<-species

erModel<-fitDiscrete(sqTree, limbState, model="ER")

## Warning in treedata(phy, dat): The following tips were not found in 'phy' and were dropped from 'data':
##  Gonatodes_albogularis
##  Lepidophyma_flavimaculatum
##  Trachyboa_boulengeri

symModel<-fitDiscrete(sqTree, limbState, model="SYM")

## Warning in treedata(phy, dat): The following tips were not found in 'phy' and were dropped from 'data':
##  Gonatodes_albogularis
##  Lepidophyma_flavimaculatum
##  Trachyboa_boulengeri

ardModel<-fitDiscrete(sqTree, limbState, model="ARD")

## Warning in treedata(phy, dat): The following tips were not found in 'phy' and were dropped from 'data':
##  Gonatodes_albogularis
##  Lepidophyma_flavimaculatum
##  Trachyboa_boulengeri

Now, let’s create a design matrix for the character evolution. This is reasonably straightforward except we need to watch out for one transition that is from TRUE FALSE to FALSE TRUE - a gain accompanied by a loss - which we will not allow.

r1Matrix<-cbind(c(NA, 1, 2, 3), c(0, NA, 0, 4), c(0, 0, NA, 5), c(0, 0, 0, NA))
r1Model<-fitDiscrete(sqTree, limbState, model=r1Matrix)

## Warning in treedata(phy, dat): The following tips were not found in 'phy' and were dropped from 'data':
##  Gonatodes_albogularis
##  Lepidophyma_flavimaculatum
##  Trachyboa_boulengeri

## Warning in fitDiscrete(sqTree, limbState, model = r1Matrix): Parameter estimates appear at bounds:
##  q12
##  q13
##  q14
##  q23
##  q24
##  q32
##  q34

r2Matrix<-r1Matrix
r2Matrix[4,1]<-0
r2Model<-fitDiscrete(sqTree, limbState, model=r2Matrix)

## Warning in treedata(phy, dat): The following tips were not found in 'phy' and were dropped from 'data':
##  Gonatodes_albogularis
##  Lepidophyma_flavimaculatum
##  Trachyboa_boulengeri

## Warning in fitDiscrete(sqTree, limbState, model = r2Matrix): Parameter estimates appear at bounds:
##  q12
##  q13
##  q14
##  q23
##  q24
##  q32
##  q34
##  q41

# compare all


aicc <- c(erModel$opt$aicc, symModel$opt$aicc, ardModel$opt$aicc, r1Model$opt$aicc, r2Model$opt$aicc)
aiccD <- aicc - min(aicc)
aw <- exp(-0.5 * aiccD)
aiccW <- aw/sum(aw)
aiccW

## [1] 3.258643e-05 5.016437e-02 3.062099e-03 9.080245e-01 3.871642e-02

90% AIC weight for loss only model, and about 94% for the two new models!