The problem was basically required that we use the same data set for hindfoot
digit evolution from the excersise 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)
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 depending on the speed of your computer, you may need to use a lower niter (number of iterations) value than the default. In an empirical study we might set this higher to ensure convergence of each model to its genuine ML solution.
## ER model (1 parameter)
fitER<-fitDiscrete(squamates,hf.digits,model="ER",
control=list(niter=100))
## 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 = 100
## 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=100))
## 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 -7.599320e-17 7.559417e-17 2.824527e-44 3.990300e-19
## 1 1.424534e-02 -1.424534e-02 6.945340e-20 1.166216e-20
## 2 5.532104e-03 2.395673e-02 -3.388997e-02 3.366999e-18
## 3 8.582673e-42 4.744476e-28 1.000410e-02 -2.925922e+00
## 4 1.274108e-02 5.879914e-33 4.269602e-24 1.792738e+00
## 5 2.682728e-23 1.367758e-03 7.104514e-04 2.103568e-03
## 4 5
## 0 1.986762e-25 3.529999e-30
## 1 2.711761e-21 3.811274e-24
## 2 2.597660e-47 4.401145e-03
## 3 2.842752e+00 7.316567e-02
## 4 -1.805479e+00 3.653471e-20
## 5 5.212739e-19 -4.181778e-03
##
## model summary:
## log-likelihood = -194.746243
## AIC = 449.492486
## AICc = 457.686318
## free parameters = 30
##
## Convergence diagnostics:
## optimization iterations = 100
## failed iterations = 0
## frequency of best fit = 0.05
##
## 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=100))
## 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 = 100)): 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 = 100
## 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=100))
## 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.02314691 -0.02314691 0.0000000 0.0000000 0.000000000
## 2 0.00000000 0.17375095 -0.1737510 0.0000000 0.000000000
## 3 0.00000000 0.00000000 0.2123602 -0.2123602 0.000000000
## 4 0.00000000 0.00000000 0.0000000 0.1118386 -0.111838628
## 5 0.00000000 0.00000000 0.0000000 0.0000000 0.004129307
## 5
## 0 0.000000000
## 1 0.000000000
## 2 0.000000000
## 3 0.000000000
## 4 0.000000000
## 5 -0.004129307
##
## model summary:
## log-likelihood = -211.796485
## AIC = 433.592970
## AICc = 433.831065
## free parameters = 5
##
## Convergence diagnostics:
## optimization iterations = 100
## failed iterations = 0
## frequency of best fit = 0.07
##
## 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)")
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 81.73611 29 6.434794e-07
## Directional-vs-ARD 34.10048 25 1.057774e-01
## Ordered-vs-ARD 18.33720 20 5.652057e-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 457.6863
and compute AIC weights, which give the weight of evidence of each model:
aic.w(aicc)
## ER Directional Ordered ARD
## 0.00000000 0.07206820 0.92793132 0.00000048
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 2017. Last updated 2 August 2017.