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$BodyMass_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.422139858 NA NA NA 0.1285498
## AR2 0.525841697 -0.20565252 NA NA 0.1305855
## ARMA1 -0.138127147 NA 0.7297148 NA 0.1273114
## ARMA21 -0.004075822 -0.13405632 0.7297199 NA 0.1273115
## ARMA22 -0.008338243 -0.12979018 -0.1282900 0.8580067 0.1273112
## Weighted Estimate 0.259674763 -0.04957914 0.2223633 0.0257402 0.1271893
diurnal.empmodelout$fitset
## modelset paranumset loglikset aicset aicweightset rankset
## 1 AR1 2 -20.84 45.68 0.48 1
## 2 AR2 3 -20.85 47.71 0.17 3
## 3 ARMA1 3 -20.58 47.17 0.23 2
## 4 ARMA21 4 -20.58 49.17 0.08 4
## 5 ARMA22 5 -20.58 51.17 0.03 5
#empmodelout$parammleset[1,][c("phi","sigsq")]
xtable(diurnal.empmodelout$parammleset,digits=3)
## % latex table generated in R 4.3.2 by xtable 1.8-4 package
## % Wed May 1 09:21:31 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrr}
## \hline
## & phi & phi2 & theta & theta2 & sigsq \\
## \hline
## AR1 & 0.422 & & & & 0.129 \\
## AR2 & 0.526 & -0.206 & & & 0.131 \\
## ARMA1 & -0.138 & & 0.730 & & 0.127 \\
## ARMA21 & -0.004 & -0.134 & 0.730 & & 0.127 \\
## ARMA22 & -0.008 & -0.130 & -0.128 & 0.858 & 0.127 \\
## Weighted Estimate & 0.260 & -0.050 & 0.222 & 0.026 & 0.127 \\
## \hline
## \end{tabular}
## \end{table}
xtable(diurnal.empmodelout$fitset,digits=3)
## % latex table generated in R 4.3.2 by xtable 1.8-4 package
## % Wed May 1 09:21:31 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rlrrrrr}
## \hline
## & modelset & paranumset & loglikset & aicset & aicweightset & rankset \\
## \hline
## 1 & AR1 & 2.000 & -20.840 & 45.680 & 0.480 & 1.000 \\
## 2 & AR2 & 3.000 & -20.850 & 47.710 & 0.170 & 3.000 \\
## 3 & ARMA1 & 3.000 & -20.580 & 47.170 & 0.230 & 2.000 \\
## 4 & ARMA21 & 4.000 & -20.580 & 49.170 & 0.080 & 4.000 \\
## 5 & ARMA22 & 5.000 & -20.580 & 51.170 & 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.236677
##
## $pval
## pval
## 2.268521e-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
y<-nocturnal_BodyMass$BodyMass_kg
names(y)<-nocturnal_BodyMass$Species
length(y)
## [1] 47
#######################################
############ 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.2128315 NA NA NA 0.003452367
## AR2 -0.2301223 -0.080081194 NA NA 0.003563155
## ARMA1 0.5932837 NA -0.90119170 NA 0.003720689
## ARMA21 0.6473108 0.075365497 -0.99998166 NA 0.003694149
## ARMA22 -0.1069107 0.610399203 -0.24303692 -0.7568841 0.003866292
## Weighted Estimate -0.1805653 -0.001648781 -0.03703557 0.0000000 0.003466158
nocturnal.empmodelout$fitset
## modelset paranumset loglikset aicset aicweightset rankset
## 1 AR1 2 129.73 -255.45 0.93 1
## 2 AR2 3 127.35 -248.71 0.03 2
## 3 ARMA1 3 127.34 -248.69 0.03 3
## 4 ARMA21 4 127.29 -246.58 0.01 4
## 5 ARMA22 5 123.66 -237.32 0.00 5
xtable(nocturnal.empmodelout$parammleset,digits=3)
## % latex table generated in R 4.3.2 by xtable 1.8-4 package
## % Wed May 1 09:21:36 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrr}
## \hline
## & phi & phi2 & theta & theta2 & sigsq \\
## \hline
## AR1 & -0.213 & & & & 0.003 \\
## AR2 & -0.230 & -0.080 & & & 0.004 \\
## ARMA1 & 0.593 & & -0.901 & & 0.004 \\
## ARMA21 & 0.647 & 0.075 & -1.000 & & 0.004 \\
## ARMA22 & -0.107 & 0.610 & -0.243 & -0.757 & 0.004 \\
## Weighted Estimate & -0.181 & -0.002 & -0.037 & 0.000 & 0.003 \\
## \hline
## \end{tabular}
## \end{table}
xtable(nocturnal.empmodelout$fitset,digits=3)
## % latex table generated in R 4.3.2 by xtable 1.8-4 package
## % Wed May 1 09:21:36 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rlrrrrr}
## \hline
## & modelset & paranumset & loglikset & aicset & aicweightset & rankset \\
## \hline
## 1 & AR1 & 2.000 & 129.730 & -255.450 & 0.930 & 1.000 \\
## 2 & AR2 & 3.000 & 127.350 & -248.710 & 0.030 & 2.000 \\
## 3 & ARMA1 & 3.000 & 127.340 & -248.690 & 0.030 & 3.000 \\
## 4 & ARMA21 & 4.000 & 127.290 & -246.580 & 0.010 & 4.000 \\
## 5 & ARMA22 & 5.000 & 123.660 & -237.320 & 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
## -2.762642
##
## $pval
## pval
## 0.005733566
####################################
#########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$BodyMass_kg
names(y)<-BodyMass$Species
y
## Alouatta_macconnelli Ateles_marginatus Brachyteles_arachnoides
## 6.3700000 6.9600000 8.8400000
## Cacajao_calvus Cebus_kaapori Cebus_olivaceus
## 3.1700000 2.4250000 2.8800000
## Cercocebus_agilis Chiropotes_albinasus Chiropotes_chiropotes
## 7.5800000 2.8200000 2.7400000
## Colobus_angolensis Eulemur_mongoz Hapalemur_aureus
## 8.6000000 1.3750000 1.5150000
## Hapalemur_occidentalis Lemur_catta Macaca_sylvanus
## 0.9315000 2.2100000 13.5000000
## Mandrillus_sphinx Papio_kindae Propithecus_coquereli
## 23.6000000 13.0000000 3.7850000
## Propithecus_deckenii Propithecus_diadema Propithecus_edwardsi
## 2.7800000 6.1050000 5.7400000
## Propithecus_perrieri Propithecus_tattersalli Propithecus_verreauxi
## 4.4100000 3.5000000 2.9100000
## Saimiri_ustus Semnopithecus_entellus Varecia_rubra
## 0.8600000 11.4400000 3.2650000
## Varecia_variegata Arctocebus_aureus Arctocebus_calabarensis
## 3.5500000 0.2100000 0.3100000
## Cheirogaleus_major Cheirogaleus_medius Euoticus_elegantulus
## 0.5050000 0.1766667 0.2700000
## Galago_gallarum Galago_matschiei Galago_moholi
## 0.2000000 0.2100000 0.2000000
## Galago_senegalensis Galagoides_cocos Galagoides_orinus
## 0.2500000 0.1400000 0.0900000
## Galagoides_rondoensis Galagoides_thomasi Galagoides_zanzibaricus
## 0.6750000 0.1200000 0.1400000
## Lepilemur_aeeclis Lepilemur_ahmansonorum Lepilemur_ankaranensis
## 0.8600000 0.5950000 0.7600000
## Lepilemur_betsileo Lepilemur_dorsalis Lepilemur_edwardsi
## 1.1400000 0.7950000 0.9750000
## Lepilemur_fleuretae Lepilemur_grewcockorum Lepilemur_jamesorum
## 0.8900000 0.8200000 0.8250000
## Lepilemur_leucopus Lepilemur_microdon Lepilemur_milanoii
## 0.5850000 1.0350000 0.7150000
## Lepilemur_mittermeieri Lepilemur_mustelinus Lepilemur_otto
## 0.8600000 0.9800000 0.9390000
## Lepilemur_petteri Lepilemur_randrianasoloi Lepilemur_ruficaudatus
## 0.6200000 0.8550000 0.7800000
## Lepilemur_sahamalazensis Lepilemur_seali Lepilemur_septentrionalis
## 0.6900000 0.9500000 0.6275000
## Lepilemur_tymerlachsoni Lepilemur_wrightae Loris_lydekkerianus
## 0.8800000 1.4350000 0.2700000
## Microcebus_myoxinus Nycticebus_bengalensis Nycticebus_coucang
## 0.0400000 1.2600000 0.6500000
## Nycticebus_javanicus Nycticebus_menagensis Nycticebus_pygmaeus
## 0.7300000 0.2850000 0.4200000
## Otolemur_crassicaudatus Otolemur_garnettii Perodicticus_potto
## 1.1500000 0.7600000 0.8300000
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 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 TRUE TRUE
## [61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE 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 TRUE TRUE
## [61] TRUE TRUE 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.3715772 NA NA NA 1.170002
## AR2 0.3359825 0.06322020 NA NA 1.168561
## ARMA1 0.8747640 NA -0.6519170 NA 1.103444
## ARMA21 0.7839125 0.09085427 -0.6519269 NA 1.103433
## ARMA22 0.4684150 0.40634466 -0.2245233 -0.42738781 1.103447
## Weighted Estimate 0.7793255 0.05656506 -0.5655744 -0.03419102 1.108694
comb.empmodelout$fitset
## modelset paranumset loglikset aicset aicweightset rankset
## 1 AR1 2 -209.64 423.28 0.03 5
## 2 AR2 3 -208.06 422.12 0.05 4
## 3 ARMA1 3 -205.54 417.08 0.61 1
## 4 ARMA21 4 -205.54 419.08 0.23 2
## 5 ARMA22 5 -205.54 421.08 0.08 3
xtable(comb.empmodelout$parammleset,digits=3)
## % latex table generated in R 4.3.2 by xtable 1.8-4 package
## % Wed May 1 09:21:42 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrr}
## \hline
## & phi & phi2 & theta & theta2 & sigsq \\
## \hline
## AR1 & 0.372 & & & & 1.170 \\
## AR2 & 0.336 & 0.063 & & & 1.169 \\
## ARMA1 & 0.875 & & -0.652 & & 1.103 \\
## ARMA21 & 0.784 & 0.091 & -0.652 & & 1.103 \\
## ARMA22 & 0.468 & 0.406 & -0.225 & -0.427 & 1.103 \\
## Weighted Estimate & 0.779 & 0.057 & -0.566 & -0.034 & 1.109 \\
## \hline
## \end{tabular}
## \end{table}
xtable(comb.empmodelout$fitset,digits=3)
## % latex table generated in R 4.3.2 by xtable 1.8-4 package
## % Wed May 1 09:21:42 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rlrrrrr}
## \hline
## & modelset & paranumset & loglikset & aicset & aicweightset & rankset \\
## \hline
## 1 & AR1 & 2.000 & -209.640 & 423.280 & 0.030 & 5.000 \\
## 2 & AR2 & 3.000 & -208.060 & 422.120 & 0.050 & 4.000 \\
## 3 & ARMA1 & 3.000 & -205.540 & 417.080 & 0.610 & 1.000 \\
## 4 & ARMA21 & 4.000 & -205.540 & 419.080 & 0.230 & 2.000 \\
## 5 & ARMA22 & 5.000 & -205.540 & 421.080 & 0.080 & 3.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.228716
##
## $pval
## pval
## 2.350288e-05
##############
####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.228716 2.350288e-05 Y
## diurnal tree 4.236677 2.268521e-05 Y
## nocturnal tree -2.762642 5.733566e-03 Y
xtable(ar1testtb,digits=3)
## % latex table generated in R 4.3.2 by xtable 1.8-4 package
## % Wed May 1 09:21:48 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrl}
## \hline
## & zval & pval & sig \\
## \hline
## combined tree & 4.229 & 0.000 & Y \\
## diurnal tree & 4.237 & 0.000 & Y \\
## nocturnal tree & -2.763 & 0.006 & 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] 637.06 629.12 624.60 624.50 617.24
#names(test_statistic)<-c(paste(LL_comb$modelset,".chi2",sep=""))
test_statistic
## [1] 637.06 629.12 624.60 624.50 617.24
# 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] -20.84 -20.85 -20.58 -20.58 -20.58
LL_nocturnal$loglikset
## [1] 129.73 127.35 127.34 127.29 123.66
LL_comb$loglikset
## [1] -209.64 -208.06 -205.54 -205.54 -205.54
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 -20.84 129.73 -209.64 637.06 0 Y
## AR2 -20.85 127.35 -208.06 629.12 0 Y
## ARMA1 -20.58 127.34 -205.54 624.60 0 Y
## ARMA21 -20.58 127.29 -205.54 624.50 0 Y
## ARMA22 -20.58 123.66 -205.54 617.24 0 Y
xtable(df2rates1rate)
## % latex table generated in R 4.3.2 by xtable 1.8-4 package
## % Wed May 1 09:21:48 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrrl}
## \hline
## & log L\_di & log L\_no & log L\_dino & chi2 & pval & NA \\
## \hline
## AR1 & -20.84 & 129.73 & -209.64 & 637.06 & 0.00 & Y \\
## AR2 & -20.85 & 127.35 & -208.06 & 629.12 & 0.00 & Y \\
## ARMA1 & -20.58 & 127.34 & -205.54 & 624.60 & 0.00 & Y \\
## ARMA21 & -20.58 & 127.29 & -205.54 & 624.50 & 0.00 & Y \\
## ARMA22 & -20.58 & 123.66 & -205.54 & 617.24 & 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)