The problem was basically required that we use the same data set for hindfoot
digit evolution from Exercise 18 tofit a directional ordered model for
toe loss (i.e., 5→4→3→2→1→0) and compare this model
to the other three models we fit: "ER", "ARD", and
a model for bi-directional, ordered, toe loss & gain.
The data were as follows:
First, let's load our packages & data:
library(phytools)
## Loading required package: ape
## Loading required package: maps
library(geiger)
alldata<-read.csv("brandley_table.csv")
squamates<-read.nexus("squamate.tre")
plotTree(squamates,type="fan",lwd=1,fsize=0.4,ftype="i")
Next, let's pull out the variable hf.digits, just as we did
in the previous exercise:
hf.digits<-setNames(as.factor(round(alldata$Toes)),
sub(" ","_",alldata$Species))
hf.digits
## Agamodon_anguliceps Amphisbaena_alba
## 0 0
## Bipes_biporus Bipes_caniculatus
## 0 0
## Bipes_tridactylus Blanus_cinereus
## 0 0
## Chirindia_swimmertoni Diplometopon_zarudnyi
## 0 0
## Geocalamus_acutus Monopeltis_capensis
## 0 0
## Rhineura_floridana Trogonophis_wiegmanni
## 0 0
## Abronia_graminea Anguis_fragilis
## 5 0
## Anniella_geronimensis Anniella_pulchra
## 0 0
## Barisia_imbricata Celestus_enneagrammus
## 5 5
## Diploglossus_bilobatus Diploglossus_pleei
## 5 5
## Elgaria_coerulea Elgaria_kingii
## 5 5
## Elgaria_multicarinata Elgaria_panamintina
## 5 5
## Elgaria_paucicarinata Gerrhonotus_liocephalus
## 5 5
## Mesaspis_moreleti Ophiodes_striatus
## 5 1
## Ophisaurus_apodus Ophisaurus_attenuatus
## 1 0
## Ophisaurus_harti Ophisaurus_koellikeri
## 0 1
## Ophisaurus_ventralis Sauresia_agasepsoides
## 0 4
## Wetmorena_haetiana Chamaesaura_anguina
## 4 1
## Cordylus_cataphractus Cordylus_cordylus
## 5 5
## Cordylus_jordani Cordylus_warreni
## 5 5
## Platysaurus_rhodiseinsis Pseudocordylus_microlepidotus
## 5 5
## Dibamus_novaeguineae Coleonyx_elegans
## 0 5
## Diplodactylus_damaeus Gekko_gecko
## 5 5
## Gonatodes_albogularis Oedura_coggeri
## 5 5
## Teratoscincus_scincus Angolosaurus_skoogi
## 5 5
## Cordylosaurus_trivittata Gerrhosaurus_major
## 5 5
## Gerrhosaurus_nigrolineatus Tetradactylus_africanus
## 5 1
## Tetradactylus_seps Tetradactylus_tetradactylus
## 5 4
## Tracheloptychus_madagascariensis Zonosaurus_ornatus
## 5 5
## Alopoglossus_atriventris Alopoglossus_carinicaudatus
## 5 5
## Alopoglossus_copii Arthrosaura_kockii
## 5 5
## Arthrosaura_reticulata Bachia_bresslaui
## 5 1
## Bachia_dorbignyi Bachia_flavescens
## 1 1
## Calyptommatus_leiolepis Calyptommatus_nicterus
## 1 1
## Calyptommatus_sinebrachiatus Cercosaura_argulus
## 1 5
## Cercosaura_eigenmanni Cercosaura_ocellata
## 5 5
## Cercosaura_quadrilineata Cercosaura_schreibersii
## 5 5
## Colobodactylus_dalcyanus Colobodactylus_taunayi
## 5 5
## Colobosaura_modesta Ecpleopus_affinis
## 5 5
## Gymnophthalmus_leucomystax Heterodactylus_imbricatus
## 5 5
## Iphisa_elegans Leposoma_percarinatum
## 5 5
## Micrablepharus_maximiliani Neusticurus_ecpleopus
## 5 5
## Neusticurus_rudis Nothobachia_ablephara
## 5 2
## Pholidobolus_montium Placosoma_glabellum
## 5 5
## Procellosaurinus_erythrocercus Procellosaurinus_tetradactylus
## 5 5
## Proctoporus_bolivianus Proctoporus_simoterus
## 5 5
## Psilophthalmus_paeminosus Ptychoglossus_brevifrontalis
## 5 5
## Rhachisaurus_brachylepis Tretioscincus_agilis
## 4 5
## Vanzosaura_rubricauda Heloderma_suspectum
## 5 5
## Basiliscus_basiliscus Calotes_versicolor
## 5 5
## Dipsosaurus_dorsalis Enyaloides_laticeps
## 5 5
## Gambelia_wislizenii Leiocephalus_carinatus
## 5 5
## Leiolepis_belliana Phrynocephalus_versicolor
## 5 5
## Polychrus_marmoratus Urosaurus_ornatus
## 5 5
## Meroles_cuneirostris Podarcis_sicula
## 5 5
## Psammodromus_hispanicus Takydromus_smaragdinus
## 5 5
## Lanthanotus_borneensis Aprasia_pseudopulchella
## 5 0
## Aprasia_pulchella Aprasia_repens
## 0 0
## Aprasia_striolata Delma_australis
## 1 1
## Delma_borea Delma_fraseri
## 1 1
## Delma_grayi Delma_impar
## 1 1
## Delma_inornata Delma_molleri
## 1 1
## Delma_nasuta Delma_tincta
## 1 1
## Lialis_burtoni Lialis_jicari
## 0 0
## Pletholax_gracilis Pygopus_lepidopus
## 1 1
## Pygopus_nigriceps Acontias_litoralis
## 1 0
## Acontias_meleagris Acontias_percivali
## 0 0
## Acontophiops_lineatus Amphiglossus_astrolabi
## 0 5
## Amphiglossus_igneocaudatus Amphiglossus_intermedius
## 5 5
## Amphiglossus_macrocercus Amphiglossus_melanopleura
## 5 5
## Amphiglossus_melanurus Amphiglossus_mouroundavae
## 5 5
## Amphiglossus_ornaticeps Amphiglossus_punctatus
## 5 5
## Amphiglossus_stumpffi Amphiglossus_tsaratananensis
## 5 5
## Amphiglossus_waterloti Anomalopus_mackayi
## 5 1
## Anomalopus_swansoni Brachymeles_gracilis
## 0 5
## Brachymeles_talinis Calyptotis_scutirostrum
## 5 5
## Chalcides_chalcides Chalcides_mionecton
## 3 4
## Chalcides_ocellatus Chalcides_polylepis
## 5 5
## Coeranoscincus_reticulatus Ctenotus_leonhardii
## 3 5
## Ctenotus_pantherinus Ctenotus_robustus
## 5 5
## Egernia_whitii Eremiascincus_richardsonii
## 5 5
## Eugongylus_rufescens Eulamprus_amplus
## 5 5
## Eulamprus_murrayi Eulamprus_quoyii
## 5 5
## Eumeces_schneideri Feylinia_polylepis
## 5 0
## Glaphyromorphus_gracilipes Glaphyromorphus_isolepis
## 5 5
## Gnypetoscincus_queenslandiae Gongylomorphus_bojeri
## 5 5
## Hakaria_simonyi Hemiergis_peroni
## 5 4
## Janetaescincus_braueri Lamprolepis_smaragdina
## 5 5
## Lerista_bipes Lerista_bougainvillii
## 2 5
## Melanoseps_occidentalis Mesoscincus_managuae
## 0 5
## Mesoscincus_schwartzei Morethia_butleri
## 5 5
## Nangura_spinosa Notoscincus_ornatus
## 5 5
## Ophiomorus_punctatissimus Ophioscincus_ophioscincus
## 0 0
## Pamelaescincus_gardineri Paracontias_brocchii
## 5 0
## Paracontias_hildebrandti Paracontias_holomelas
## 0 0
## Plestiodon_egregius Plestiodon_elegans
## 5 5
## Plestiodon_fasciatus Plestiodon_inexpectatus
## 5 5
## Plestiodon_laticeps Plestiodon_longirostris
## 5 5
## Plestiodon_lynxe Plestiodon_obsoletus
## 5 5
## Plestiodon_reynoldsi Prasinohaema_virens
## 2 5
## Proscelotes_eggeli Pygomeles_braconnieri
## 5 1
## Saiphos_equalis Scelotes_anguineus
## 3 0
## Scelotes_arenicolus Scelotes_bipes
## 0 2
## Scelotes_caffer Scelotes_gronovii
## 3 1
## Scelotes_kasneri Scelotes_mirus
## 2 5
## Scelotes_sexlineatus Scincopus_fasciatus
## 2 5
## Scincus_scincus Sepsina_angolensis
## 5 3
## Sphenops_boulengeri Sphenops_sphenopsiformis
## 5 4
## Tiliqua_adelaidensis Tribolonotus_gracilis
## 5 5
## Typhlacontias_brevipes Typhlacontias_punctatissimus
## 0 0
## Typhlosaurus_caecus Typhlosaurus_lineatus
## 0 0
## Voeltzkowia_fierinensis Voeltzkowia_lineata
## 2 0
## Acrochordus_granulatus Agkistrodon_contortrix
## 0 0
## Anilius_scytale Calabaria_reinhardtii
## 0 1
## Candoia_carinata Casarea_dussumieri
## 0 0
## Causus_maculatus Cereberus_rhynchops
## 0 0
## Charina_bottae Chondropython_viridis
## 1 1
## Cylindrophis_maculatus Elaphe_guttata
## 0 0
## Epicrates_striatus Eryx_johni
## 1 1
## Exiliboa_plicata Farancia_abacura
## 1 0
## Leptotyphlops_humilis Loxocemus_bicolor
## 0 1
## Micrurus_fulvius Nerodia_sipedon
## 0 0
## Pareas_kuangtungensis Trachyboa_boulengeri
## 0 0
## Tropidophis_haitianus Typhlops_angolensis
## 1 0
## Ungaliophis_continentalis Uropeltis_ceylonicus
## 0 0
## Xenopeltis_unicolor Shinisaurus_crocodiluris
## 0 5
## Sphenodon_punctatus Ameiva_ameiva
## 5 5
## Aspidoscelis_sexlineatus Dicrodon_guttulatum
## 5 5
## Teius_teyou Varanus_griseus
## 4 5
## Lepidophyma_flavimaculatum Xantusia_vigilis
## 5 5
## Xenosaurus_grandis
## 5
## Levels: 0 1 2 3 4 5
Now, let's fit all our models in turn. Note that the exercise also requires that we take efforts to evaluate convergence of each model to it's genuine ML solution.
## ER model (1 parameter)
fitER<-fitDiscrete(squamates,hf.digits,model="ER",
control=list(niter=500))
## Warning in treedata(phy, dat): The following tips were not found in 'phy' and were dropped from 'data':
## Gonatodes_albogularis
## Lepidophyma_flavimaculatum
## Trachyboa_boulengeri
fitER
## GEIGER-fitted comparative model of discrete data
## fitted Q matrix:
## 0 1 2 3
## 0 -0.0045795845 0.0009159169 0.0009159169 0.0009159169
## 1 0.0009159169 -0.0045795845 0.0009159169 0.0009159169
## 2 0.0009159169 0.0009159169 -0.0045795845 0.0009159169
## 3 0.0009159169 0.0009159169 0.0009159169 -0.0045795845
## 4 0.0009159169 0.0009159169 0.0009159169 0.0009159169
## 5 0.0009159169 0.0009159169 0.0009159169 0.0009159169
## 4 5
## 0 0.0009159169 0.0009159169
## 1 0.0009159169 0.0009159169
## 2 0.0009159169 0.0009159169
## 3 0.0009159169 0.0009159169
## 4 -0.0045795845 0.0009159169
## 5 0.0009159169 -0.0045795845
##
## model summary:
## log-likelihood = -235.614300
## AIC = 473.228600
## AICc = 473.244225
## free parameters = 1
##
## Convergence diagnostics:
## optimization iterations = 500
## failed iterations = 0
## frequency of best fit = 1.00
##
## object summary:
## 'lik' -- likelihood function
## 'bnd' -- bounds for likelihood search
## 'res' -- optimization iteration summary
## 'opt' -- maximum likelihood parameter estimates
plot(fitER,signif=5)
title(main="ER model\nfor hindfoot digit loss/gain")
## ARD model (30 parameters!)
fitARD<-fitDiscrete(squamates,hf.digits,model="ARD",
control=list(niter=500))
## Warning in treedata(phy, dat): The following tips were not found in 'phy' and were dropped from 'data':
## Gonatodes_albogularis
## Lepidophyma_flavimaculatum
## Trachyboa_boulengeri
fitARD
## GEIGER-fitted comparative model of discrete data
## fitted Q matrix:
## 0 1 2 3
## 0 -3.264874e-17 3.175482e-17 7.264919e-23 4.566045e-45
## 1 1.505938e-02 -1.505938e-02 4.824941e-24 2.957816e-17
## 2 2.088159e-18 3.078712e-02 -3.078712e-02 4.372003e-136
## 3 2.681855e-02 8.229762e-22 9.008409e-03 -4.997548e-02
## 4 1.357745e-20 6.756395e-16 5.494325e-33 1.678902e-02
## 5 9.496594e-23 1.389924e-03 7.313331e-04 9.780807e-05
## 4 5
## 0 8.938406e-19 3.500030e-25
## 1 9.458480e-19 5.257957e-39
## 2 5.160562e-25 1.557695e-24
## 3 2.372426e-16 1.414853e-02
## 4 -5.293565e-02 3.614663e-02
## 5 2.063583e-03 -4.282648e-03
##
## model summary:
## log-likelihood = -191.369910
## AIC = 442.739820
## AICc = 450.933652
## free parameters = 30
##
## Convergence diagnostics:
## optimization iterations = 500
## failed iterations = 0
## frequency of best fit = 0.00
##
## object summary:
## 'lik' -- likelihood function
## 'bnd' -- bounds for likelihood search
## 'res' -- optimization iteration summary
## 'opt' -- maximum likelihood parameter estimates
plot(fitARD,signif=5,show.zeros=FALSE)
title(main="ARD model\nfor hindfoot digit loss/gain")
## ordered model (10 parameters)
model.ordered<-matrix(c(
0,1,0,0,0,0,
2,0,3,0,0,0,
0,4,0,5,0,0,
0,0,6,0,7,0,
0,0,0,8,0,9,
0,0,0,0,10,0),byrow=TRUE,6,6)
colnames(model.ordered)<-rownames(model.ordered)<-
levels(hf.digits)
model.ordered
## 0 1 2 3 4 5
## 0 0 1 0 0 0 0
## 1 2 0 3 0 0 0
## 2 0 4 0 5 0 0
## 3 0 0 6 0 7 0
## 4 0 0 0 8 0 9
## 5 0 0 0 0 10 0
fitOrdered<-fitDiscrete(squamates,hf.digits,
model=model.ordered,control=list(niter=500))
## 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(squamates, hf.digits, model = model.ordered, control = list(niter = 500)): Parameter estimates appear at bounds:
## q13
## q14
## q15
## q16
## q24
## q25
## q26
## q31
## q35
## q36
## q41
## q42
## q46
## q51
## q52
## q53
## q61
## q62
## q63
## q64
fitOrdered
## GEIGER-fitted comparative model of discrete data
## fitted Q matrix:
## 0 1 2 3 4
## 0 -5.339848e-17 5.339848e-17 0.000000e+00 0.00000000 0.00000000
## 1 2.106530e-02 -2.106530e-02 1.089522e-15 0.00000000 0.00000000
## 2 0.000000e+00 8.020167e-02 -3.064538e+00 2.98433653 0.00000000
## 3 0.000000e+00 0.000000e+00 4.223240e+00 -4.25269934 0.02945968
## 4 0.000000e+00 0.000000e+00 0.000000e+00 0.06492344 -4.28813991
## 5 0.000000e+00 0.000000e+00 0.000000e+00 0.00000000 0.24108513
## 5
## 0 0.0000000
## 1 0.0000000
## 2 0.0000000
## 3 0.0000000
## 4 4.2232165
## 5 -0.2410851
##
## model summary:
## log-likelihood = -203.914844
## AIC = 427.829687
## AICc = 428.720376
## free parameters = 10
##
## Convergence diagnostics:
## optimization iterations = 500
## failed iterations = 0
## frequency of best fit = 0.02
##
## object summary:
## 'lik' -- likelihood function
## 'bnd' -- bounds for likelihood search
## 'res' -- optimization iteration summary
## 'opt' -- maximum likelihood parameter estimates
plot(fitOrdered,signif=5,show.zeros=FALSE)
title(main="Ordered model\nfor hindfoot digit loss/gain")
## ordered directional model (5 parameters)
model.directional<-matrix(c(
0,0,0,0,0,0,
1,0,0,0,0,0,
0,2,0,0,0,0,
0,0,3,0,0,0,
0,0,0,4,0,0,
0,0,0,0,5,0),byrow=TRUE,6,6)
colnames(model.directional)<-rownames(model.directional)<-
levels(hf.digits)
model.directional
## 0 1 2 3 4 5
## 0 0 0 0 0 0 0
## 1 1 0 0 0 0 0
## 2 0 2 0 0 0 0
## 3 0 0 3 0 0 0
## 4 0 0 0 4 0 0
## 5 0 0 0 0 5 0
fitDirectional<-fitDiscrete(squamates,hf.digits,
model=model.directional,control=list(niter=500))
## 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(squamates, hf.digits, model = model.directional, : Parameter estimates appear at bounds:
## q12
## q13
## q14
## q15
## q16
## q23
## q24
## q25
## q26
## q31
## q34
## q35
## q36
## q41
## q42
## q45
## q46
## q51
## q52
## q53
## q56
## q61
## q62
## q63
## q64
fitDirectional
## GEIGER-fitted comparative model of discrete data
## fitted Q matrix:
## 0 1 2 3 4
## 0 0.00000000 0.00000000 0.0000000 0.0000000 0.000000000
## 1 0.02314705 -0.02314705 0.0000000 0.0000000 0.000000000
## 2 0.00000000 0.17375142 -0.1737514 0.0000000 0.000000000
## 3 0.00000000 0.00000000 0.2123592 -0.2123592 0.000000000
## 4 0.00000000 0.00000000 0.0000000 0.1118384 -0.111838407
## 5 0.00000000 0.00000000 0.0000000 0.0000000 0.004129312
## 5
## 0 0.000000000
## 1 0.000000000
## 2 0.000000000
## 3 0.000000000
## 4 0.000000000
## 5 -0.004129312
##
## model summary:
## log-likelihood = -211.796485
## AIC = 433.592970
## AICc = 433.831065
## free parameters = 5
##
## Convergence diagnostics:
## optimization iterations = 500
## failed iterations = 0
## frequency of best fit = 0.03
##
## object summary:
## 'lik' -- likelihood function
## 'bnd' -- bounds for likelihood search
## 'res' -- optimization iteration summary
## 'opt' -- maximum likelihood parameter estimates
plot(fitDirectional,signif=5,show.zeros=FALSE)
title(main="Directional model\nfor hindfoot digit loss (only)")
fitDiscrete's convergence diagnostics suggest that we have
converged on the genuine ML solutions; however we can also check our
likelihoods (model fits), knowing the in the case of nested models, the
model with more parameters (which has the second model as a special case)
should always have a likelihood score that is the same or better. In
particular, that means the "ARD" (which has as a special case
all other models), should have the best fit, but also that our ordered,
bi-directional model, which has the ordered, toe-loss, to-loss model as a
special case, should in turn have a better likelihood that that model.
Let's check (from least to most parameters, although not always nested):
setNames(c(logLik(fitER),logLik(fitDirectional),
logLik(fitOrdered),logLik(fitARD)),
c("ER","Directional","Ordered","ARD"))
## ER Directional Ordered ARD
## -235.6143 -211.7965 -203.9148 -191.3699
We should see that "ARD" has the highest likelihood, & that
"Ordered" has a higher likelihood than
"Directional".
To compare the models, we have multiple options. For instance, we can compare the models that are nested using likelihood-ratio tests. For instance:
models<-c("ER","Directional","Ordered","ARD")
LR.test<-matrix(NA,4,3,dimnames=list(paste(models,"-vs-ARD",sep=""),
c("LR","df","P(X2)")))
allfits<-list(fitER,fitDirectional,fitOrdered,fitARD)
for(i in 1:4){
LR.test[i,1]<--2*(logLik(allfits[[i]])-logLik(allfits[[4]]))
LR.test[i,2]<-allfits[[4]]$opt$k-allfits[[i]]$opt$k
LR.test[i,3]<-pchisq(LR.test[i,1],df=LR.test[i,2],
lower.tail=FALSE)
}
LR.test
## LR df P(X2)
## ER-vs-ARD 88.48878 29 6.207996e-08
## Directional-vs-ARD 40.85315 25 2.378670e-02
## Ordered-vs-ARD 25.08987 20 1.980146e-01
## ARD-vs-ARD 0.00000 0 1.000000e+00
However, we can also compared non-nested models. For instance:
aicc<-setNames(c(fitER$opt$aicc,fitDirectional$opt$aicc,
fitOrdered$opt$aicc,fitARD$opt$aicc),
c("ER","Directional","Ordered","ARD"))
aicc
## ER Directional Ordered ARD
## 473.2442 433.8311 428.7204 450.9337
and compute AIC weights, which give the weight of evidence of each model:
aic.w(aicc)
## ER Directional Ordered ARD
## 0.00000000 0.07206723 0.92791884 0.00001393
Hopefully these results show us the an ordered model fits best - but not a model in which digits can be lost but not re-acquired!
Written by Liam J. Revell 2016. Last updated 3 December 2016.