R.Version()$version.string; R.Version()$platform
## [1] "R version 4.4.1 (2024-06-14)"
## [1] "aarch64-apple-darwin20"
library(geomorph)
packageVersion("geomorph")
## [1] '4.0.8'
library(treeplyr)
packageVersion("treeplyr")
## [1] '0.1.11'
library(Rphylopars)
packageVersion("Rphylopars")
## [1] '0.3.10'
library(mvMORPH)
packageVersion("mvMORPH")
## [1] '1.2.1'
library(naniar)
packageVersion("naniar")
## [1] '1.1.0'
library(phytools)
packageVersion("phytools")
## [1] '2.3.0'
library(tidyverse)
packageVersion("tidyverse")
## [1] '2.0.0'
library(psych)
packageVersion("psych")
## [1] '2.4.6.26'
morph_data <- read.table('pca_caltomontana.txt', row.names = 1)
glimpse(morph_data)
## Rows: 14
## Columns: 16
## $ BL <dbl> 6.89, 8.80, 9.65, 9.84, 2.73, 6.96, 6.81, 5.30, 6.60, 6.75, 5…
## $ TAL <dbl> 10.51, 11.20, 13.43, 13.45, 5.80, 13.39, 11.53, 8.79, 11.90, …
## $ TL <dbl> 17.40, 18.60, 22.55, 22.98, 8.52, 18.49, 18.35, 14.17, 18.50,…
## $ ED <dbl> 0.95, 1.10, 0.84, 0.91, NA, 0.99, 0.90, 0.86, 0.90, 0.98, 0.8…
## $ BH <dbl> 4.27, NA, 5.87, 5.92, NA, NA, 4.35, 3.60, 3.90, NA, 3.30, NA,…
## $ MTH <dbl> 4.95, 6.40, 7.08, 6.97, 1.45, 3.77, 4.30, 4.60, 5.10, 3.68, 3…
## $ TMH <dbl> 1.71, 2.20, 2.75, 2.83, NA, NA, 2.16, 1.72, 1.80, 1.90, 1.40,…
## $ DFH <dbl> 1.63, NA, 2.44, 2.56, NA, NA, NA, 1.37, 1.40, NA, 1.32, NA, N…
## $ VFH <dbl> 2.00, NA, 2.72, 2.93, NA, NA, NA, 1.72, 2.00, NA, 1.48, NA, N…
## $ IOD <dbl> 5.26, 6.20, 6.83, 6.81, 0.93, 5.75, NA, 4.24, 5.50, 3.88, 2.9…
## $ BW <dbl> 5.70, NA, 6.78, 6.79, NA, 5.47, 5.06, 4.37, 5.50, NA, 4.00, N…
## $ TMW <dbl> NA, 2.10, 1.44, 1.70, 0.45, NA, NA, NA, 1.40, NA, NA, NA, 0.4…
## $ ESD <dbl> 2.46, NA, 4.39, 4.43, NA, 3.69, NA, NA, 3.50, NA, 1.60, NA, N…
## $ Subgenus <chr> "Chiasmocleis", "Chiasmocleis", "Chiasmocleis", "Chiasmocleis…
## $ Biome <chr> "AF", "CH_CE", "AM", "AF", "AM", "AF", "AF", "AF", "AF", "AM"…
## $ Stage <int> NA, NA, 36, NA, NA, NA, NA, NA, NA, NA, NA, NA, 18, NA
head(morph_data)
## BL TAL TL ED BH MTH TMH DFH VFH IOD BW TMW ESD
## alagoana 6.89 10.51 17.40 0.95 4.27 4.95 1.71 1.63 2.00 5.26 5.70 NA 2.46
## albopunctata 8.80 11.20 18.60 1.10 NA 6.40 2.20 NA NA 6.20 NA 2.10 NA
## anatipes 9.65 13.43 22.55 0.84 5.87 7.08 2.75 2.44 2.72 6.83 6.78 1.44 4.39
## altomontana 9.84 13.45 22.98 0.91 5.92 6.97 2.83 2.56 2.93 6.81 6.79 1.70 4.43
## hudsoni 2.73 5.80 8.52 NA NA 1.45 NA NA NA 0.93 NA 0.45 NA
## lacrimae 6.96 13.39 18.49 0.99 NA 3.77 NA NA NA 5.75 5.47 NA 3.69
## Subgenus Biome Stage
## alagoana Chiasmocleis AF NA
## albopunctata Chiasmocleis CH_CE NA
## anatipes Chiasmocleis AM 36
## altomontana Chiasmocleis AF NA
## hudsoni Syncope AM NA
## lacrimae Chiasmocleis AF NA
morph_data$Biome <- factor(morph_data$Biome)
biome <- factor(c("AF", "AF", "AF", "CH_CE", "AF","AF", "AF","AM","AM","AM","AM","AM","AM","AM"))
# Importing the phylogeny
phylo_chiasm <- read.nexus("Fig_SI_2Chiasmocleis_MrBayes2 - editada.tre")
is.ultrametric(phylo_chiasm)
## [1] FALSE
#match phylo data
tree_data <- make.treedata(phylo_chiasm, morph_data)
summary(tree_data)
## A treeplyr treedata object
## The dataset contains 16 traits
## Continuous traits: BL TAL TL ED BH MTH TMH DFH VFH IOD BW TMW ESD Subgenus Biome Stage
## Discrete traits:
## The following traits have missing values: ED, BH, MTH, TMH, DFH, VFH, IOD, BW, TMW, ESD, Stage
## These taxa were dropped from the tree:
## These taxa were dropped from the data: Xenopus_laevis, Amietia_angolensis, Tomopterna_tuberculosa, Lithobates_spp, Polypedates_leucomystax, Ptychadena_anchietae, Breviceps_mossambicus, Callulina_kisiwamsitu, Spelaeophryne_methneri, Hemisus_marmoratus, Altigius_alios, Hamptophryne_boliviana, Arcovomer_passarellii, Dermatonotus_muelleri, Elachistocleis_bicolor, Elachistocleis_helianneae, Elachistocleis_ovalis, Elachistocleis_pearsei, Elachistocleis_panamensis, Gastrophryne_carolinensis, Gastrophryne_olivacea, Gastrophryne_elegans, Hypopachus_barberi, Hypopachus_variolosus, Hypopachus_pictiventris, Hypopachus_ustus, Dasypops_schirchi, Myersiella_microps, Stereocyclops_histrio, Stereocyclops_incrassatus, alagoanus_2683, cordeiroi_CFBH15784_IuriC_BA_Ilh, cordeiroi_CTA1935_BA_Itubera, cordeiroi_IuriD_BA_Igrapiuna, cordeiroi_MTR22122_BA_WenceslauG, cordeiroi_MTR22123_BA_WenceslauG, cordeiroi_CTMZ7023_BA_Amargosa, cordeiroi_PEU137_BA_Jaguaripe, cordeiroi_PEU146_BA_Jaguaripe, schubarti_CFBH9027_ES_Linhares_V, schubarti_CFBH9060_ES_Linhares_P, schubarti_CTMZ6906_CTA1937_ES_Li, schubarti_MTR12266_ES_Linhares_G, schubarti_CTMZ6915_ES_Linhares, schubarti_LGA3600_ES_SantaTeresa, schubarti_CTA1860_ES_Cariacica, schubarti_CTA1887_ES_Conceicaoda, schubarti_LGA2630_LGA2606_ES_Pin, schubarti_CTMZ6910_ES_Linhares, schubarti_MTR17524_MG_Marlieria, schubarti_MTR17755_MG_Marlieria, crucis_CFBH11434_BA_Una, crucis_MTR16118_BA_Camacan, crucis_CTMZ6901_BA_Ilheus, crucis_MTR16070_BA_Camacan, crucis_MTR16087_BA_Camacan, crucis_CFBH11437_BA_Una, crucis_MTR16106_BA_Camacan, crucis_CTMZ6900_BA_Ilheus, crucis_MTR16198_BA_Camacan, crucis_CTMZ6898_BA_Jussari, capixaba_CFBH22990, capixaba_CTA1875_ES_Serra, capixaba_CTA1870_ES_Cariacica, capixaba_CTA1874_ES_Serra, capixaba_MTR12076_ES_Linhares_Va, capixaba_CTMZ6908_ES_Linhares_Go, capixaba_MTR12297_ES_Linhares_Go, capixaba_CTA1861_ES_Cariacica, capixaba_CTMZ6939_ES_Guarapari, lacrimae_MNRJ66494_ES_MimosodoSu, quilombola_CFBH1437, quilombola_CFBH9082_ES_Linhares, quilombola_LGA3267_ES_Pinheiros, quilombola_CFBH1438, quilombola_CTA1909_ES_Conceicaod, veracruz_CFBH15818_BA_PortoSegur, veracruz_CTMZ6924_BA_Trancoso, veracruz_CTMZ6935_BA_Trancoso, veracruz_CTMZ6930_BA_Trancoso, veracruz_CTMZ6923_BA_Trancoso, veracruz_CTMZ6927_BA_Trancoso, lacrimae_47477_RJ_Macae, lacrimae_48415, lacrimae_CFBH73_SP_Ubatuba, lacrimae_CFBH76_SP_Ubatuba, lacrimae_MNRJ49302_RJ_Cachoeiras, lacrimae_MNRJ60744_RJ_DuquedeCax, lacrimae_A12281, lacrimae_A12282, lacrimae_CFBH7361_SP_IlhaBela, lacrimae_CTMZ6946_SP_Bertioga, lacrimae_CTRN175_RJ_AngradosReis, lacrimae_ZUEC_DCC3433, albopunctata_C572_SP_Itirapina, albopunctata_C621_SP_RioClaro, albopunctata_CHUNB44451_MG_Burit, mehelyi_MAP680, mehelyi_MAP681, albopunctata_CHUNB45602_Caseara_TO, albopunctata_CTMZ5444_Tesouro_MT, albopunctata_CTMZ5445_Tesouro_MT, albopunctata_CTMZ5489_MT_Tangara, albopunctata_CTMZ6872_Pocone_MT, albopunctata_JMP218_Bolivia, albopunctata_LECOH0001_Argentina, albopunctata_CHUNB51830_MA_Carol, albopunctata_MPEG29414_TO_Aragua, albopunctata_VLF076_Jardim_MS, albopunctata_IIBPH74_Paraguay, altomontani_CTMZ1172_SP_Bananal, mantiqueira_CTMZ6902_SP_Camposdo, mantiqueira_CTMZ6894_SP_Camposdo, mantiqueira_CFBH13670_SP_Desenga, mantiqueira_CTMZ6891_MG_SerradoB, mantiqueira_CTMZ6893_MG_SerradoB, mantiqueira_UFMGt1810_UFMGt1801, mantiqueira_UFMGt1815_MG_OuroBra, leucosticta_CTMZ2484_CTMZ2485_SP, leucosticta_CTMZ6943_CTMZ6942_SP, leucosticta_CTMZ6944_SP_Piedade, leucosticta_CTMZ2164_SP_Jacupira, leucosticta_CTMZ2166_SP_Jacupira, leucosticta_CFBH8594_CFBH8595_SP, leucosticta_CFBH8622_SP_Cananeia, anatipes_MNCN27460_Peru_Loreto, anatipes_QCAZ51041_Ecuador_Napo, anatipes_QCAZ51042_Ecuador_Napo, anatipes_QCAZ44341_Ecuador_Orell, devriesi_MHNSM21540, anatipes_QCAZ44342_Ecuador_Orell, papachibe_MPEG27788_PA_Barcarena, papachibe_MPEG30683, papachibe_MPEG30684_PA_Paragomin, ventrimaculata_JMP1944, ventrimaculata_MHNSM21539, royi_KU215540_Peru_MadredeDios, royi_KU215542_Peru_MadredeDios, royi_ROM40139_Peru_MadredeDios, ventrimaculata_MNCN10054, avilapiresae_CTMZ6650_MT_Paranai, avilapiresae_CTMZ6974_PA_NovoPro, avilapiresae_MPEG27768_AM_Maues, avilapiresae_MPEG27769_AM_Maues, avilapiresae_HJ445_RO_Abuna, avilapiresae_MTR12730_AM_SaoSeba, avilapiresae_HJ587_RO_Mutum, avilapiresae_MPEG22787_PA_Portel, avilapiresae_MPEG28121_PA_Caraja, avilapiresae_MPEG18571_PA_Itaitu, shudikarensis_MW5645_FrenchGuyan, shudikarensis_AMNHFS20018_others, shudikarensis_INPA31274_AM_Tamag, shudikarensis_INPA31275_AM_Tamag, shudikarensis_MVZ247574_Guyana, shudikarensis_MNRJ53427_MPEG2832, shudikarensis_JIW458_Suriname, sp_QCAZ41441, sp_QCAZ51773, antenori_MTR28373_AC_Divisor, antenori_MTR28416_AC_Divisor, tingomaria_MHNC5461, urubamba_MHNC10266, urubamba_MHNC10304, antenori_QCAZ38719_Ecuador_Pasta, antenori_QCAZ56275, antenori_QCAZ23824_Ecuador_Orell, antenori_QCAZ55408_Ecuador_Tambococha, tridactyla_QCAZ55408_Ecuador_Tambococha, magnova_NMPGV71196, carvalhoi_31_Peru_Loreto, carvalhoi_63_Peru_Loreto, carvalhoi_MNCN27392_Peru_Loreto, carvalhoi_MNCN27503_Peru_Loreto, carvalhoi_MNCN3713, carvalhoi_MNCN27507, carvalhoi_MNCN27505, tridactyla_MNCN27490, carvalhoi_MNCN27506, tridactyla_MNCN46870, carvalhoi_MNCN26570_Peru_Loreto, carvalhoi_MNCN26618, carvalhoi_MNCN46994, tridactyla_MNCN46820, tridactyla_JMP1948, tridactyla_MNCN46711, tridactyla_MNCN46605, tridactyla_MNCN46828, tridactyla_JMP2231, tridactyla_JMP2057, sp_MNCN47179, haddadi_28bm_FrenchGuiana_MontBa, haddadi_MNHN20110139_FrenchGuian, haddadi_PG445_FrenchGuiana_MontB, haddadi_MNHN20110140_FrenchGuian, haddadi_PG446_FrenchGuiana_MontB, hudsoni_INPA31276_RO, hudsoni_HJ55_RO_Jirau, hudsoni_MAD116_Guyana_Iwokrama, hudsoni_MPEG23206_PA_Faro, hudsoni_FPR96_MPEG27763_AM_Maues, hudsoni_MTR13195_AM_Pacamiri, hudsoni_MTR12728_AM_IgarapeAcu, hudsoni_MTR12941_AM_SaoSebastiao, hudsoni_MTR13197_AM_Pacamiri, hudsoni_MPEG23285_PA_RioXingu, hudsoni_MPEG24527_PA_Portel, hudsoni_T275_PA_Tapirape, hudsoni_T81_PA_Tapirape, hudsoni_JMP2286_Colombia_Amazona, hudsoni_MPEG18456, hudsoni_MPEG18900, hudsoni_MPEG18545_PA_Itaituba, hudsoni_MPEG18547, bassleri_ALP14940_AM, supercilialba_MNCN62999, supercilialba_MNCN63003, supercilialba_MNCN63001, supercilialba_MNCN63007, bassleri_CTMZ6764_MT_Paranaita, bassleri_CTMZ6877_MT_Aripuana, bassleri_CTMZ6882_MT_Apiacas, bassleri_HJ112_HJ232_RO_Jirau, bassleri_HJ596_RO_Mutum, bassleri_MPEG28326_PA_Juriti, bassleri_MPEG27764_AM_Maues, bassleri_MPEG27765_AM_Maues, bassleri_HJ650_RO_Caicara, bassleri_MNCN47166_Colombia_Vaup, gnoma_CFBH15820, Ctenophryne_aterrima, Ctenophryne_equatorialis, Ctenophryne_geayi, Ctenophryne_barbatula, Anodonthohyla_spp, Platypelis_grandis, Plethodontohyla_spp, Rhombophryne_testudo, Dyscophus_antongilii, Dyscophus_guineti, Glyphoglossus_molossus, Microhyla_heymonsi, Microhyla_ornata, Micryletta_inornata, Kaloula_picta, Kaloula_pulchra, Phrynomantis_microps, Ramanella_variegata, Uperodon_spp, Scaphiophryne_calcarata, Scaphiophryne_madagascariensis, Kalophrynus_interlineatus, Callulops_spp, Metamagnusia_slateri, Xenorhina_obesa, Sphenophryne_spp, Cophixalus_spp, Oreophryne_spp, Otophryne_robusta, Otophryne_steyermarki, Synapturanus_spp, Hoplophryne_rogersi
## $phy
##
## Phylogenetic tree with 14 tips and 13 internal nodes.
##
## Tip labels:
## alagoana, schubarti, lacrimae, albopunctata, altomontana, mantiqueira, ...
##
## Rooted; includes branch lengths.
##
## $dat
## # A tibble: 14 × 16
## BL TAL TL ED BH MTH TMH DFH VFH IOD BW TMW ESD
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 6.89 10.5 17.4 0.95 4.27 4.95 1.71 1.63 2 5.26 5.7 NA 2.46
## 2 6.6 11.9 18.5 0.9 3.9 5.1 1.8 1.4 2 5.5 5.5 1.4 3.5
## 3 6.96 13.4 18.5 0.99 NA 3.77 NA NA NA 5.75 5.47 NA 3.69
## 4 8.8 11.2 18.6 1.1 NA 6.4 2.2 NA NA 6.2 NA 2.1 NA
## 5 9.84 13.4 23.0 0.91 5.92 6.97 2.83 2.56 2.93 6.81 6.79 1.7 4.43
## 6 5.3 8.79 14.2 0.86 3.6 4.6 1.72 1.37 1.72 4.24 4.37 NA NA
## 7 6.81 11.5 18.4 0.9 4.35 4.3 2.16 NA NA NA 5.06 NA NA
## 8 9.65 13.4 22.6 0.84 5.87 7.08 2.75 2.44 2.72 6.83 6.78 1.44 4.39
## 9 7 10 17 NA NA NA NA NA NA NA NA NA NA
## 10 6.75 8.9 15.7 0.98 NA 3.68 1.9 NA NA 3.88 NA NA NA
## 11 3.41 8.31 11.7 0.37 2.06 2.47 0.91 NA NA NA 2.66 0.48 NA
## 12 5.8 10.8 16.6 0.84 3.3 3.7 1.4 1.32 1.48 2.9 4 NA 1.6
## 13 2.88 5.49 8.37 0.45 1.82 1.13 0.7 NA NA 0.64 2 NA NA
## 14 2.73 5.8 8.52 NA NA 1.45 NA NA NA 0.93 NA 0.45 NA
## # ℹ 3 more variables: Subgenus <fct>, Biome <fct>, Stage <int>
tree_data$phy
##
## Phylogenetic tree with 14 tips and 13 internal nodes.
##
## Tip labels:
## alagoana, schubarti, lacrimae, albopunctata, altomontana, mantiqueira, ...
##
## Rooted; includes branch lengths.
plotTree(tree_data$phy)
str(tree_data$dat)
## tibble [14 × 16] (S3: tbl_df/tbl/data.frame)
## $ BL : num [1:14] 6.89 6.6 6.96 8.8 9.84 5.3 6.81 9.65 7 6.75 ...
## $ TAL : num [1:14] 10.5 11.9 13.4 11.2 13.4 ...
## $ TL : num [1:14] 17.4 18.5 18.5 18.6 23 ...
## $ ED : num [1:14] 0.95 0.9 0.99 1.1 0.91 0.86 0.9 0.84 NA 0.98 ...
## $ BH : num [1:14] 4.27 3.9 NA NA 5.92 3.6 4.35 5.87 NA NA ...
## $ MTH : num [1:14] 4.95 5.1 3.77 6.4 6.97 4.6 4.3 7.08 NA 3.68 ...
## $ TMH : num [1:14] 1.71 1.8 NA 2.2 2.83 1.72 2.16 2.75 NA 1.9 ...
## $ DFH : num [1:14] 1.63 1.4 NA NA 2.56 1.37 NA 2.44 NA NA ...
## $ VFH : num [1:14] 2 2 NA NA 2.93 1.72 NA 2.72 NA NA ...
## $ IOD : num [1:14] 5.26 5.5 5.75 6.2 6.81 4.24 NA 6.83 NA 3.88 ...
## $ BW : num [1:14] 5.7 5.5 5.47 NA 6.79 4.37 5.06 6.78 NA NA ...
## $ TMW : num [1:14] NA 1.4 NA 2.1 1.7 NA NA 1.44 NA NA ...
## $ ESD : num [1:14] 2.46 3.5 3.69 NA 4.43 NA NA 4.39 NA NA ...
## $ Subgenus: Factor w/ 2 levels "Chiasmocleis",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Biome : Factor w/ 3 levels "AF","AM","CH_CE": 1 1 1 3 1 1 1 2 2 2 ...
## $ Stage : int [1:14] NA NA NA NA NA NA NA 36 NA NA ...
only_morph <- tree_data$dat %>%
select(BL:ESD)
nearly_complete_morph <- tree_data$dat %>%
select(BL:TMH, IOD, BW)
vis_miss(only_morph)#13 variables | 27.5% missing
vis_miss(nearly_complete_morph)#9 variables | 14.3% missing overall
bm <- mvBM(tree_data$phy, nearly_complete_morph, model = "BM1")
## species in the matrix are assumed to be in the same order as in the phylogeny, otherwise specify rownames of 'data'
##
## convergence of the optimizer has not been reached, try simpler model
## unreliable solution has been reached, check hessian eigenvectors or try simpler model
##
## -- Summary results for multiple rate BM1 model --
## LogLikelihood: -66.64756
## AIC: 241.2951
## AICc: 353.3706
## 54 parameters
##
## Estimated rate matrix
## ______________________
## BL TAL TL ED BH MTH
## BL 0.124865848 0.077545880 0.070666033 0.0048441534 0.034291783 0.037758287
## TAL 0.077545880 0.130759729 0.114479757 0.0029027838 0.028718470 0.029645921
## TL 0.070666033 0.114479757 0.381013944 0.0029519766 0.078947185 0.099835835
## ED 0.004844153 0.002902784 0.002951977 0.0013068885 0.001037524 0.001367513
## BH 0.034291783 0.028718470 0.078947185 0.0010375242 0.023713003 0.026164974
## MTH 0.037758287 0.029645921 0.099835835 0.0013675129 0.026164974 0.048127152
## TMH 0.014495315 0.015454489 0.041953710 0.0003909764 0.010937150 0.011792067
## IOD 0.032680985 0.048239901 0.104923921 0.0016024752 0.025319757 0.034720740
## BW 0.033532795 0.035553940 0.090502782 0.0017408108 0.023623557 0.030140673
## TMH IOD BW
## BL 0.0144953148 0.032680985 0.033532795
## TAL 0.0154544888 0.048239901 0.035553940
## TL 0.0419537105 0.104923921 0.090502782
## ED 0.0003909764 0.001602475 0.001740811
## BH 0.0109371497 0.025319757 0.023623557
## MTH 0.0117920666 0.034720740 0.030140673
## TMH 0.0063151415 0.013254381 0.010739972
## IOD 0.0132543809 0.039339368 0.030759554
## BW 0.0107399719 0.030759554 0.027387189
##
## Estimated root state
## ______________________
## BL TAL TL ED BH MTH TMH IOD
## theta: 5.798242 9.344886 14.91867 0.77963 3.494701 3.728357 1.589682 3.5984
## BW
## theta: 4.221421
eb <- mvEB(tree_data$phy, nearly_complete_morph)
## species in the matrix are assumed to be in the same order as in the phylogeny, otherwise specify rownames of 'data'
## No lower bound provided. Use of default setting " -0.2952306 "
## maximum limit iteration has been reached, please consider increase maxit
## unreliable solution has been reached, check hessian eigenvectors or try simpler model
##
## -- Summary results for Early Burst or ACDC model --
## LogLikelihood: -106.0753
## AIC: 322.1506
## AICc: 440.6121
## 55 parameters
## Rate change:
## ______________________
## [1] 0
##
## Estimated rate matrix
## ______________________
## BL TAL TL ED BH
## BL 0.141766582 0.086224164 0.07683940 0.0264262235 0.1068051038
## TAL 0.086224164 0.130865312 0.11036246 0.0145375088 0.0603950028
## TL 0.076839397 0.110362461 0.30208651 -0.0141148675 0.0127294129
## ED 0.026426223 0.014537509 -0.01411487 0.0104983282 0.0256241715
## BH 0.106805104 0.060395003 0.01272941 0.0256241715 0.0906566652
## MTH 0.056502516 0.052863056 0.04557451 0.0044322657 0.0488790314
## TMH 0.003309587 0.007582614 0.03577143 -0.0055359317 -0.0007944269
## IOD -0.050479828 -0.041144823 0.19137947 -0.0351972725 -0.0881591732
## BW 0.036866532 0.019442771 0.08684366 -0.0005339748 0.0131612219
## MTH TMH IOD BW
## BL 0.056502516 0.0033095870 -0.050479828 0.0368665315
## TAL 0.052863056 0.0075826141 -0.041144823 0.0194427710
## TL 0.045574509 0.0357714325 0.191379468 0.0868436569
## ED 0.004432266 -0.0055359317 -0.035197272 -0.0005339748
## BH 0.048879031 -0.0007944269 -0.088159173 0.0131612219
## MTH 0.100567585 0.0270421592 -0.078467069 0.0027268417
## TMH 0.027042159 0.0143333217 0.007458263 0.0050687303
## IOD -0.078467069 0.0074582630 0.343702105 0.0789436309
## BW 0.002726842 0.0050687303 0.078943631 0.0381457800
##
## Estimated root states
## ______________________
## BL TAL TL ED BH MTH TMH IOD
## theta 5.798242 9.344886 14.91867 0.796422 3.665043 3.721601 1.547495 3.73808
## BW
## theta 4.282402
ou <- mvOU(tree_data$phy, nearly_complete_morph, model = "OU1", param = list(method="rpf") )
## species in the matrix are assumed to be in the same order as in the phylogeny, otherwise specify rownames of 'data'
## successful convergence of the optimizer
## unreliable solution has been reached, check hessian eigenvectors or try simpler model
##
## -- Summary results --
## LogLikelihood: -14.36973
## AIC: 226.7395
## AICc: 2701.739
## 99 parameters
##
## Estimated theta values
## ______________________
## BL TAL TL ED BH MTH TMH IOD
## OU1 5.79212 9.556684 15.16282 0.7670049 3.5444 3.889587 1.570293 3.754317
## BW
## OU1 4.344308
##
## ML alpha values
## ______________________
## BL TAL TL ED BH MTH
## BL 0.393805220 0.07446245 -0.13809035 -0.007885661 -0.092419756 -0.06496476
## TAL 0.074462448 0.26379836 -0.21376872 0.016545133 0.064779360 0.03654304
## TL -0.138090354 -0.21376872 0.25540206 -0.021205996 -0.044885231 -0.03397251
## ED -0.007885661 0.01654513 -0.02120600 0.431444905 -0.009196911 0.05747504
## BH -0.092419756 0.06477936 -0.04488523 -0.009196911 0.372859755 -0.02427106
## MTH -0.064964760 0.03654304 -0.03397251 0.057475042 -0.024271065 0.22600322
## TMH -0.062854513 0.02790761 -0.02422277 0.003282572 -0.069178023 -0.01036645
## IOD -0.041858776 -0.02042531 -0.01231794 -0.015327825 -0.029950549 -0.04617092
## BW -0.048254796 0.00691237 -0.03809273 -0.045366651 -0.050544791 -0.02038833
## TMH IOD BW
## BL -0.062854513 -0.04185878 -0.04825480
## TAL 0.027907613 -0.02042531 0.00691237
## TL -0.024222765 -0.01231794 -0.03809273
## ED 0.003282572 -0.01532782 -0.04536665
## BH -0.069178023 -0.02995055 -0.05054479
## MTH -0.010366448 -0.04617092 -0.02038833
## TMH 0.431084697 -0.05001726 0.06837168
## IOD -0.050017257 0.27802009 -0.13181903
## BW 0.068371684 -0.13181903 0.38761488
##
## ML sigma values
## ______________________
## BL TAL TL ED BH MTH
## BL 0.24302768 0.08519265 0.22466682 -0.006850490 0.098010788 0.04285332
## TAL 0.08519265 0.24827921 0.28786978 -0.023473135 0.014168550 0.03235283
## TL 0.22466682 0.28786978 0.90645930 -0.036690491 0.204871642 0.30764605
## ED -0.00685049 -0.02347313 -0.03669049 0.007136362 -0.006211753 -0.01158154
## BH 0.09801079 0.01416855 0.20487164 -0.006211753 0.080086098 0.08760113
## MTH 0.04285332 0.03235283 0.30764605 -0.011581540 0.087601127 0.26097718
## TMH 0.02495125 0.02067983 0.11124768 -0.006868312 0.032875774 0.04272427
## IOD 0.01294291 0.09811343 0.25801942 -0.010667386 0.042067851 0.14107626
## BW 0.06740360 0.05559559 0.23986585 -0.004009099 0.063080358 0.11569625
## TMH IOD BW
## BL 0.024951250 0.01294291 0.067403602
## TAL 0.020679831 0.09811343 0.055595587
## TL 0.111247676 0.25801942 0.239865853
## ED -0.006868312 -0.01066739 -0.004009099
## BH 0.032875774 0.04206785 0.063080358
## MTH 0.042724269 0.14107626 0.115696251
## TMH 0.021411533 0.03280402 0.027008128
## IOD 0.032804017 0.20500991 0.125667007
## BW 0.027008128 0.12566701 0.105576633
mvMORPH::aicw(list(bm, eb, ou), aicc=TRUE)
## -- Akaike weights --
## Rank AIC diff wi AICw
## BM1 default 1 1 353 0.0 1 1
## ACDC 2 2 441 87.2 0 0
## OU1 3 3 2702 2348.4 0 0
BM is the best model.
species <- rownames(morph_data)
tadpoles_reord <- cbind(species, nearly_complete_morph)
head(tadpoles_reord)
## species BL TAL TL ED BH MTH TMH IOD BW
## 1 alagoana 6.89 10.51 17.40 0.95 4.27 4.95 1.71 5.26 5.70
## 2 albopunctata 6.60 11.90 18.50 0.90 3.90 5.10 1.80 5.50 5.50
## 3 anatipes 6.96 13.39 18.49 0.99 NA 3.77 NA 5.75 5.47
## 4 altomontana 8.80 11.20 18.60 1.10 NA 6.40 2.20 6.20 NA
## 5 hudsoni 9.84 13.45 22.98 0.91 5.92 6.97 2.83 6.81 6.79
## 6 lacrimae 5.30 8.79 14.17 0.86 3.60 4.60 1.72 4.24 4.37
impt_data <- phylopars(tadpoles_reord, tree_data$phy, model = "BM")
impt_data_phypars <-as.data.frame(impt_data$anc_recon[-c(15:27),])
head(impt_data_phypars)
## BL TAL TL ED BH MTH TMH IOD
## alagoana 6.89 10.51 17.40 0.950000 4.270000 4.950000 1.710000 5.260000
## schubarti 7.00 10.00 17.00 0.952175 4.521072 5.112425 1.946353 5.070827
## lacrimae 5.30 8.79 14.17 0.860000 3.600000 4.600000 1.720000 4.240000
## albopunctata 6.60 11.90 18.50 0.900000 3.900000 5.100000 1.800000 5.500000
## altomontana 8.80 11.20 18.60 1.100000 5.268059 6.400000 2.200000 6.200000
## mantiqueira 9.65 13.43 22.55 0.840000 5.870000 7.080000 2.750000 6.830000
## BW
## alagoana 5.700000
## schubarti 5.476825
## lacrimae 4.370000
## albopunctata 5.500000
## altomontana 5.821802
## mantiqueira 6.780000
Let’s look at the data:
GGally::ggpairs(impt_data_phypars)
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
new_complete <- rownames_to_column(morph_data, var = "species")
impt_data_new <- phylopars(new_complete[,1:14], tree_data$phy, model = "BM")
impt_data_phypars_new <-as.data.frame(impt_data_new$anc_recon[-c(15:27),])
Computing the geometric mean for obtaining size for the reduced dataset (variables < 57% missing)
size <- geometric.mean(t(impt_data_phypars))
Doing the same with the full dataset
size2 <- geometric.mean(t(impt_data_phypars_new))
##computing the log shape ratios
Reduced dataset
LSR <- log(impt_data_phypars/ size)
LSR[1:4,1:4]
## BL TAL TL ED
## alagoana 0.3661081 0.7883642 1.292507 -1.615256
## schubarti 0.3722310 0.7289060 1.259534 -1.622686
## lacrimae 0.2665661 0.7724740 1.249986 -1.551964
## albopunctata 0.3133454 0.9028142 1.344046 -1.679085
GGally::ggpairs(LSR)
Full dataset
LSR2 <- log(impt_data_phypars_new/ size2)
phylomorpho<-gm.prcomp(phy=tree_data$phy, LSR)#PCA
phylomorpho$rotation[,1:2]
## Comp1 Comp2
## BL 0.09410904 -0.06940087
## TAL 0.17319816 -0.38730930
## TL 0.09747095 -0.28036400
## ED 0.13869545 -0.34521475
## BH -0.19680186 -0.04663830
## MTH 0.04835839 0.30814735
## TMH -0.85947307 0.15306714
## IOD 0.36494697 0.72367569
## BW 0.13949597 -0.05596296
phylomorpho$x#species scores
## Comp1 Comp2 Comp3 Comp4 Comp5
## alagoana 0.04989487 0.24644226 -0.089303996 0.01083508 -0.014790659
## schubarti -0.09907199 0.27316945 -0.116066758 -0.00209459 0.030219424
## lacrimae -0.13027083 0.27092178 -0.128279521 0.16736399 0.051534442
## albopunctata 0.05233370 0.25836155 0.041119870 0.12101382 -0.085354349
## altomontana -0.07174865 0.36126560 -0.134098575 -0.03348880 0.120962592
## mantiqueira -0.20187600 0.44537359 0.215577458 -0.12636799 0.014573265
## leucosticta -0.21452409 0.11578701 -0.014733830 0.01809835 -0.138139838
## anatipes 1.75365054 -0.14270227 0.005575633 -0.04660236 -0.003669904
## ventrimaculata -0.04140366 -0.21582677 0.033209711 0.08282697 0.022632649
## shudikarensis -0.23277846 0.04545131 -0.266509188 -0.09624347 -0.051799473
## antenori -0.37041438 -0.99625923 -0.050709409 -0.13628266 -0.000486511
## carvalhoi -0.14341435 -0.32932527 0.479022390 0.07081653 0.023558230
## magnova -0.13559890 -0.74013917 -0.129857234 0.08460951 0.025496265
## hudsoni -0.21477781 0.40748016 0.155053446 -0.11448437 0.005263868
## Comp6 Comp7 Comp8
## alagoana 0.069579757 0.007866697 0.0004854996
## schubarti 0.047190411 0.000558463 0.0080201394
## lacrimae 0.024488109 -0.010449794 0.0034612375
## albopunctata 0.011237415 0.023878480 -0.0124611504
## altomontana -0.049480152 -0.003522770 -0.0167910483
## mantiqueira 0.009020579 -0.005393212 -0.0020605167
## leucosticta -0.041648613 -0.021733142 -0.0068016507
## anatipes -0.002548512 -0.003772060 0.0006248414
## ventrimaculata -0.087045108 0.018501798 0.0116566400
## shudikarensis -0.017123040 -0.003617091 0.0136894242
## antenori 0.024643986 0.009266596 -0.0078844483
## carvalhoi 0.015371428 -0.009918836 0.0039104196
## magnova 0.007103168 -0.009606702 -0.0012781728
## hudsoni -0.010789429 0.007941575 0.0054287854
Loadings: PC1 : TMH neg, IOD pos, PC2: IOD pos TAL neg
plot(phylomorpho, pch = 16, col = biome, cex = 2, phylo = TRUE,
phylo.par = list(edge.color = "grey60", edge.width = 1.5,
tip.txt.cex = 1, node.labels = F,
anc.states = F))