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$BodyMassFemale_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.4207331 NA NA NA 0.03117165
## AR2 0.4524003 -0.06377864 NA NA 0.03175290
## ARMA1 -0.1605595 NA 0.7353059 NA 0.03087015
## ARMA21 0.1340958 -0.29525485 0.7358054 NA 0.03087196
## ARMA22 -0.4133221 0.25222509 0.6718619 0.063899637 0.03087223
## Weighted Estimate 0.2232337 -0.02474626 0.2702048 0.001916989 0.03080092
diurnal.empmodelout$fitset
## modelset paranumset loglikset aicset aicweightset rankset
## 1 AR1 2 16.71 -29.42 0.53 1
## 2 AR2 3 15.92 -25.84 0.09 4
## 3 ARMA1 3 16.98 -27.96 0.25 2
## 4 ARMA21 4 16.98 -25.96 0.09 3
## 5 ARMA22 5 16.98 -23.96 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:19:44 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrr}
## \hline
## & phi & phi2 & theta & theta2 & sigsq \\
## \hline
## AR1 & 0.421 & & & & 0.031 \\
## AR2 & 0.452 & -0.064 & & & 0.032 \\
## ARMA1 & -0.161 & & 0.735 & & 0.031 \\
## ARMA21 & 0.134 & -0.295 & 0.736 & & 0.031 \\
## ARMA22 & -0.413 & 0.252 & 0.672 & 0.064 & 0.031 \\
## Weighted Estimate & 0.223 & -0.025 & 0.270 & 0.002 & 0.031 \\
## \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:19:44 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rlrrrrr}
## \hline
## & modelset & paranumset & loglikset & aicset & aicweightset & rankset \\
## \hline
## 1 & AR1 & 2.000 & 16.710 & -29.420 & 0.530 & 1.000 \\
## 2 & AR2 & 3.000 & 15.920 & -25.840 & 0.090 & 4.000 \\
## 3 & ARMA1 & 3.000 & 16.980 & -27.960 & 0.250 & 2.000 \\
## 4 & ARMA21 & 4.000 & 16.980 & -25.960 & 0.090 & 3.000 \\
## 5 & ARMA22 & 5.000 & 16.980 & -23.960 & 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
## 3.968701
##
## $pval
## pval
## 7.226562e-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$BodyMassFemale_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_microdon"
## [11] "Lepilemur_milanoii" "Lepilemur_mittermeieri"
## [13] "Lepilemur_mustelinus" "Lepilemur_otto"
## [15] "Lepilemur_randrianasoloi" "Lepilemur_seali"
## [17] "Lepilemur_septentrionalis" "Lepilemur_tymerlachsoni"
## [19] "Nycticebus_menagensis"
# Drop the tips
nocturnal_subtree_with_bodymass <- drop.tip(nocturnal_subtree_with_bodymass, tips_to_drop)
nocturnal_subtree_with_bodymass
##
## Phylogenetic tree with 28 tips and 27 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.3060 0.4430
## Cheirogaleus_medius Euoticus_elegantulus Galago_gallarum
## 0.1560 0.2610 0.2000
## Galago_matschiei Galago_moholi Galago_senegalensis
## 0.2120 0.1940 0.2250
## Galagoides_cocos Galagoides_orinus Galagoides_rondoensis
## 0.1300 NA 0.6600
## Galagoides_thomasi Galagoides_zanzibaricus Lepilemur_aeeclis
## 0.1300 0.1370 0.9090
## Lepilemur_ahmansonorum Lepilemur_ankaranensis Lepilemur_betsileo
## NA NA NA
## Lepilemur_dorsalis Lepilemur_edwardsi Lepilemur_fleuretae
## NA 0.9340 NA
## Lepilemur_grewcockorum Lepilemur_jamesorum Lepilemur_leucopus
## NA NA 0.5940
## Lepilemur_microdon Lepilemur_milanoii Lepilemur_mittermeieri
## NA NA NA
## Lepilemur_mustelinus Lepilemur_otto Lepilemur_petteri
## NA NA 0.5850
## Lepilemur_randrianasoloi Lepilemur_ruficaudatus Lepilemur_sahamalazensis
## NA 0.7185 0.6780
## Lepilemur_seali Lepilemur_septentrionalis Lepilemur_tymerlachsoni
## NA NA NA
## Lepilemur_wrightae Loris_lydekkerianus Microcebus_myoxinus
## 1.8500 0.2690 0.0307
## Nycticebus_bengalensis Nycticebus_coucang Nycticebus_javanicus
## 1.4000 0.6260 0.6260
## Nycticebus_menagensis Nycticebus_pygmaeus Otolemur_crassicaudatus
## NA 0.3760 1.1100
## Otolemur_garnettii Perodicticus_potto
## 0.7340 0.8360
y<-y[!is.na(y)]
y
## Arctocebus_calabarensis Cheirogaleus_major Cheirogaleus_medius
## 0.3060 0.4430 0.1560
## Euoticus_elegantulus Galago_gallarum Galago_matschiei
## 0.2610 0.2000 0.2120
## Galago_moholi Galago_senegalensis Galagoides_cocos
## 0.1940 0.2250 0.1300
## Galagoides_rondoensis Galagoides_thomasi Galagoides_zanzibaricus
## 0.6600 0.1300 0.1370
## Lepilemur_aeeclis Lepilemur_edwardsi Lepilemur_leucopus
## 0.9090 0.9340 0.5940
## Lepilemur_petteri Lepilemur_ruficaudatus Lepilemur_sahamalazensis
## 0.5850 0.7185 0.6780
## Lepilemur_wrightae Loris_lydekkerianus Microcebus_myoxinus
## 1.8500 0.2690 0.0307
## Nycticebus_bengalensis Nycticebus_coucang Nycticebus_javanicus
## 1.4000 0.6260 0.6260
## Nycticebus_pygmaeus Otolemur_crassicaudatus Otolemur_garnettii
## 0.3760 1.1100 0.7340
## Perodicticus_potto
## 0.8360
length(y)
## [1] 28
#######################################
############ 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.30078711 NA NA NA 0.006945290
## AR2 0.39367210 0.22040630 NA NA 0.007010080
## ARMA1 0.43457561 NA -0.99991393 NA 0.007331100
## ARMA21 0.47066231 -0.09348048 -0.99991863 NA 0.007474127
## ARMA22 0.05428971 0.37455739 0.07855326 -0.221988872 0.006931994
## Weighted Estimate 0.30751609 0.02854113 -0.06685194 -0.008879555 0.006979161
nocturnal.empmodelout$fitset
## modelset paranumset loglikset aicset aicweightset rankset
## 1 AR1 2 56.82 -109.65 0.82 1
## 2 AR2 3 55.34 -104.69 0.07 2
## 3 ARMA1 3 55.10 -104.19 0.05 3
## 4 ARMA21 4 54.84 -101.68 0.02 5
## 5 ARMA22 5 56.85 -103.69 0.04 4
xtable(nocturnal.empmodelout$parammleset,digits=3)
## % latex table generated in R 4.4.0 by xtable 1.8-4 package
## % Sat May 4 14:19:47 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrr}
## \hline
## & phi & phi2 & theta & theta2 & sigsq \\
## \hline
## AR1 & 0.301 & & & & 0.007 \\
## AR2 & 0.394 & 0.220 & & & 0.007 \\
## ARMA1 & 0.435 & & -1.000 & & 0.007 \\
## ARMA21 & 0.471 & -0.093 & -1.000 & & 0.007 \\
## ARMA22 & 0.054 & 0.375 & 0.079 & -0.222 & 0.007 \\
## Weighted Estimate & 0.308 & 0.029 & -0.067 & -0.009 & 0.007 \\
## \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:19:47 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rlrrrrr}
## \hline
## & modelset & paranumset & loglikset & aicset & aicweightset & rankset \\
## \hline
## 1 & AR1 & 2.000 & 56.820 & -109.650 & 0.820 & 1.000 \\
## 2 & AR2 & 3.000 & 55.340 & -104.690 & 0.070 & 2.000 \\
## 3 & ARMA1 & 3.000 & 55.100 & -104.190 & 0.050 & 3.000 \\
## 4 & ARMA21 & 4.000 & 54.840 & -101.680 & 0.020 & 5.000 \\
## 5 & ARMA22 & 5.000 & 56.850 & -103.690 & 0.040 & 4.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
## 3.053179
##
## $pval
## pval
## 0.00226431
####################################
#########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$BodyMassFemale_kg
names(y)<-BodyMass$Species
y
## Alouatta_macconnelli Ateles_marginatus Brachyteles_arachnoides
## 5.5500 5.8000 8.0700
## Cacajao_calvus Cebus_kaapori Cebus_olivaceus
## 2.8800 2.4000 2.5200
## Cercocebus_agilis Chiropotes_albinasus Chiropotes_chiropotes
## 5.6600 2.4900 2.5800
## Colobus_angolensis Eulemur_mongoz Hapalemur_aureus
## 7.4000 1.5600 1.3900
## Hapalemur_occidentalis Lemur_catta Macaca_sylvanus
## 1.1800 2.2100 11.0000
## Mandrillus_sphinx Papio_kindae Propithecus_coquereli
## 12.8000 9.9000 3.7560
## Propithecus_deckenii Propithecus_diadema Propithecus_edwardsi
## 2.6300 6.2100 5.8950
## Propithecus_perrieri Propithecus_tattersalli Propithecus_verreauxi
## 4.4400 3.5900 2.8700
## Saimiri_ustus Semnopithecus_entellus Varecia_rubra
## 0.7990 9.8900 3.4700
## Varecia_variegata Arctocebus_aureus Arctocebus_calabarensis
## 3.5200 NA 0.3060
## Cheirogaleus_major Cheirogaleus_medius Euoticus_elegantulus
## 0.4430 0.1560 0.2610
## Galago_gallarum Galago_matschiei Galago_moholi
## 0.2000 0.2120 0.1940
## Galago_senegalensis Galagoides_cocos Galagoides_orinus
## 0.2250 0.1300 NA
## Galagoides_rondoensis Galagoides_thomasi Galagoides_zanzibaricus
## 0.6600 0.1300 0.1370
## Lepilemur_aeeclis Lepilemur_ahmansonorum Lepilemur_ankaranensis
## 0.9090 NA NA
## Lepilemur_betsileo Lepilemur_dorsalis Lepilemur_edwardsi
## NA NA 0.9340
## Lepilemur_fleuretae Lepilemur_grewcockorum Lepilemur_jamesorum
## NA NA NA
## Lepilemur_leucopus Lepilemur_microdon Lepilemur_milanoii
## 0.5940 NA NA
## Lepilemur_mittermeieri Lepilemur_mustelinus Lepilemur_otto
## NA NA NA
## Lepilemur_petteri Lepilemur_randrianasoloi Lepilemur_ruficaudatus
## 0.5850 NA 0.7185
## Lepilemur_sahamalazensis Lepilemur_seali Lepilemur_septentrionalis
## 0.6780 NA NA
## Lepilemur_tymerlachsoni Lepilemur_wrightae Loris_lydekkerianus
## NA 1.8500 0.2690
## Microcebus_myoxinus Nycticebus_bengalensis Nycticebus_coucang
## 0.0307 1.4000 0.6260
## Nycticebus_javanicus Nycticebus_menagensis Nycticebus_pygmaeus
## 0.6260 NA 0.3760
## Otolemur_crassicaudatus Otolemur_garnettii Perodicticus_potto
## 1.1100 0.7340 0.8360
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 FALSE FALSE FALSE FALSE FALSE TRUE FALSE TRUE
## [61] TRUE FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE FALSE 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
#######################################
############ 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.37239902 NA NA NA 1.065007
## AR2 0.30392716 0.211847665 NA NA 1.037273
## ARMA1 0.74939706 NA -0.4197426 NA 1.034466
## ARMA21 0.75069148 -0.001296982 -0.4197400 NA 1.034466
## ARMA22 -0.01971449 0.769109963 0.1875373 -0.60727814 1.034465
## Weighted Estimate 0.40379798 0.150886886 -0.0927899 -0.01214556 1.028972
comb.empmodelout$fitset
## modelset paranumset loglikset aicset aicweightset rankset
## 1 AR1 2 -152.30 308.59 0.10 3
## 2 AR2 3 -149.47 304.94 0.64 1
## 3 ARMA1 3 -150.77 307.54 0.17 2
## 4 ARMA21 4 -150.77 309.54 0.06 4
## 5 ARMA22 5 -150.77 311.54 0.02 5
xtable(comb.empmodelout$parammleset,digits=3)
## % latex table generated in R 4.4.0 by xtable 1.8-4 package
## % Sat May 4 14:19:50 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrr}
## \hline
## & phi & phi2 & theta & theta2 & sigsq \\
## \hline
## AR1 & 0.372 & & & & 1.065 \\
## AR2 & 0.304 & 0.212 & & & 1.037 \\
## ARMA1 & 0.749 & & -0.420 & & 1.034 \\
## ARMA21 & 0.751 & -0.001 & -0.420 & & 1.034 \\
## ARMA22 & -0.020 & 0.769 & 0.188 & -0.607 & 1.034 \\
## Weighted Estimate & 0.404 & 0.151 & -0.093 & -0.012 & 1.029 \\
## \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:19:50 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rlrrrrr}
## \hline
## & modelset & paranumset & loglikset & aicset & aicweightset & rankset \\
## \hline
## 1 & AR1 & 2.000 & -152.300 & 308.590 & 0.100 & 3.000 \\
## 2 & AR2 & 3.000 & -149.470 & 304.940 & 0.640 & 1.000 \\
## 3 & ARMA1 & 3.000 & -150.770 & 307.540 & 0.170 & 2.000 \\
## 4 & ARMA21 & 4.000 & -150.770 & 309.540 & 0.060 & 4.000 \\
## 5 & ARMA22 & 5.000 & -150.770 & 311.540 & 0.020 & 5.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
## 3.395437
##
## $pval
## pval
## 0.0006851917
##############
####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 3.395437 6.851917e-04 Y
## diurnal tree 3.968701 7.226562e-05 Y
## nocturnal tree 3.053179 2.264310e-03 Y
xtable(ar1testtb,digits=3)
## % latex table generated in R 4.4.0 by xtable 1.8-4 package
## % Sat May 4 14:19:52 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrl}
## \hline
## & zval & pval & sig \\
## \hline
## combined tree & 3.395 & 0.001 & Y \\
## diurnal tree & 3.969 & 0.000 & Y \\
## nocturnal tree & 3.053 & 0.002 & Y \\
## \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] 451.66 441.46 445.70 445.18 449.20
#names(test_statistic)<-c(paste(LL_comb$modelset,".chi2",sep=""))
test_statistic
## [1] 451.66 441.46 445.70 445.18 449.20
# 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] 16.71 15.92 16.98 16.98 16.98
LL_nocturnal$loglikset
## [1] 56.82 55.34 55.10 54.84 56.85
LL_comb$loglikset
## [1] -152.30 -149.47 -150.77 -150.77 -150.77
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 16.71 56.82 -152.30 451.66 0 Y
## AR2 15.92 55.34 -149.47 441.46 0 Y
## ARMA1 16.98 55.10 -150.77 445.70 0 Y
## ARMA21 16.98 54.84 -150.77 445.18 0 Y
## ARMA22 16.98 56.85 -150.77 449.20 0 Y
xtable(df2rates1rate)
## % latex table generated in R 4.4.0 by xtable 1.8-4 package
## % Sat May 4 14:19:52 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrrl}
## \hline
## & log L\_di & log L\_no & log L\_dino & chi2 & pval & NA \\
## \hline
## AR1 & 16.71 & 56.82 & -152.30 & 451.66 & 0.00 & Y \\
## AR2 & 15.92 & 55.34 & -149.47 & 441.46 & 0.00 & Y \\
## ARMA1 & 16.98 & 55.10 & -150.77 & 445.70 & 0.00 & Y \\
## ARMA21 & 16.98 & 54.84 & -150.77 & 445.18 & 0.00 & Y \\
## ARMA22 & 16.98 & 56.85 & -150.77 & 449.20 & 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)