rm(list=ls())
setwd("~/Dropbox/FCU/Teaching/Mentoring/2022Spring-Current_Professor/2022Fall/YourueiMin/rcode_dcj/empricaldataanalysis/matchedtreetraitTony/HomeRange/")
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_HomeRange<-read.csv("diurnal_HomeRange.csv")
diurnal_HomeRange
## Species HomeRange_ha High_range Low_range
## 1 Alouatta_macconnelli 25.810000 NA NA
## 2 Brachyteles_arachnoides 205.000000 NA NA
## 3 Cacajao_calvus 1746.500000 NA NA
## 4 Cebus_kaapori 62.400000 NA NA
## 5 Cebus_olivaceus 257.000000 NA NA
## 6 Cercocebus_agilis 251.500000 NA NA
## 7 Chiropotes_chiropotes 142.000000 NA NA
## 8 Colobus_angolensis 1405.500000 NA NA
## 9 Eulemur_mongoz 3.925000 2.90 2.80
## 10 Hapalemur_aureus 26.000000 NA NA
## 11 Hapalemur_occidentalis 26.000000 NA NA
## 12 Lemur_catta 15.600000 35.00 16.00
## 13 Macaca_sylvanus 537.500000 700.50 253.50
## 14 Mandrillus_sphinx 6160.500000 2850.00 500.00
## 15 Papio_kindae NA NA NA
## 16 Pithecia_hirsuta 24.900000 41.00 9.70
## 17 Propithecus_candidus 38.850000 40.20 37.50
## 18 Propithecus_coquereli 7.625000 8.50 6.75
## 19 Propithecus_diadema 46.305000 62.59 27.11
## 20 Propithecus_edwardsi 89.580000 46.13 23.26
## 21 Propithecus_perrieri 15.520000 1.07 1.01
## 22 Propithecus_tattersalli 9.425000 12.15 6.70
## 23 Propithecus_verreauxi 6.083333 8.50 1.00
## 24 Saimiri_ustus 50.000000 NA NA
## 25 Semnopithecus_entellus 78.100000 106.00 45.00
## 26 Varecia_rubra 17.525000 17.55 14.36
## 27 Varecia_variegata 127.835000 150.00 100.00
length(diurnal_HomeRange$Species)
## [1] 27
same<-intersect(tiplabel,diurnal_HomeRange$Species)
diurnal_HomeRange$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
# Get the names of the tips to drop
tips_to_drop <- setdiff(tiplabel, diurnal_HomeRange$Species)
tips_to_drop<-c(tips_to_drop, "Papio_kindae")# "Pithecia_hirsuta","Pithecia_vanzolinii","Propithecus_candidus")
# Drop the tips
diurnal_subtree_with_HomeRange <- drop.tip(diurnal_subtree, tips_to_drop)
plot(diurnal_subtree_with_HomeRange)
diurnal_HomeRange <- diurnal_HomeRange[!(diurnal_HomeRange$Species %in% c("Papio_kindae")), ]
length(diurnal_subtree_with_HomeRange$tip.label)
## [1] 26
dim(diurnal_HomeRange)
## [1] 26 4
y<-diurnal_HomeRange$HomeRange_ha
names(y)<-diurnal_HomeRange$Species
y
## Alouatta_macconnelli Brachyteles_arachnoides Cacajao_calvus
## 25.810000 205.000000 1746.500000
## Cebus_kaapori Cebus_olivaceus Cercocebus_agilis
## 62.400000 257.000000 251.500000
## Chiropotes_chiropotes Colobus_angolensis Eulemur_mongoz
## 142.000000 1405.500000 3.925000
## Hapalemur_aureus Hapalemur_occidentalis Lemur_catta
## 26.000000 26.000000 15.600000
## Macaca_sylvanus Mandrillus_sphinx Pithecia_hirsuta
## 537.500000 6160.500000 24.900000
## Propithecus_candidus Propithecus_coquereli Propithecus_diadema
## 38.850000 7.625000 46.305000
## Propithecus_edwardsi Propithecus_perrieri Propithecus_tattersalli
## 89.580000 15.520000 9.425000
## Propithecus_verreauxi Saimiri_ustus Semnopithecus_entellus
## 6.083333 50.000000 78.100000
## Varecia_rubra Varecia_variegata
## 17.525000 127.835000
#######################################
############ diurnal Analysis #######
#######################################
mfit<-modelfit(trait=y, tree=diurnal_subtree_with_HomeRange)
diurnal.empmodelout<-empmodel(mfit)
diurnal.empmodelout$parammleset
## phi phi2 theta theta2 sigsq
## AR1 0.36623078 NA NA NA 19385.97
## AR2 0.42133197 -0.13740745 NA NA 15699.47
## ARMA1 0.26802555 NA 0.116330168 NA 19385.24
## ARMA21 0.28115410 -0.01311716 0.116317460 NA 19385.22
## ARMA22 -0.04043964 0.30840495 0.035641588 0.08079197 12419.87
## Weighted Estimate 0.42291020 -0.13465930 0.001163302 0.00000000 15967.06
diurnal.empmodelout$fitset
## modelset paranumset loglikset aicset aicweightset rankset
## 1 AR1 2 -311.40 626.80 0.02 2
## 2 AR2 3 -306.23 618.46 0.98 1
## 3 ARMA1 3 -311.40 628.80 0.01 3
## 4 ARMA21 4 -311.40 630.80 0.00 4
## 5 ARMA22 5 -314.23 638.46 0.00 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
## % Sun May 5 09:00:43 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrr}
## \hline
## & phi & phi2 & theta & theta2 & sigsq \\
## \hline
## AR1 & 0.366 & & & & 19385.975 \\
## AR2 & 0.421 & -0.137 & & & 15699.473 \\
## ARMA1 & 0.268 & & 0.116 & & 19385.244 \\
## ARMA21 & 0.281 & -0.013 & 0.116 & & 19385.215 \\
## ARMA22 & -0.040 & 0.308 & 0.036 & 0.081 & 12419.868 \\
## Weighted Estimate & 0.423 & -0.135 & 0.001 & 0.000 & 15967.056 \\
## \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
## % Sun May 5 09:00:43 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rlrrrrr}
## \hline
## & modelset & paranumset & loglikset & aicset & aicweightset & rankset \\
## \hline
## 1 & AR1 & 2.000 & -311.400 & 626.800 & 0.020 & 2.000 \\
## 2 & AR2 & 3.000 & -306.230 & 618.460 & 0.980 & 1.000 \\
## 3 & ARMA1 & 3.000 & -311.400 & 628.800 & 0.010 & 3.000 \\
## 4 & ARMA21 & 4.000 & -311.400 & 630.800 & 0.000 & 4.000 \\
## 5 & ARMA22 & 5.000 & -314.230 & 638.460 & 0.000 & 5.000 \\
## \hline
## \end{tabular}
## \end{table}
## add up test.phi here
diurnal_betas.rrphylo<-compute_betas_rrphylo(diurnal_subtree_with_HomeRange, y)$betas.rrphylo
test.phi(param=diurnal.empmodelout$parammleset["AR1",][c("phi","sigsq")], betas.rrphylo=diurnal_betas.rrphylo,tree=diurnal_subtree_with_HomeRange)
## $z
## phi
## 3.58281
##
## $pval
## pval
## 0.0003399173
###########################
##### 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_HomeRange<-read.csv("nocturnal_HomeRange.csv")
length(nocturnal_HomeRange$Species)
## [1] 42
same<-intersect(tiplabel,nocturnal_HomeRange$Species)
nocturnal_HomeRange$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
# Get the names of the tips to drop
tips_to_drop <- setdiff(tiplabel, nocturnal_HomeRange$Species)
tips_to_drop
## [1] "Nycticebus_javanicus" "Nycticebus_menagensis"
## [3] "Galagoides_demidoff" "Galago_granti"
## [5] "Galago_alleni" "Galago_gabonensis"
## [7] "Otolemur_monteiri" "Lepilemur_manasamody"
## [9] "Lepilemur_tymerlachsoni" "Lepilemur_septentrionalis"
## [11] "Lepilemur_hubbardorum" "Lepilemur_randrianasoloi"
# Drop the tips
nocturnal_subtree_with_HomeRange <- drop.tip(nocturnal_subtree, tips_to_drop)
plot(nocturnal_subtree_with_HomeRange)
nocturnal_HomeRange
## Species HomeRange_ha High_range Low_range
## 1 Arctocebus_aureus NA NA NA
## 2 Arctocebus_calabarensis 10.000000 NA NA
## 3 Cheirogaleus_major 3.000000 NA NA
## 4 Cheirogaleus_medius 1.560000 NA NA
## 5 Euoticus_elegantulus 6.700000 NA NA
## 6 Galago_gallarum 2.000000 3.00 1.00
## 7 Galago_matschiei NA NA NA
## 8 Galago_moholi 6.050000 14.00 4.60
## 9 Galago_senegalensis 8.166667 NA NA
## 10 Galagoides_cocos 1.950000 2.60 1.30
## 11 Galagoides_orinus NA NA NA
## 12 Galagoides_rondoensis NA NA NA
## 13 Galagoides_thomasi 1.600000 NA NA
## 14 Galagoides_zanzibaricus 2.150000 3.10 0.90
## 15 Lepilemur_aeeclis NA NA NA
## 16 Lepilemur_ahmansonorum NA NA NA
## 17 Lepilemur_ankaranensis 3.000000 NA NA
## 18 Lepilemur_betsileo NA NA NA
## 19 Lepilemur_dorsalis NA NA NA
## 20 Lepilemur_edwardsi 1.090000 1.70 0.67
## 21 Lepilemur_fleuretae 5.160000 7.28 3.04
## 22 Lepilemur_grewcockorum NA NA NA
## 23 Lepilemur_jamesorum NA NA NA
## 24 Lepilemur_leucopus 0.500000 NA NA
## 25 Lepilemur_microdon 0.500000 NA NA
## 26 Lepilemur_milanoii NA NA NA
## 27 Lepilemur_mittermeieri 3.600000 3.60 1.40
## 28 Lepilemur_mustelinus 0.900000 NA NA
## 29 Lepilemur_otto NA NA NA
## 30 Lepilemur_petteri 0.300000 NA NA
## 31 Lepilemur_ruficaudatus 0.800000 NA NA
## 32 Lepilemur_sahamalazensis 0.840000 4.04 0.18
## 33 Lepilemur_seali NA NA NA
## 34 Lepilemur_wrightae NA NA NA
## 35 Loris_lydekkerianus 14.900000 20.55 7.10
## 36 Microcebus_myoxinus NA NA NA
## 37 Nycticebus_bengalensis 15.600000 20.68 12.94
## 38 Nycticebus_coucang 6.300000 25.00 0.40
## 39 Nycticebus_pygmaeus 9.300000 12.55 6.05
## 40 Otolemur_crassicaudatus 8.500000 10.00 7.00
## 41 Otolemur_garnettii 14.170000 17.80 10.80
## 42 Perodicticus_potto 7.500000 NA NA
y<-nocturnal_HomeRange$HomeRange_ha
names(y)<-nocturnal_HomeRange$Species
length(y)
## [1] 42
length(nocturnal_subtree_with_HomeRange$tip.label)
## [1] 42
tip.to.drop<-names(y)[is.na(y)]
y
## Arctocebus_aureus Arctocebus_calabarensis Cheirogaleus_major
## NA 10.000000 3.000000
## Cheirogaleus_medius Euoticus_elegantulus Galago_gallarum
## 1.560000 6.700000 2.000000
## Galago_matschiei Galago_moholi Galago_senegalensis
## NA 6.050000 8.166667
## Galagoides_cocos Galagoides_orinus Galagoides_rondoensis
## 1.950000 NA NA
## Galagoides_thomasi Galagoides_zanzibaricus Lepilemur_aeeclis
## 1.600000 2.150000 NA
## Lepilemur_ahmansonorum Lepilemur_ankaranensis Lepilemur_betsileo
## NA 3.000000 NA
## Lepilemur_dorsalis Lepilemur_edwardsi Lepilemur_fleuretae
## NA 1.090000 5.160000
## Lepilemur_grewcockorum Lepilemur_jamesorum Lepilemur_leucopus
## NA NA 0.500000
## Lepilemur_microdon Lepilemur_milanoii Lepilemur_mittermeieri
## 0.500000 NA 3.600000
## Lepilemur_mustelinus Lepilemur_otto Lepilemur_petteri
## 0.900000 NA 0.300000
## Lepilemur_ruficaudatus Lepilemur_sahamalazensis Lepilemur_seali
## 0.800000 0.840000 NA
## Lepilemur_wrightae Loris_lydekkerianus Microcebus_myoxinus
## NA 14.900000 NA
## Nycticebus_bengalensis Nycticebus_coucang Nycticebus_pygmaeus
## 15.600000 6.300000 9.300000
## Otolemur_crassicaudatus Otolemur_garnettii Perodicticus_potto
## 8.500000 14.170000 7.500000
y<-y[names(y)[!is.na(y)]]
# Drop the tips
tips_to_drop <- setdiff(nocturnal_subtree_with_HomeRange$tip.label,names(y) )
tips_to_drop
## [1] "Arctocebus_aureus" "Galagoides_orinus" "Galagoides_rondoensis"
## [4] "Galago_matschiei" "Microcebus_myoxinus" "Lepilemur_otto"
## [7] "Lepilemur_grewcockorum" "Lepilemur_milanoii" "Lepilemur_ahmansonorum"
## [10] "Lepilemur_dorsalis" "Lepilemur_aeeclis" "Lepilemur_betsileo"
## [13] "Lepilemur_jamesorum" "Lepilemur_wrightae" "Lepilemur_seali"
sort(names(y))
## [1] "Arctocebus_calabarensis" "Cheirogaleus_major"
## [3] "Cheirogaleus_medius" "Euoticus_elegantulus"
## [5] "Galago_gallarum" "Galago_moholi"
## [7] "Galago_senegalensis" "Galagoides_cocos"
## [9] "Galagoides_thomasi" "Galagoides_zanzibaricus"
## [11] "Lepilemur_ankaranensis" "Lepilemur_edwardsi"
## [13] "Lepilemur_fleuretae" "Lepilemur_leucopus"
## [15] "Lepilemur_microdon" "Lepilemur_mittermeieri"
## [17] "Lepilemur_mustelinus" "Lepilemur_petteri"
## [19] "Lepilemur_ruficaudatus" "Lepilemur_sahamalazensis"
## [21] "Loris_lydekkerianus" "Nycticebus_bengalensis"
## [23] "Nycticebus_coucang" "Nycticebus_pygmaeus"
## [25] "Otolemur_crassicaudatus" "Otolemur_garnettii"
## [27] "Perodicticus_potto"
sort(nocturnal_subtree_with_HomeRange$tip.label)
## [1] "Arctocebus_aureus" "Arctocebus_calabarensis"
## [3] "Cheirogaleus_major" "Cheirogaleus_medius"
## [5] "Euoticus_elegantulus" "Galago_gallarum"
## [7] "Galago_matschiei" "Galago_moholi"
## [9] "Galago_senegalensis" "Galagoides_cocos"
## [11] "Galagoides_orinus" "Galagoides_rondoensis"
## [13] "Galagoides_thomasi" "Galagoides_zanzibaricus"
## [15] "Lepilemur_aeeclis" "Lepilemur_ahmansonorum"
## [17] "Lepilemur_ankaranensis" "Lepilemur_betsileo"
## [19] "Lepilemur_dorsalis" "Lepilemur_edwardsi"
## [21] "Lepilemur_fleuretae" "Lepilemur_grewcockorum"
## [23] "Lepilemur_jamesorum" "Lepilemur_leucopus"
## [25] "Lepilemur_microdon" "Lepilemur_milanoii"
## [27] "Lepilemur_mittermeieri" "Lepilemur_mustelinus"
## [29] "Lepilemur_otto" "Lepilemur_petteri"
## [31] "Lepilemur_ruficaudatus" "Lepilemur_sahamalazensis"
## [33] "Lepilemur_seali" "Lepilemur_wrightae"
## [35] "Loris_lydekkerianus" "Microcebus_myoxinus"
## [37] "Nycticebus_bengalensis" "Nycticebus_coucang"
## [39] "Nycticebus_pygmaeus" "Otolemur_crassicaudatus"
## [41] "Otolemur_garnettii" "Perodicticus_potto"
tips_to_drop
## [1] "Arctocebus_aureus" "Galagoides_orinus" "Galagoides_rondoensis"
## [4] "Galago_matschiei" "Microcebus_myoxinus" "Lepilemur_otto"
## [7] "Lepilemur_grewcockorum" "Lepilemur_milanoii" "Lepilemur_ahmansonorum"
## [10] "Lepilemur_dorsalis" "Lepilemur_aeeclis" "Lepilemur_betsileo"
## [13] "Lepilemur_jamesorum" "Lepilemur_wrightae" "Lepilemur_seali"
nocturnal_subtree_with_HomeRange <- drop.tip(nocturnal_subtree_with_HomeRange, tips_to_drop)
y
## Arctocebus_calabarensis Cheirogaleus_major Cheirogaleus_medius
## 10.000000 3.000000 1.560000
## Euoticus_elegantulus Galago_gallarum Galago_moholi
## 6.700000 2.000000 6.050000
## Galago_senegalensis Galagoides_cocos Galagoides_thomasi
## 8.166667 1.950000 1.600000
## Galagoides_zanzibaricus Lepilemur_ankaranensis Lepilemur_edwardsi
## 2.150000 3.000000 1.090000
## Lepilemur_fleuretae Lepilemur_leucopus Lepilemur_microdon
## 5.160000 0.500000 0.500000
## Lepilemur_mittermeieri Lepilemur_mustelinus Lepilemur_petteri
## 3.600000 0.900000 0.300000
## Lepilemur_ruficaudatus Lepilemur_sahamalazensis Loris_lydekkerianus
## 0.800000 0.840000 14.900000
## Nycticebus_bengalensis Nycticebus_coucang Nycticebus_pygmaeus
## 15.600000 6.300000 9.300000
## Otolemur_crassicaudatus Otolemur_garnettii Perodicticus_potto
## 8.500000 14.170000 7.500000
length(y)
## [1] 27
length(nocturnal_subtree_with_HomeRange$tip.label)
## [1] 27
#######################################
############ nocturnal Analysis #######
#######################################
#compute_betas_rrphylo(nocturnal_subtree_with_HomeRange, y)$betas.rrphylo
mfit<-modelfit(trait=y, tree=nocturnal_subtree_with_HomeRange)
nocturnal.empmodelout<-empmodel(mfit)
nocturnal.empmodelout$parammleset
## phi phi2 theta theta2 sigsq
## AR1 0.23354275 NA NA NA 0.9599596
## AR2 0.27737680 -0.20538307 NA NA 0.9788471
## ARMA1 -0.05165768 NA 0.31376346 NA 0.9599020
## ARMA21 -0.06462954 0.01293143 0.31381363 NA 0.9599012
## ARMA22 0.19614925 -0.24780564 -0.87172312 1.18547835 0.9599022
## Weighted Estimate 0.19132480 -0.08030128 0.04532074 0.02370957 0.9669353
nocturnal.empmodelout$fitset
## modelset paranumset loglikset aicset aicweightset rankset
## 1 AR1 2 -71.32 146.65 0.41 1
## 2 AR2 3 -70.41 146.83 0.37 2
## 3 ARMA1 3 -71.32 148.64 0.15 3
## 4 ARMA21 4 -71.32 150.64 0.05 4
## 5 ARMA22 5 -71.32 152.64 0.02 5
xtable(nocturnal.empmodelout$parammleset,digits=3)
## % latex table generated in R 4.4.0 by xtable 1.8-4 package
## % Sun May 5 09:00:44 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrr}
## \hline
## & phi & phi2 & theta & theta2 & sigsq \\
## \hline
## AR1 & 0.234 & & & & 0.960 \\
## AR2 & 0.277 & -0.205 & & & 0.979 \\
## ARMA1 & -0.052 & & 0.314 & & 0.960 \\
## ARMA21 & -0.065 & 0.013 & 0.314 & & 0.960 \\
## ARMA22 & 0.196 & -0.248 & -0.872 & 1.185 & 0.960 \\
## Weighted Estimate & 0.191 & -0.080 & 0.045 & 0.024 & 0.967 \\
## \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
## % Sun May 5 09:00:44 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rlrrrrr}
## \hline
## & modelset & paranumset & loglikset & aicset & aicweightset & rankset \\
## \hline
## 1 & AR1 & 2.000 & -71.320 & 146.650 & 0.410 & 1.000 \\
## 2 & AR2 & 3.000 & -70.410 & 146.830 & 0.370 & 2.000 \\
## 3 & ARMA1 & 3.000 & -71.320 & 148.640 & 0.150 & 3.000 \\
## 4 & ARMA21 & 4.000 & -71.320 & 150.640 & 0.050 & 4.000 \\
## 5 & ARMA22 & 5.000 & -71.320 & 152.640 & 0.020 & 5.000 \\
## \hline
## \end{tabular}
## \end{table}
## add up test.phi here
nocturnal_betas.rrphylo<-compute_betas_rrphylo(nocturnal_subtree_with_HomeRange, y)$betas.rrphylo
test.phi(param=nocturnal.empmodelout$parammleset["AR1",][c("phi","sigsq")], betas.rrphylo=nocturnal_betas.rrphylo,tree=nocturnal_subtree_with_HomeRange)
## $z
## phi
## 2.347543
##
## $pval
## pval
## 0.0188977
####################################
#########Combined tree##############
####################################
tree_with_HomeRange<-bind.tree(diurnal_subtree_with_HomeRange, nocturnal_subtree_with_HomeRange, where = "root", position = 0, interactive = FALSE)
tree_with_HomeRange<-chronopl(tree_with_HomeRange, 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_HomeRange, 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_HomeRange, show.tip.label = FALSE)
HomeRange<-rbind(diurnal_HomeRange,nocturnal_HomeRange)
y<-HomeRange$HomeRange_ha
names(y)<-HomeRange$Species
y
## Alouatta_macconnelli Brachyteles_arachnoides Cacajao_calvus
## 25.810000 205.000000 1746.500000
## Cebus_kaapori Cebus_olivaceus Cercocebus_agilis
## 62.400000 257.000000 251.500000
## Chiropotes_chiropotes Colobus_angolensis Eulemur_mongoz
## 142.000000 1405.500000 3.925000
## Hapalemur_aureus Hapalemur_occidentalis Lemur_catta
## 26.000000 26.000000 15.600000
## Macaca_sylvanus Mandrillus_sphinx Pithecia_hirsuta
## 537.500000 6160.500000 24.900000
## Propithecus_candidus Propithecus_coquereli Propithecus_diadema
## 38.850000 7.625000 46.305000
## Propithecus_edwardsi Propithecus_perrieri Propithecus_tattersalli
## 89.580000 15.520000 9.425000
## Propithecus_verreauxi Saimiri_ustus Semnopithecus_entellus
## 6.083333 50.000000 78.100000
## Varecia_rubra Varecia_variegata Arctocebus_aureus
## 17.525000 127.835000 NA
## Arctocebus_calabarensis Cheirogaleus_major Cheirogaleus_medius
## 10.000000 3.000000 1.560000
## Euoticus_elegantulus Galago_gallarum Galago_matschiei
## 6.700000 2.000000 NA
## Galago_moholi Galago_senegalensis Galagoides_cocos
## 6.050000 8.166667 1.950000
## Galagoides_orinus Galagoides_rondoensis Galagoides_thomasi
## NA NA 1.600000
## Galagoides_zanzibaricus Lepilemur_aeeclis Lepilemur_ahmansonorum
## 2.150000 NA NA
## Lepilemur_ankaranensis Lepilemur_betsileo Lepilemur_dorsalis
## 3.000000 NA NA
## Lepilemur_edwardsi Lepilemur_fleuretae Lepilemur_grewcockorum
## 1.090000 5.160000 NA
## Lepilemur_jamesorum Lepilemur_leucopus Lepilemur_microdon
## NA 0.500000 0.500000
## Lepilemur_milanoii Lepilemur_mittermeieri Lepilemur_mustelinus
## NA 3.600000 0.900000
## Lepilemur_otto Lepilemur_petteri Lepilemur_ruficaudatus
## NA 0.300000 0.800000
## Lepilemur_sahamalazensis Lepilemur_seali Lepilemur_wrightae
## 0.840000 NA NA
## Loris_lydekkerianus Microcebus_myoxinus Nycticebus_bengalensis
## 14.900000 NA 15.600000
## Nycticebus_coucang Nycticebus_pygmaeus Otolemur_crassicaudatus
## 6.300000 9.300000 8.500000
## Otolemur_garnettii Perodicticus_potto
## 14.170000 7.500000
length(y)
## [1] 68
names(y)%in% tree_with_HomeRange$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 FALSE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE
## [37] FALSE FALSE TRUE TRUE FALSE FALSE TRUE FALSE FALSE TRUE TRUE FALSE
## [49] FALSE TRUE TRUE FALSE TRUE TRUE FALSE TRUE TRUE TRUE FALSE FALSE
## [61] TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE
tree_with_HomeRange$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
#######################################
############ Full Analysis ##########
#######################################
#compute_betas_rrphylo(tree_with_HomeRange, y)$betas.rrphylo
mfit<-modelfit(trait=y, tree=tree_with_HomeRange)
comb.empmodelout<-empmodel(mfit)
comb.empmodelout$parammleset
## phi phi2 theta theta2 sigsq
## AR1 0.3471208 NA NA NA 1270.459
## AR2 0.2868200 0.1416716 NA NA 1254.376
## ARMA1 1.0878505 NA -0.88061273 NA 1271.323
## ARMA21 0.2703806 0.8174643 -0.88060586 NA 1240.228
## ARMA22 0.3432883 0.7445734 -0.28408022 -0.5965398 1236.744
## Weighted Estimate 0.3084214 0.1413460 -0.03522444 0.0000000 1242.360
comb.empmodelout$fitset
## modelset paranumset loglikset aicset aicweightset rankset
## 1 AR1 2 -494.26 992.52 0.01 4
## 2 AR2 3 -488.64 983.29 0.94 1
## 3 ARMA1 3 -491.98 989.97 0.03 2
## 4 ARMA21 4 -491.94 991.88 0.01 3
## 5 ARMA22 5 -491.94 993.87 0.00 5
xtable(comb.empmodelout$parammleset,digits=3)
## % latex table generated in R 4.4.0 by xtable 1.8-4 package
## % Sun May 5 09:00:46 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrr}
## \hline
## & phi & phi2 & theta & theta2 & sigsq \\
## \hline
## AR1 & 0.347 & & & & 1270.459 \\
## AR2 & 0.287 & 0.142 & & & 1254.376 \\
## ARMA1 & 1.088 & & -0.881 & & 1271.323 \\
## ARMA21 & 0.270 & 0.817 & -0.881 & & 1240.228 \\
## ARMA22 & 0.343 & 0.745 & -0.284 & -0.597 & 1236.744 \\
## Weighted Estimate & 0.308 & 0.141 & -0.035 & 0.000 & 1242.360 \\
## \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
## % Sun May 5 09:00:46 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rlrrrrr}
## \hline
## & modelset & paranumset & loglikset & aicset & aicweightset & rankset \\
## \hline
## 1 & AR1 & 2.000 & -494.260 & 992.520 & 0.010 & 4.000 \\
## 2 & AR2 & 3.000 & -488.640 & 983.290 & 0.940 & 1.000 \\
## 3 & ARMA1 & 3.000 & -491.980 & 989.970 & 0.030 & 2.000 \\
## 4 & ARMA21 & 4.000 & -491.940 & 991.880 & 0.010 & 3.000 \\
## 5 & ARMA22 & 5.000 & -491.940 & 993.870 & 0.000 & 5.000 \\
## \hline
## \end{tabular}
## \end{table}
## add up test.phi here
comb_betas.rrphylo<-compute_betas_rrphylo(tree_with_HomeRange, y)$betas.rrphylo
test.phi(param=comb.empmodelout$parammleset["AR1",][c("phi","sigsq")], betas.rrphylo=comb_betas.rrphylo,tree=tree_with_HomeRange)
## $z
## phi
## 2.741864
##
## $pval
## pval
## 0.006109161
##############
####AR1 test##
##############
#########Combined tree##############
comb_betas.rrphylo<-compute_betas_rrphylo(tree_with_HomeRange, y)$betas.rrphylo
combAr1test <- test.phi(param=comb.empmodelout$parammleset["AR1",][c("phi","sigsq")], betas.rrphylo=comb_betas.rrphylo,tree=tree_with_HomeRange)
#########diurnal tree##############
diurnal_betas.rrphylo<-compute_betas_rrphylo(diurnal_subtree_with_HomeRange, y)$betas.rrphylo
diurnalAr1test <- test.phi(param=diurnal.empmodelout$parammleset["AR1",][c("phi","sigsq")], betas.rrphylo=diurnal_betas.rrphylo,tree=diurnal_subtree_with_HomeRange)
#########nocturnal tree##############
nocturnal_betas.rrphylo<-compute_betas_rrphylo(nocturnal_subtree_with_HomeRange, y)$betas.rrphylo
nocturnalAr1test <- test.phi(param=nocturnal.empmodelout$parammleset["AR1",][c("phi","sigsq")], betas.rrphylo=nocturnal_betas.rrphylo,tree=nocturnal_subtree_with_HomeRange)
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 2.741864 0.0061091606 Y
## diurnal tree 3.582810 0.0003399173 Y
## nocturnal tree 2.347543 0.0188977035 Y
xtable(ar1testtb,digits=3)
## % latex table generated in R 4.4.0 by xtable 1.8-4 package
## % Sun May 5 09:00:47 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrl}
## \hline
## & zval & pval & sig \\
## \hline
## combined tree & 2.742 & 0.006 & Y \\
## diurnal tree & 3.583 & 0.000 & Y \\
## nocturnal tree & 2.348 & 0.019 & 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] 223.08 224.00 218.52 218.44 212.78
#names(test_statistic)<-c(paste(LL_comb$modelset,".chi2",sep=""))
test_statistic
## [1] 223.08 224.00 218.52 218.44 212.78
# 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] -311.40 -306.23 -311.40 -311.40 -314.23
LL_nocturnal$loglikset
## [1] -71.32 -70.41 -71.32 -71.32 -71.32
LL_comb$loglikset
## [1] -494.26 -488.64 -491.98 -491.94 -491.94
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","diff.rate")
df2rates1rate
## log L_di log L_no log L_dino chi2 pval diff.rate
## AR1 -311.40 -71.32 -494.26 223.08 0 Y
## AR2 -306.23 -70.41 -488.64 224.00 0 Y
## ARMA1 -311.40 -71.32 -491.98 218.52 0 Y
## ARMA21 -311.40 -71.32 -491.94 218.44 0 Y
## ARMA22 -314.23 -71.32 -491.94 212.78 0 Y
xtable(df2rates1rate)
## % latex table generated in R 4.4.0 by xtable 1.8-4 package
## % Sun May 5 09:00:47 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrrl}
## \hline
## & log L\_di & log L\_no & log L\_dino & chi2 & pval & diff.rate \\
## \hline
## AR1 & -311.40 & -71.32 & -494.26 & 223.08 & 0.00 & Y \\
## AR2 & -306.23 & -70.41 & -488.64 & 224.00 & 0.00 & Y \\
## ARMA1 & -311.40 & -71.32 & -491.98 & 218.52 & 0.00 & Y \\
## ARMA21 & -311.40 & -71.32 & -491.94 & 218.44 & 0.00 & Y \\
## ARMA22 & -314.23 & -71.32 & -491.94 & 212.78 & 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_HomeRange, y=HomeRange)
# tworatedf<- treeARanova(testrate2subtrees_output)
# tworatedf
# xtable(tworatedf,digits=3)