R.Version()$version.string; R.Version()$platform
## [1] "R version 4.4.1 (2024-06-14)"
## [1] "aarch64-apple-darwin20"

Loading packages

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'

Importing the data

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

Processing the data

#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 ...

Exploring Missing data

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

Evolutionary model fitting

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.

Missing data imputation using Rphylopars

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

Phylopars with full dataset

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),])

Log-shape ratios

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)

Phylomorphospace

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 phyllomorphospace

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))