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:
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!