rm(list=ls())
setwd("~/Dropbox/FCU/Teaching/Mentoring/2022Spring-Current_Professor/2022Fall/YourueiMin/rcode_dcj/empricaldataanalysis/matchedtreetraitTony/BodyMass")
source("~/Dropbox/FCU/Teaching/Mentoring/2022Spring-Current_Professor/2022Fall/YourueiMin/rcode_dcj/empricaldataanalysis/mainPhyarima.R")
## Loading required package: ape
## Loading required package: phytools
## Loading required package: maps
##
## Attaching package: 'extraDistr'
## The following object is masked from 'package:geiger':
##
## dtpois
##
## Attaching package: 'LaplacesDemon'
## The following objects are masked from 'package:extraDistr':
##
## dbern, dcat, ddirichlet, dgpd, dgpois, dinvchisq, dinvgamma,
## dlaplace, dpareto, pbern, plaplace, ppareto, qbern, qcat, qlaplace,
## qpareto, rbern, rcat, rdirichlet, rgpd, rinvchisq, rinvgamma,
## rlaplace, rpareto
## The following object is masked from 'package:geiger':
##
## is.constrained
## Loading required package: emmeans
##
## Attaching package: 'RRphylo'
## The following object is masked from 'package:geiger':
##
## tips
## The following object is masked from 'package:phytools':
##
## node.paths
## Loading required package: MASS
## Loading required package: minpack.lm
## Loading required package: rgl
## Loading required package: robustbase
## Loading required package: Matrix
###########################
##### diurnal_subtree #####
###########################
diurnal_subtree<-read.tree("diurnal_subtree.nwk")
tiplabel<-diurnal_subtree$tip.label
tiplabel
## [1] "Propithecus_deckenii" "Propithecus_verreauxi"
## [3] "Propithecus_candidus" "Propithecus_coquereli"
## [5] "Propithecus_tattersalli" "Propithecus_edwardsi"
## [7] "Propithecus_diadema" "Propithecus_perrieri"
## [9] "Varecia_variegata" "Varecia_rubra"
## [11] "Hapalemur_occidentalis" "Hapalemur_aureus"
## [13] "Lemur_catta" "Eulemur_mongoz"
## [15] "Pithecia_vanzolinii" "Pithecia_hirsuta"
## [17] "Cacajao_rubicundus" "Cacajao_calvus"
## [19] "Chiropotes_albinasus" "Chiropotes_chiropotes"
## [21] "Cebus_olivaceus" "Cebus_kaapori"
## [23] "Saimiri_ustus" "Ateles_marginatus"
## [25] "Brachyteles_arachnoides" "Lagothrix_lugens"
## [27] "Alouatta_macconnelli" "Papio_kindae"
## [29] "Mandrillus_sphinx" "Cercocebus_agilis"
## [31] "Macaca_sylvanus" "Semnopithecus_entellus"
## [33] "Trachypithecus_johnii" "Colobus_angolensis"
length(tiplabel)
## [1] 34
diurnal_BodyMass<-read.csv("diurnal_BodyMass.csv")
diurnal_BodyMass
## Species BodyMass_kg BodyMassMale_kg BodyMassFemale_kg
## 1 Alouatta_macconnelli 6.3700 7.200 5.550
## 2 Ateles_marginatus 6.9600 10.400 5.800
## 3 Brachyteles_arachnoides 8.8400 9.610 8.070
## 4 Cacajao_calvus 3.1700 3.450 2.880
## 5 Cebus_kaapori 2.4250 3.000 2.400
## 6 Cebus_olivaceus 2.8800 3.240 2.520
## 7 Cercocebus_agilis 7.5800 9.500 5.660
## 8 Chiropotes_albinasus 2.8200 3.150 2.490
## 9 Chiropotes_chiropotes 2.7400 2.900 2.580
## 10 Colobus_angolensis 8.6000 9.800 7.400
## 11 Eulemur_mongoz 1.3750 1.410 1.560
## 12 Hapalemur_aureus 1.5150 1.520 1.390
## 13 Hapalemur_occidentalis 0.9315 0.846 1.180
## 14 Lemur_catta 2.2100 2.210 2.210
## 15 Macaca_sylvanus 13.5000 16.000 11.000
## 16 Mandrillus_sphinx 23.6000 34.400 12.800
## 17 Papio_kindae 13.0000 16.000 9.900
## 18 Pithecia_hirsuta NA NA NA
## 19 Pithecia_vanzolinii NA NA NA
## 20 Propithecus_candidus 5.6850 NA NA
## 21 Propithecus_coquereli 3.7850 3.703 3.756
## 22 Propithecus_deckenii 2.7800 2.930 2.630
## 23 Propithecus_diadema 6.1050 6.050 6.210
## 24 Propithecus_edwardsi 5.7400 5.590 5.895
## 25 Propithecus_perrieri 4.4100 4.220 4.440
## 26 Propithecus_tattersalli 3.5000 3.390 3.590
## 27 Propithecus_verreauxi 2.9100 2.885 2.870
## 28 Saimiri_ustus 0.8600 0.921 0.799
## 29 Semnopithecus_entellus 11.4400 13.000 9.890
## 30 Varecia_rubra 3.2650 3.550 3.470
## 31 Varecia_variegata 3.5500 3.630 3.520
length(diurnal_BodyMass$Species)
## [1] 31
same<-intersect(tiplabel,diurnal_BodyMass$Species)
diurnal_BodyMass$Species%in%same
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [31] TRUE
# Get the names of the tips to drop
tips_to_drop <- setdiff(tiplabel, diurnal_BodyMass$Species)
tips_to_drop<-c(tips_to_drop,"Pithecia_hirsuta","Pithecia_vanzolinii","Propithecus_candidus")
# Drop the tips
diurnal_subtree_with_bodymass <- drop.tip(diurnal_subtree, tips_to_drop)
plot(diurnal_subtree_with_bodymass)
diurnal_BodyMass <- diurnal_BodyMass[!(diurnal_BodyMass$Species %in% c("Pithecia_hirsuta", "Pithecia_vanzolinii", "Propithecus_candidus")), ]
dim(diurnal_BodyMass)
## [1] 28 4
y<-diurnal_BodyMass$BodyMassMale_kg
names(y)<-diurnal_BodyMass$Species
#######################################
############ diurnal Analysis #######
#######################################
mfit<-modelfit(trait=y, tree=diurnal_subtree_with_bodymass)
diurnal.empmodelout<-empmodel(mfit)
diurnal.empmodelout$parammleset
## phi phi2 theta theta2 sigsq
## AR1 0.42588485 NA NA NA 0.2916463
## AR2 0.56653671 -0.27656899 NA NA 0.2959332
## ARMA1 -0.16885716 NA 0.7704224 NA 0.2885269
## ARMA21 -0.03872185 -0.13013469 0.7704221 NA 0.2885269
## ARMA22 0.02187839 -0.19074593 -0.2113302 0.98175575 0.2885287
## Weighted Estimate 0.28545674 -0.08250971 0.2170826 0.02945267 0.2916770
diurnal.empmodelout$fitset
## modelset paranumset loglikset aicset aicweightset rankset
## 1 AR1 2 -42.55 89.10 0.44 1
## 2 AR2 3 -42.13 90.25 0.24 2
## 3 ARMA1 3 -42.26 90.53 0.21 3
## 4 ARMA21 4 -42.26 92.53 0.08 4
## 5 ARMA22 5 -42.26 94.53 0.03 5
#empmodelout$parammleset[1,][c("phi","sigsq")]
xtable(diurnal.empmodelout$parammleset,digits=3)
## % latex table generated in R 4.4.0 by xtable 1.8-4 package
## % Sat May 4 14:12:45 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrr}
## \hline
## & phi & phi2 & theta & theta2 & sigsq \\
## \hline
## AR1 & 0.426 & & & & 0.292 \\
## AR2 & 0.567 & -0.277 & & & 0.296 \\
## ARMA1 & -0.169 & & 0.770 & & 0.289 \\
## ARMA21 & -0.039 & -0.130 & 0.770 & & 0.289 \\
## ARMA22 & 0.022 & -0.191 & -0.211 & 0.982 & 0.289 \\
## Weighted Estimate & 0.285 & -0.083 & 0.217 & 0.029 & 0.292 \\
## \hline
## \end{tabular}
## \end{table}
xtable(diurnal.empmodelout$fitset,digits=3)
## % latex table generated in R 4.4.0 by xtable 1.8-4 package
## % Sat May 4 14:12:45 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rlrrrrr}
## \hline
## & modelset & paranumset & loglikset & aicset & aicweightset & rankset \\
## \hline
## 1 & AR1 & 2.000 & -42.550 & 89.100 & 0.440 & 1.000 \\
## 2 & AR2 & 3.000 & -42.130 & 90.250 & 0.240 & 2.000 \\
## 3 & ARMA1 & 3.000 & -42.260 & 90.530 & 0.210 & 3.000 \\
## 4 & ARMA21 & 4.000 & -42.260 & 92.530 & 0.080 & 4.000 \\
## 5 & ARMA22 & 5.000 & -42.260 & 94.530 & 0.030 & 5.000 \\
## \hline
## \end{tabular}
## \end{table}
## add up test.phi here
diurnal_betas.rrphylo<-compute_betas_rrphylo(diurnal_subtree_with_bodymass, y)$betas.rrphylo
test.phi(param=diurnal.empmodelout$parammleset["AR1",][c("phi","sigsq")], betas.rrphylo=diurnal_betas.rrphylo,tree=diurnal_subtree_with_bodymass)
## $z
## phi
## 4.303365
##
## $pval
## pval
## 1.682235e-05
###########################
##### nocturnal_subtree ###
###########################
nocturnal_subtree<-read.tree("nocturnal_subtree.nwk")
tiplabel<-nocturnal_subtree$tip.label
tiplabel
## [1] "Perodicticus_potto" "Arctocebus_calabarensis"
## [3] "Arctocebus_aureus" "Nycticebus_pygmaeus"
## [5] "Nycticebus_bengalensis" "Nycticebus_coucang"
## [7] "Nycticebus_javanicus" "Nycticebus_menagensis"
## [9] "Loris_lydekkerianus" "Euoticus_elegantulus"
## [11] "Galagoides_demidoff" "Galagoides_thomasi"
## [13] "Galagoides_zanzibaricus" "Galagoides_cocos"
## [15] "Galago_granti" "Galagoides_orinus"
## [17] "Galagoides_rondoensis" "Galago_senegalensis"
## [19] "Galago_gallarum" "Galago_moholi"
## [21] "Galago_matschiei" "Galago_alleni"
## [23] "Galago_gabonensis" "Otolemur_monteiri"
## [25] "Otolemur_garnettii" "Otolemur_crassicaudatus"
## [27] "Microcebus_myoxinus" "Cheirogaleus_major"
## [29] "Cheirogaleus_medius" "Lepilemur_otto"
## [31] "Lepilemur_grewcockorum" "Lepilemur_manasamody"
## [33] "Lepilemur_edwardsi" "Lepilemur_microdon"
## [35] "Lepilemur_tymerlachsoni" "Lepilemur_milanoii"
## [37] "Lepilemur_ankaranensis" "Lepilemur_ahmansonorum"
## [39] "Lepilemur_sahamalazensis" "Lepilemur_mittermeieri"
## [41] "Lepilemur_dorsalis" "Lepilemur_septentrionalis"
## [43] "Lepilemur_hubbardorum" "Lepilemur_ruficaudatus"
## [45] "Lepilemur_randrianasoloi" "Lepilemur_aeeclis"
## [47] "Lepilemur_petteri" "Lepilemur_leucopus"
## [49] "Lepilemur_betsileo" "Lepilemur_jamesorum"
## [51] "Lepilemur_mustelinus" "Lepilemur_fleuretae"
## [53] "Lepilemur_wrightae" "Lepilemur_seali"
length(tiplabel)
## [1] 54
nocturnal_BodyMass<-read.csv("nocturnal_BodyMass.csv")
length(nocturnal_BodyMass$Species)
## [1] 47
same<-intersect(tiplabel,nocturnal_BodyMass$Species)
nocturnal_BodyMass$Species%in%same
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [46] TRUE TRUE
# Get the names of the tips to drop
tips_to_drop <- setdiff(tiplabel, nocturnal_BodyMass$Species)
tips_to_drop
## [1] "Galagoides_demidoff" "Galago_granti" "Galago_alleni"
## [4] "Galago_gabonensis" "Otolemur_monteiri" "Lepilemur_manasamody"
## [7] "Lepilemur_hubbardorum"
# Drop the tips
nocturnal_subtree_with_bodymass <- drop.tip(nocturnal_subtree, tips_to_drop)
plot(nocturnal_subtree_with_bodymass)
nocturnal_BodyMass
## Species BodyMass_kg BodyMassMale_kg BodyMassFemale_kg
## 1 Arctocebus_aureus 0.2100000 NA NA
## 2 Arctocebus_calabarensis 0.3100000 0.3120 0.3060
## 3 Cheirogaleus_major 0.5050000 0.5740 0.4430
## 4 Cheirogaleus_medius 0.1766667 0.1640 0.1560
## 5 Euoticus_elegantulus 0.2700000 0.2870 0.2610
## 6 Galago_gallarum 0.2000000 0.2350 0.2000
## 7 Galago_matschiei 0.2100000 0.2070 0.2120
## 8 Galago_moholi 0.2000000 0.2140 0.1940
## 9 Galago_senegalensis 0.2500000 0.2710 0.2250
## 10 Galagoides_cocos 0.1400000 0.1500 0.1300
## 11 Galagoides_orinus 0.0900000 NA NA
## 12 Galagoides_rondoensis 0.6750000 0.6900 0.6600
## 13 Galagoides_thomasi 0.1200000 0.1030 0.1300
## 14 Galagoides_zanzibaricus 0.1400000 0.1500 0.1370
## 15 Lepilemur_aeeclis 0.8600000 0.8680 0.9090
## 16 Lepilemur_ahmansonorum 0.5950000 NA NA
## 17 Lepilemur_ankaranensis 0.7600000 NA NA
## 18 Lepilemur_betsileo 1.1400000 NA NA
## 19 Lepilemur_dorsalis 0.7950000 NA NA
## 20 Lepilemur_edwardsi 0.9750000 0.9280 0.9340
## 21 Lepilemur_fleuretae 0.8900000 NA NA
## 22 Lepilemur_grewcockorum 0.8200000 NA NA
## 23 Lepilemur_jamesorum 0.8250000 NA NA
## 24 Lepilemur_leucopus 0.5850000 0.6170 0.5940
## 25 Lepilemur_microdon 1.0350000 0.9700 NA
## 26 Lepilemur_milanoii 0.7150000 NA NA
## 27 Lepilemur_mittermeieri 0.8600000 NA NA
## 28 Lepilemur_mustelinus 0.9800000 NA NA
## 29 Lepilemur_otto 0.9390000 NA NA
## 30 Lepilemur_petteri 0.6200000 0.6550 0.5850
## 31 Lepilemur_randrianasoloi 0.8550000 NA NA
## 32 Lepilemur_ruficaudatus 0.7800000 0.7290 0.7185
## 33 Lepilemur_sahamalazensis 0.6900000 0.7100 0.6780
## 34 Lepilemur_seali 0.9500000 NA NA
## 35 Lepilemur_septentrionalis 0.6275000 NA NA
## 36 Lepilemur_tymerlachsoni 0.8800000 NA NA
## 37 Lepilemur_wrightae 1.4350000 1.8000 1.8500
## 38 Loris_lydekkerianus 0.2700000 0.2640 0.2690
## 39 Microcebus_myoxinus 0.0400000 0.0307 0.0307
## 40 Nycticebus_bengalensis 1.2600000 1.1300 1.4000
## 41 Nycticebus_coucang 0.6500000 0.6790 0.6260
## 42 Nycticebus_javanicus 0.7300000 0.6860 0.6260
## 43 Nycticebus_menagensis 0.2850000 0.2850 NA
## 44 Nycticebus_pygmaeus 0.4200000 0.4620 0.3760
## 45 Otolemur_crassicaudatus 1.1500000 1.1900 1.1100
## 46 Otolemur_garnettii 0.7600000 0.7940 0.7340
## 47 Perodicticus_potto 0.8300000 0.8300 0.8360
nocturnal_subtree_with_bodymass
##
## Phylogenetic tree with 47 tips and 46 internal nodes.
##
## Tip labels:
## Perodicticus_potto, Arctocebus_calabarensis, Arctocebus_aureus, Nycticebus_pygmaeus, Nycticebus_bengalensis, Nycticebus_coucang, ...
## Node labels:
## '204', '82', '43', '13', '14', '34', ...
##
## Rooted; includes branch lengths.
y<-nocturnal_BodyMass$BodyMassMale_kg
names(y)<-nocturnal_BodyMass$Species
tips_to_drop<-names(y)[is.na(y)]
tips_to_drop
## [1] "Arctocebus_aureus" "Galagoides_orinus"
## [3] "Lepilemur_ahmansonorum" "Lepilemur_ankaranensis"
## [5] "Lepilemur_betsileo" "Lepilemur_dorsalis"
## [7] "Lepilemur_fleuretae" "Lepilemur_grewcockorum"
## [9] "Lepilemur_jamesorum" "Lepilemur_milanoii"
## [11] "Lepilemur_mittermeieri" "Lepilemur_mustelinus"
## [13] "Lepilemur_otto" "Lepilemur_randrianasoloi"
## [15] "Lepilemur_seali" "Lepilemur_septentrionalis"
## [17] "Lepilemur_tymerlachsoni"
# Drop the tips
nocturnal_subtree_with_bodymass <- drop.tip(nocturnal_subtree_with_bodymass, tips_to_drop)
nocturnal_subtree_with_bodymass
##
## Phylogenetic tree with 30 tips and 29 internal nodes.
##
## Tip labels:
## Perodicticus_potto, Arctocebus_calabarensis, Nycticebus_pygmaeus, Nycticebus_bengalensis, Nycticebus_coucang, Nycticebus_javanicus, ...
## Node labels:
## '204', '82', '43', '13', '34', '35', ...
##
## Rooted; includes branch lengths.
length(y)
## [1] 47
y
## Arctocebus_aureus Arctocebus_calabarensis Cheirogaleus_major
## NA 0.3120 0.5740
## Cheirogaleus_medius Euoticus_elegantulus Galago_gallarum
## 0.1640 0.2870 0.2350
## Galago_matschiei Galago_moholi Galago_senegalensis
## 0.2070 0.2140 0.2710
## Galagoides_cocos Galagoides_orinus Galagoides_rondoensis
## 0.1500 NA 0.6900
## Galagoides_thomasi Galagoides_zanzibaricus Lepilemur_aeeclis
## 0.1030 0.1500 0.8680
## Lepilemur_ahmansonorum Lepilemur_ankaranensis Lepilemur_betsileo
## NA NA NA
## Lepilemur_dorsalis Lepilemur_edwardsi Lepilemur_fleuretae
## NA 0.9280 NA
## Lepilemur_grewcockorum Lepilemur_jamesorum Lepilemur_leucopus
## NA NA 0.6170
## Lepilemur_microdon Lepilemur_milanoii Lepilemur_mittermeieri
## 0.9700 NA NA
## Lepilemur_mustelinus Lepilemur_otto Lepilemur_petteri
## NA NA 0.6550
## Lepilemur_randrianasoloi Lepilemur_ruficaudatus Lepilemur_sahamalazensis
## NA 0.7290 0.7100
## Lepilemur_seali Lepilemur_septentrionalis Lepilemur_tymerlachsoni
## NA NA NA
## Lepilemur_wrightae Loris_lydekkerianus Microcebus_myoxinus
## 1.8000 0.2640 0.0307
## Nycticebus_bengalensis Nycticebus_coucang Nycticebus_javanicus
## 1.1300 0.6790 0.6860
## Nycticebus_menagensis Nycticebus_pygmaeus Otolemur_crassicaudatus
## 0.2850 0.4620 1.1900
## Otolemur_garnettii Perodicticus_potto
## 0.7940 0.8300
y<-y[!is.na(y)]
y
## Arctocebus_calabarensis Cheirogaleus_major Cheirogaleus_medius
## 0.3120 0.5740 0.1640
## Euoticus_elegantulus Galago_gallarum Galago_matschiei
## 0.2870 0.2350 0.2070
## Galago_moholi Galago_senegalensis Galagoides_cocos
## 0.2140 0.2710 0.1500
## Galagoides_rondoensis Galagoides_thomasi Galagoides_zanzibaricus
## 0.6900 0.1030 0.1500
## Lepilemur_aeeclis Lepilemur_edwardsi Lepilemur_leucopus
## 0.8680 0.9280 0.6170
## Lepilemur_microdon Lepilemur_petteri Lepilemur_ruficaudatus
## 0.9700 0.6550 0.7290
## Lepilemur_sahamalazensis Lepilemur_wrightae Loris_lydekkerianus
## 0.7100 1.8000 0.2640
## Microcebus_myoxinus Nycticebus_bengalensis Nycticebus_coucang
## 0.0307 1.1300 0.6790
## Nycticebus_javanicus Nycticebus_menagensis Nycticebus_pygmaeus
## 0.6860 0.2850 0.4620
## Otolemur_crassicaudatus Otolemur_garnettii Perodicticus_potto
## 1.1900 0.7940 0.8300
length(y)
## [1] 30
#######################################
############ nocturnal Analysis #######
#######################################
#compute_betas_rrphylo(nocturnal_subtree_with_bodymass, y)$betas.rrphylo
mfit<-modelfit(trait=y, tree=nocturnal_subtree_with_bodymass)
nocturnal.empmodelout<-empmodel(mfit)
nocturnal.empmodelout$parammleset
## phi phi2 theta theta2 sigsq
## AR1 -0.088092596 NA NA NA 0.002745380
## AR2 -0.095383322 -0.088129315 NA NA 0.002832202
## ARMA1 0.789382942 NA -0.9999971 NA 0.002879928
## ARMA21 0.802306192 -0.017950138 -0.9999971 NA 0.002884236
## ARMA22 0.106274777 0.544479601 -0.3255466 -0.6744097 0.003016103
## Weighted Estimate 0.008598689 -0.003182384 -0.1099997 0.0000000 0.002762915
nocturnal.empmodelout$fitset
## modelset paranumset loglikset aicset aicweightset rankset
## 1 AR1 2 87.33 -170.66 0.86 1
## 2 AR2 3 84.92 -163.83 0.03 4
## 3 ARMA1 3 85.97 -165.93 0.08 2
## 4 ARMA21 4 85.92 -163.84 0.03 3
## 5 ARMA22 5 84.53 -159.06 0.00 5
xtable(nocturnal.empmodelout$parammleset,digits=3)
## % latex table generated in R 4.4.0 by xtable 1.8-4 package
## % Sat May 4 14:12:48 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrr}
## \hline
## & phi & phi2 & theta & theta2 & sigsq \\
## \hline
## AR1 & -0.088 & & & & 0.003 \\
## AR2 & -0.095 & -0.088 & & & 0.003 \\
## ARMA1 & 0.789 & & -1.000 & & 0.003 \\
## ARMA21 & 0.802 & -0.018 & -1.000 & & 0.003 \\
## ARMA22 & 0.106 & 0.544 & -0.326 & -0.674 & 0.003 \\
## Weighted Estimate & 0.009 & -0.003 & -0.110 & 0.000 & 0.003 \\
## \hline
## \end{tabular}
## \end{table}
xtable(nocturnal.empmodelout$fitset,digits=3)
## % latex table generated in R 4.4.0 by xtable 1.8-4 package
## % Sat May 4 14:12:48 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rlrrrrr}
## \hline
## & modelset & paranumset & loglikset & aicset & aicweightset & rankset \\
## \hline
## 1 & AR1 & 2.000 & 87.330 & -170.660 & 0.860 & 1.000 \\
## 2 & AR2 & 3.000 & 84.920 & -163.830 & 0.030 & 4.000 \\
## 3 & ARMA1 & 3.000 & 85.970 & -165.930 & 0.080 & 2.000 \\
## 4 & ARMA21 & 4.000 & 85.920 & -163.840 & 0.030 & 3.000 \\
## 5 & ARMA22 & 5.000 & 84.530 & -159.060 & 0.000 & 5.000 \\
## \hline
## \end{tabular}
## \end{table}
## add up test.phi here
nocturnal_betas.rrphylo<-compute_betas_rrphylo(nocturnal_subtree_with_bodymass, y)$betas.rrphylo
test.phi(param=nocturnal.empmodelout$parammleset["AR1",][c("phi","sigsq")], betas.rrphylo=nocturnal_betas.rrphylo,tree=nocturnal_subtree_with_bodymass)
## $z
## phi
## -0.9013365
##
## $pval
## pval
## 0.3674094
####################################
#########Combined tree##############
####################################
tree_with_bodymass<-bind.tree(diurnal_subtree_with_bodymass, nocturnal_subtree_with_bodymass, where = "root", position = 0, interactive = FALSE)
tree_with_bodymass<-chronopl(tree_with_bodymass, 0, age.min = 1, age.max = NULL,node = "root", S = 1, tol = 1e-8,CV = FALSE, eval.max = 500, iter.max = 500)
## Warning in chronopl(tree_with_bodymass, 0, age.min = 1, age.max = NULL, : at least one internal branch is of length zero:
## it was collapsed and some nodes have been deleted.
plot(tree_with_bodymass, show.tip.label = FALSE)
BodyMass<-rbind(diurnal_BodyMass,nocturnal_BodyMass)
y<-BodyMass$BodyMassMale_kg
names(y)<-BodyMass$Species
y
## Alouatta_macconnelli Ateles_marginatus Brachyteles_arachnoides
## 7.2000 10.4000 9.6100
## Cacajao_calvus Cebus_kaapori Cebus_olivaceus
## 3.4500 3.0000 3.2400
## Cercocebus_agilis Chiropotes_albinasus Chiropotes_chiropotes
## 9.5000 3.1500 2.9000
## Colobus_angolensis Eulemur_mongoz Hapalemur_aureus
## 9.8000 1.4100 1.5200
## Hapalemur_occidentalis Lemur_catta Macaca_sylvanus
## 0.8460 2.2100 16.0000
## Mandrillus_sphinx Papio_kindae Propithecus_coquereli
## 34.4000 16.0000 3.7030
## Propithecus_deckenii Propithecus_diadema Propithecus_edwardsi
## 2.9300 6.0500 5.5900
## Propithecus_perrieri Propithecus_tattersalli Propithecus_verreauxi
## 4.2200 3.3900 2.8850
## Saimiri_ustus Semnopithecus_entellus Varecia_rubra
## 0.9210 13.0000 3.5500
## Varecia_variegata Arctocebus_aureus Arctocebus_calabarensis
## 3.6300 NA 0.3120
## Cheirogaleus_major Cheirogaleus_medius Euoticus_elegantulus
## 0.5740 0.1640 0.2870
## Galago_gallarum Galago_matschiei Galago_moholi
## 0.2350 0.2070 0.2140
## Galago_senegalensis Galagoides_cocos Galagoides_orinus
## 0.2710 0.1500 NA
## Galagoides_rondoensis Galagoides_thomasi Galagoides_zanzibaricus
## 0.6900 0.1030 0.1500
## Lepilemur_aeeclis Lepilemur_ahmansonorum Lepilemur_ankaranensis
## 0.8680 NA NA
## Lepilemur_betsileo Lepilemur_dorsalis Lepilemur_edwardsi
## NA NA 0.9280
## Lepilemur_fleuretae Lepilemur_grewcockorum Lepilemur_jamesorum
## NA NA NA
## Lepilemur_leucopus Lepilemur_microdon Lepilemur_milanoii
## 0.6170 0.9700 NA
## Lepilemur_mittermeieri Lepilemur_mustelinus Lepilemur_otto
## NA NA NA
## Lepilemur_petteri Lepilemur_randrianasoloi Lepilemur_ruficaudatus
## 0.6550 NA 0.7290
## Lepilemur_sahamalazensis Lepilemur_seali Lepilemur_septentrionalis
## 0.7100 NA NA
## Lepilemur_tymerlachsoni Lepilemur_wrightae Loris_lydekkerianus
## NA 1.8000 0.2640
## Microcebus_myoxinus Nycticebus_bengalensis Nycticebus_coucang
## 0.0307 1.1300 0.6790
## Nycticebus_javanicus Nycticebus_menagensis Nycticebus_pygmaeus
## 0.6860 0.2850 0.4620
## Otolemur_crassicaudatus Otolemur_garnettii Perodicticus_potto
## 1.1900 0.7940 0.8300
length(y)
## [1] 75
names(y)%in% tree_with_bodymass$tip.label
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [13] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [25] TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [37] TRUE TRUE FALSE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE TRUE
## [49] FALSE FALSE FALSE TRUE TRUE FALSE FALSE FALSE FALSE TRUE FALSE TRUE
## [61] TRUE FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [73] TRUE TRUE TRUE
tree_with_bodymass$tip.label %in%names(y)
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
#######################################
############ Full Analysis ##########
#######################################
#compute_betas_rrphylo(tree_with_bodymass, y)$betas.rrphylo
mfit<-modelfit(trait=y, tree=tree_with_bodymass)
comb.empmodelout<-empmodel(mfit)
comb.empmodelout$parammleset
## phi phi2 theta theta2 sigsq
## AR1 0.4486049 NA NA NA 2.828582
## AR2 0.3622290 0.1532215 NA NA 2.708466
## ARMA1 0.8009025 NA -0.46561495 NA 2.633089
## ARMA21 0.6988761 0.1020271 -0.46561862 NA 2.633084
## ARMA22 0.3582694 0.4426323 -0.08695956 -0.37865766 2.633076
## Weighted Estimate 0.5489580 0.1050271 -0.21300555 -0.01514631 2.673933
comb.empmodelout$fitset
## modelset paranumset loglikset aicset aicweightset rankset
## 1 AR1 2 -209.39 422.79 0.02 5
## 2 AR2 3 -205.13 416.27 0.49 1
## 3 ARMA1 3 -205.53 417.05 0.33 2
## 4 ARMA21 4 -205.53 419.05 0.12 3
## 5 ARMA22 5 -205.53 421.05 0.04 4
xtable(comb.empmodelout$parammleset,digits=3)
## % latex table generated in R 4.4.0 by xtable 1.8-4 package
## % Sat May 4 14:12:52 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrr}
## \hline
## & phi & phi2 & theta & theta2 & sigsq \\
## \hline
## AR1 & 0.449 & & & & 2.829 \\
## AR2 & 0.362 & 0.153 & & & 2.708 \\
## ARMA1 & 0.801 & & -0.466 & & 2.633 \\
## ARMA21 & 0.699 & 0.102 & -0.466 & & 2.633 \\
## ARMA22 & 0.358 & 0.443 & -0.087 & -0.379 & 2.633 \\
## Weighted Estimate & 0.549 & 0.105 & -0.213 & -0.015 & 2.674 \\
## \hline
## \end{tabular}
## \end{table}
xtable(comb.empmodelout$fitset,digits=3)
## % latex table generated in R 4.4.0 by xtable 1.8-4 package
## % Sat May 4 14:12:52 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rlrrrrr}
## \hline
## & modelset & paranumset & loglikset & aicset & aicweightset & rankset \\
## \hline
## 1 & AR1 & 2.000 & -209.390 & 422.790 & 0.020 & 5.000 \\
## 2 & AR2 & 3.000 & -205.130 & 416.270 & 0.490 & 1.000 \\
## 3 & ARMA1 & 3.000 & -205.530 & 417.050 & 0.330 & 2.000 \\
## 4 & ARMA21 & 4.000 & -205.530 & 419.050 & 0.120 & 3.000 \\
## 5 & ARMA22 & 5.000 & -205.530 & 421.050 & 0.040 & 4.000 \\
## \hline
## \end{tabular}
## \end{table}
## add up test.phi here
comb_betas.rrphylo<-compute_betas_rrphylo(tree_with_bodymass, y)$betas.rrphylo
test.phi(param=comb.empmodelout$parammleset["AR1",][c("phi","sigsq")], betas.rrphylo=comb_betas.rrphylo,tree=tree_with_bodymass)
## $z
## phi
## 4.498322
##
## $pval
## pval
## 6.849175e-06
##############
####AR1 test##
##############
#########Combined tree##############
comb_betas.rrphylo<-compute_betas_rrphylo(tree_with_bodymass, y)$betas.rrphylo
combAr1test <- test.phi(param=comb.empmodelout$parammleset["AR1",][c("phi","sigsq")], betas.rrphylo=comb_betas.rrphylo,tree=tree_with_bodymass)
#########diurnal tree##############
diurnal_betas.rrphylo<-compute_betas_rrphylo(diurnal_subtree_with_bodymass, y)$betas.rrphylo
diurnalAr1test <- test.phi(param=diurnal.empmodelout$parammleset["AR1",][c("phi","sigsq")], betas.rrphylo=diurnal_betas.rrphylo,tree=diurnal_subtree_with_bodymass)
#########nocturnal tree##############
nocturnal_betas.rrphylo<-compute_betas_rrphylo(nocturnal_subtree_with_bodymass, y)$betas.rrphylo
nocturnalAr1test <- test.phi(param=nocturnal.empmodelout$parammleset["AR1",][c("phi","sigsq")], betas.rrphylo=nocturnal_betas.rrphylo,tree=nocturnal_subtree_with_bodymass)
treecase<-c("combined tree", "diurnal tree", "nocturnal tree")
zval<-c(combAr1test$z,diurnalAr1test$z,nocturnalAr1test$z)
pval<-c(combAr1test$pval,diurnalAr1test$pval,nocturnalAr1test$pval)
sig<-ifelse(pval<0.05,"Y","N")
ar1testtb<-data.frame(zval=zval,pval=pval,sig=sig)
rownames(ar1testtb)<-treecase
ar1testtb
## zval pval sig
## combined tree 4.4983225 6.849175e-06 Y
## diurnal tree 4.3033649 1.682235e-05 Y
## nocturnal tree -0.9013365 3.674094e-01 N
xtable(ar1testtb,digits=3)
## % latex table generated in R 4.4.0 by xtable 1.8-4 package
## % Sat May 4 14:12:54 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrl}
## \hline
## & zval & pval & sig \\
## \hline
## combined tree & 4.498 & 0.000 & Y \\
## diurnal tree & 4.303 & 0.000 & Y \\
## nocturnal tree & -0.901 & 0.367 & N \\
## \hline
## \end{tabular}
## \end{table}
#################################
#### LRT across models ##
#################################
LL_diurnal <- diurnal.empmodelout$fitset
LL_nocturnal <- nocturnal.empmodelout$fitset
LL_comb <- comb.empmodelout$fitset
# Compute the test statistic
test_statistic <- 2 * ((LL_diurnal$loglikset + LL_nocturnal$loglikset)-LL_comb$loglikset)
test_statistic
## [1] 508.34 495.84 498.48 498.38 495.60
#names(test_statistic)<-c(paste(LL_comb$modelset,".chi2",sep=""))
test_statistic
## [1] 508.34 495.84 498.48 498.38 495.60
# Perform the chi-square test
p_value.array <- 1 - pchisq(test_statistic, df=LL_comb$paranumset) # two AR(1) - one AR(1)
#names(p_value.array)<-c(paste(LL_comb$modelset,".pal",sep=""))
p_value.array # zero seems make sense as the rate is from the y simulated from
## [1] 0 0 0 0 0
tworates<-ifelse(p_value.array<0.05,"Y","N")
#names(tworates)<-rep("two rates?",5)
tworates
## [1] "Y" "Y" "Y" "Y" "Y"
LL_diurnal$loglikset
## [1] -42.55 -42.13 -42.26 -42.26 -42.26
LL_nocturnal$loglikset
## [1] 87.33 84.92 85.97 85.92 84.53
LL_comb$loglikset
## [1] -209.39 -205.13 -205.53 -205.53 -205.53
mset<-LL_comb$modelset
diurnallikeset<-LL_diurnal$loglikset
nocturnallikeset<-LL_nocturnal$loglikset
comblikeset<-LL_comb$loglikset
df2rates1rate<-data.frame(
diurnallikeset=diurnallikeset,
nocturnallikeset=nocturnallikeset,
comblikeset=comblikeset,
test_statistic=test_statistic,
p_value.array=p_value.array,tworates=tworates
)
rownames(df2rates1rate)<-LL_diurnal$modelset
colnames(df2rates1rate)<-c("log L_di","log L_no","log L_dino","chi2","pval")
df2rates1rate
## log L_di log L_no log L_dino chi2 pval NA
## AR1 -42.55 87.33 -209.39 508.34 0 Y
## AR2 -42.13 84.92 -205.13 495.84 0 Y
## ARMA1 -42.26 85.97 -205.53 498.48 0 Y
## ARMA21 -42.26 85.92 -205.53 498.38 0 Y
## ARMA22 -42.26 84.53 -205.53 495.60 0 Y
xtable(df2rates1rate)
## % latex table generated in R 4.4.0 by xtable 1.8-4 package
## % Sat May 4 14:12:54 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrrl}
## \hline
## & log L\_di & log L\_no & log L\_dino & chi2 & pval & NA \\
## \hline
## AR1 & -42.55 & 87.33 & -209.39 & 508.34 & 0.00 & Y \\
## AR2 & -42.13 & 84.92 & -205.13 & 495.84 & 0.00 & Y \\
## ARMA1 & -42.26 & 85.97 & -205.53 & 498.48 & 0.00 & Y \\
## ARMA21 & -42.26 & 85.92 & -205.53 & 498.38 & 0.00 & Y \\
## ARMA22 & -42.26 & 84.53 & -205.53 & 495.60 & 0.00 & Y \\
## \hline
## \end{tabular}
## \end{table}
# Do it separately, just call from the previous diurnal case and nocturnal case.
# testrate2subtrees_output <- testrate2subtrees(tree=tree_with_bodymass, y=BodyMass)
# tworatedf<- treeARanova(testrate2subtrees_output)
# tworatedf
# xtable(tworatedf,digits=3)