rm(list=ls())
setwd("/Users/tonyjhwueng/Dropbox/FCU/Teaching/Mentoring/2022Spring-Current_Professor/2022Fall/YourueiMin/rcode_dcj/empricaldataanalysis/henniges2023plant")
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
## Warning: package 'RRphylo' was built under R version 4.4.1
## Loading required package: emmeans
## Warning: package 'emmeans' was built under R version 4.4.1
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'
##
## 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
## Warning: package 'robustbase' was built under R version 4.4.1
## Loading required package: Matrix
#公的
set.seed(0930)
###########################
##### herba_subtree #####
###########################
herba_subtree<-read.tree("herba.nwk")
tiplabel<-herba_subtree$tip.label
tiplabel
## [1] "Viscaria_vulgaris" "Digitalis_purpurea"
## [3] "Mentha_longifolia" "Nepeta_cataria"
## [5] "Salvia_pratensis" "Verbascum_thapsus"
## [7] "Convolvulus_arvensis" "Symphytum_officinale"
## [9] "Galium_aparine" "Valeriana_officinalis"
## [11] "Angelica_archangelica" "Campanula_rotundifolia"
## [13] "Achillea_millefolium" "Tanacetum_vulgare"
## [15] "Aster_amellus" "Solidago_virgaurea"
## [17] "Eupatorium_cannabinum" "Echinacea_purpurea"
## [19] "Rudbeckia_hirta" "Geranium_robertianum"
## [21] "Chamaenerion_angustifolium" "Filipendula_ulmaria"
## [23] "Sanguisorba_officinalis" "Alchemilla_vulgaris"
## [25] "Fragaria_vesca" "Potentilla_indica"
## [27] "Geum_rivale" "Lupinus_perennis"
## [29] "Lathyrus_vernus" "Anemone_nemorosa"
## [31] "Pulsatilla_vulgaris" "Delphinium_elatum"
## [33] "Aconitum_napellus" "Aquilegia_vulgaris"
## [35] "Thalictrum_aquilegiifolium"
length(tiplabel)
## [1] 35
tiplabel <- gsub("_", " ", tiplabel)
herba_subtree$tip.label <- tiplabel
tiplabel
## [1] "Viscaria vulgaris" "Digitalis purpurea"
## [3] "Mentha longifolia" "Nepeta cataria"
## [5] "Salvia pratensis" "Verbascum thapsus"
## [7] "Convolvulus arvensis" "Symphytum officinale"
## [9] "Galium aparine" "Valeriana officinalis"
## [11] "Angelica archangelica" "Campanula rotundifolia"
## [13] "Achillea millefolium" "Tanacetum vulgare"
## [15] "Aster amellus" "Solidago virgaurea"
## [17] "Eupatorium cannabinum" "Echinacea purpurea"
## [19] "Rudbeckia hirta" "Geranium robertianum"
## [21] "Chamaenerion angustifolium" "Filipendula ulmaria"
## [23] "Sanguisorba officinalis" "Alchemilla vulgaris"
## [25] "Fragaria vesca" "Potentilla indica"
## [27] "Geum rivale" "Lupinus perennis"
## [29] "Lathyrus vernus" "Anemone nemorosa"
## [31] "Pulsatilla vulgaris" "Delphinium elatum"
## [33] "Aconitum napellus" "Aquilegia vulgaris"
## [35] "Thalictrum aquilegiifolium"
herba_Gsize<-read.csv("herba.csv")
#length(diurnal_BodyMass$Species)
#diurnal_BodyMass <- diurnal_BodyMass[!is.na(diurnal_BodyMass$BodyMassMale_kg), ]
same<-intersect(tiplabel,herba_Gsize$Species)
herba_Gsize <- herba_Gsize[herba_Gsize$Species %in% same, ]
# 更新樹中的物種列表
herba_subtree <- drop.tip(herba_subtree, setdiff(tiplabel, same))
tiplabel <- herba_subtree$tip.label
# 確保樹和體重數據中物種一致
same <- intersect(tiplabel, herba_Gsize$Species)
herba_Gsize <- herba_Gsize[herba_Gsize$Species %in% same, ]
#隨機刪除20個物種
#tips_to_drop <- sample(tiplabel, 20)
#diurnal_subtree <- drop.tip(diurnal_subtree, tips_to_drop)
#tiplabel <- diurnal_subtree$tip.label
# 更新體重數據中的物種列表
#same <- intersect(tiplabel, diurnal_BodyMass$Species)
#diurnal_BodyMass <- diurnal_BodyMass[diurnal_BodyMass$Species %in% same, ]
# 6. 繪製更新後的樹
plot(herba_subtree)
y<-herba_Gsize$GenomeSize
names(y)<-herba_Gsize$Species
# 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)
#y<-diurnal_BodyMass$BodyMassMale_kg
#names(y)<-diurnal_BodyMass$Species
#######################################
############ herba Analysis #######
#######################################
mfit<-modelfit(trait=y, tree=herba_subtree)
herba.empmodelout<-empmodel(mfit)
herba.empmodelout$parammleset
## phi phi2 theta theta2 sigsq
## AR1 0.4697693 NA NA NA 0.008960324
## AR2 0.4507656 -0.01943987 NA NA 0.009143424
## ARMA1 0.4913526 NA -0.02426855 NA 0.009019282
## ARMA21 0.1670061 0.16342028 0.14953385 NA 0.009030695
## ARMA22 -0.1275510 0.60854091 0.47593236 -0.47995182 0.009009934
## Weighted Estimate 0.4267291 0.03035786 0.02090160 -0.01439855 0.008899965
herba.empmodelout$fitset
## modelset paranumset loglikset aicset aicweightset rankset
## 1 AR1 2 57.27 -110.55 0.61 1
## 2 AR2 3 55.84 -105.69 0.05 4
## 3 ARMA1 3 57.27 -108.55 0.22 2
## 4 ARMA21 4 57.25 -106.50 0.08 3
## 5 ARMA22 5 57.27 -104.55 0.03 5
#empmodelout$parammleset[1,][c("phi","sigsq")]
xtable(herba.empmodelout$parammleset,digits=3)
## % latex table generated in R 4.4.0 by xtable 1.8-4 package
## % Sun Sep 29 21:42:05 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrr}
## \hline
## & phi & phi2 & theta & theta2 & sigsq \\
## \hline
## AR1 & 0.470 & & & & 0.009 \\
## AR2 & 0.451 & -0.019 & & & 0.009 \\
## ARMA1 & 0.491 & & -0.024 & & 0.009 \\
## ARMA21 & 0.167 & 0.163 & 0.150 & & 0.009 \\
## ARMA22 & -0.128 & 0.609 & 0.476 & -0.480 & 0.009 \\
## Weighted Estimate & 0.427 & 0.030 & 0.021 & -0.014 & 0.009 \\
## \hline
## \end{tabular}
## \end{table}
xtable(herba.empmodelout$fitset,digits=3)
## % latex table generated in R 4.4.0 by xtable 1.8-4 package
## % Sun Sep 29 21:42:05 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rlrrrrr}
## \hline
## & modelset & paranumset & loglikset & aicset & aicweightset & rankset \\
## \hline
## 1 & AR1 & 2.000 & 57.270 & -110.550 & 0.610 & 1.000 \\
## 2 & AR2 & 3.000 & 55.840 & -105.690 & 0.050 & 4.000 \\
## 3 & ARMA1 & 3.000 & 57.270 & -108.550 & 0.220 & 2.000 \\
## 4 & ARMA21 & 4.000 & 57.250 & -106.500 & 0.080 & 3.000 \\
## 5 & ARMA22 & 5.000 & 57.270 & -104.550 & 0.030 & 5.000 \\
## \hline
## \end{tabular}
## \end{table}
## add up test.phi here
herba_betas.rrphylo<-compute_betas_rrphylo(herba_subtree, y)$betas.rrphylo
test.phi(param=herba.empmodelout$parammleset["AR1",][c("phi","sigsq")], betas.rrphylo=herba_betas.rrphylo,tree=herba_subtree)
## $z
## phi
## 5.099736
##
## $pval
## pval
## 3.401275e-07
###########################
##### woody_subtree ###
###########################
woody_subtree<-read.tree("woody.nwk")
tiplabel<-woody_subtree$tip.label
tiplabel
## [1] "Liriodendron_tulipifera" "Asarum_europaeum"
## [3] "Origanum_vulgare" "Fraxinus_excelsior"
## [5] "Scabiosa_columbaria" "Artemisia_absinthium"
## [7] "Inula_helenium" "Tilia_cordata"
## [9] "Acer_rubrum" "Aesculus_hippocastanum"
## [11] "Populus_alba" "Salix_alba"
## [13] "Hypericum_perforatum" "Corylus_avellana"
## [15] "Betula_pendula" "Juglans_regia"
## [17] "Carya_ovata" "Quercus_robur"
## [19] "Castanea_sativa" "Fagus_sylvatica"
## [21] "Malus_domestica" "Gleditsia_triacanthos"
## [23] "Baptisia_australis" "Buxus_sempervirens"
## [25] "Ranunculus_acris" "Clematis_vitalba"
## [27] "Hepatica_nobilis" "Helleborus_niger"
## [29] "Actaea_racemosa" "Picea_abies"
## [31] "Pinus_sylvestris" "Larix_decidua"
length(tiplabel)
## [1] 32
tiplabel <- gsub("_", " ", tiplabel)
woody_subtree$tip.label <- tiplabel
tiplabel
## [1] "Liriodendron tulipifera" "Asarum europaeum"
## [3] "Origanum vulgare" "Fraxinus excelsior"
## [5] "Scabiosa columbaria" "Artemisia absinthium"
## [7] "Inula helenium" "Tilia cordata"
## [9] "Acer rubrum" "Aesculus hippocastanum"
## [11] "Populus alba" "Salix alba"
## [13] "Hypericum perforatum" "Corylus avellana"
## [15] "Betula pendula" "Juglans regia"
## [17] "Carya ovata" "Quercus robur"
## [19] "Castanea sativa" "Fagus sylvatica"
## [21] "Malus domestica" "Gleditsia triacanthos"
## [23] "Baptisia australis" "Buxus sempervirens"
## [25] "Ranunculus acris" "Clematis vitalba"
## [27] "Hepatica nobilis" "Helleborus niger"
## [29] "Actaea racemosa" "Picea abies"
## [31] "Pinus sylvestris" "Larix decidua"
woody_Gsize<-read.csv("woody.csv")
#length(diurnal_BodyMass$Species)
#diurnal_BodyMass <- diurnal_BodyMass[!is.na(diurnal_BodyMass$BodyMassMale_kg), ]
same<-intersect(tiplabel,woody_Gsize$Species)
woody_Gsize <- woody_Gsize[woody_Gsize$Species %in% same, ]
# 更新樹中的物種列表
woody_subtree <- drop.tip(woody_subtree, setdiff(tiplabel, same))
tiplabel <- woody_subtree$tip.label
# 確保樹和體重數據中物種一致
same <- intersect(tiplabel, woody_Gsize$Species)
woody_Gsize <- woody_Gsize[woody_Gsize$Species %in% same, ]
#隨機刪除20個物種
#tips_to_drop <- sample(tiplabel, 20)
#diurnal_subtree <- drop.tip(diurnal_subtree, tips_to_drop)
#tiplabel <- diurnal_subtree$tip.label
# 更新體重數據中的物種列表
#same <- intersect(tiplabel, diurnal_BodyMass$Species)
#diurnal_BodyMass <- diurnal_BodyMass[diurnal_BodyMass$Species %in% same, ]
# 6. 繪製更新後的樹
plot(woody_subtree)
y<-woody_Gsize$GenomeSize
names(y)<-woody_Gsize$Species
# 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)
#y<-diurnal_BodyMass$BodyMassMale_kg
#names(y)<-diurnal_BodyMass$Species
#######################################
############ diurnal Analysis #######
#######################################
mfit<-modelfit(trait=y, tree=woody_subtree)
woody.empmodelout<-empmodel(mfit)
woody.empmodelout$parammleset
## phi phi2 theta theta2 sigsq
## AR1 -0.01322364 NA NA NA 0.002350349
## AR2 -0.01739255 -0.30212951 NA NA 0.002638185
## ARMA1 0.35852288 NA -0.44616675 NA 0.002413481
## ARMA21 -0.26877030 -0.31466190 0.28250260 NA 0.002574759
## ARMA22 -0.22288296 -0.71387980 0.20666883 0.450465790 0.002497982
## Weighted Estimate 0.02314866 -0.01028542 -0.04418663 0.004504658 0.002337511
woody.empmodelout$fitset
## modelset paranumset loglikset aicset aicweightset rankset
## 1 AR1 2 98.39 -192.77 0.86 1
## 2 AR2 3 93.44 -180.88 0.00 5
## 3 ARMA1 3 97.36 -188.71 0.11 2
## 4 ARMA21 4 95.96 -183.91 0.01 4
## 5 ARMA22 5 97.16 -184.31 0.01 3
#empmodelout$parammleset[1,][c("phi","sigsq")]
xtable(woody.empmodelout$parammleset,digits=3)
## % latex table generated in R 4.4.0 by xtable 1.8-4 package
## % Sun Sep 29 21:42:10 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrr}
## \hline
## & phi & phi2 & theta & theta2 & sigsq \\
## \hline
## AR1 & -0.013 & & & & 0.002 \\
## AR2 & -0.017 & -0.302 & & & 0.003 \\
## ARMA1 & 0.359 & & -0.446 & & 0.002 \\
## ARMA21 & -0.269 & -0.315 & 0.283 & & 0.003 \\
## ARMA22 & -0.223 & -0.714 & 0.207 & 0.450 & 0.002 \\
## Weighted Estimate & 0.023 & -0.010 & -0.044 & 0.005 & 0.002 \\
## \hline
## \end{tabular}
## \end{table}
xtable(woody.empmodelout$fitset,digits=3)
## % latex table generated in R 4.4.0 by xtable 1.8-4 package
## % Sun Sep 29 21:42:10 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rlrrrrr}
## \hline
## & modelset & paranumset & loglikset & aicset & aicweightset & rankset \\
## \hline
## 1 & AR1 & 2.000 & 98.390 & -192.770 & 0.860 & 1.000 \\
## 2 & AR2 & 3.000 & 93.440 & -180.880 & 0.000 & 5.000 \\
## 3 & ARMA1 & 3.000 & 97.360 & -188.710 & 0.110 & 2.000 \\
## 4 & ARMA21 & 4.000 & 95.960 & -183.910 & 0.010 & 4.000 \\
## 5 & ARMA22 & 5.000 & 97.160 & -184.310 & 0.010 & 3.000 \\
## \hline
## \end{tabular}
## \end{table}
## add up test.phi here
woody_betas.rrphylo<-compute_betas_rrphylo(woody_subtree, y)$betas.rrphylo
test.phi(param=woody.empmodelout$parammleset["AR1",][c("phi","sigsq")], betas.rrphylo=woody_betas.rrphylo,tree=woody_subtree)
## $z
## phi
## -0.1214592
##
## $pval
## pval
## 0.9033273
####################################
#########Combined tree##############
####################################
tree_with_gsize<-bind.tree(herba_subtree, woody_subtree, where = "root", position = 0, interactive = FALSE)
tree_with_gsize<-chronopl(tree_with_gsize, 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_gsize, 0, age.min = 1, age.max = NULL, node = "root", : at least one internal branch is of length zero:
## it was collapsed and some nodes have been deleted.
plot(tree_with_gsize, show.tip.label = FALSE)
#合併日夜間體重
herba <- data.frame(Species = herba_Gsize$Species, Gsize = herba_Gsize$GenomeSize)
woody <- data.frame(Species = woody_Gsize$Species, Gsize = woody_Gsize$GenomeSize)
Gsize <- rbind(herba, woody)
# 確保合併後的體重數據和合併後的樹的端點名稱一致
common_species <- intersect(tree_with_gsize$tip.label, Gsize$Species)
Gsize <- Gsize[Gsize$Species %in% common_species, ]
y <- Gsize$Gsize
names(y) <- Gsize$Species
y
## Achillea millefolium Aconitum napellus
## 3.50 15.50
## Alchemilla vulgaris Anemone nemorosa
## 5.50 18.30
## Angelica archangelica Aquilegia vulgaris
## 7.20 6.90
## Aster amellus Campanula rotundifolia
## 8.10 2.60
## Convolvulus arvensis Delphinium elatum
## 1.43 13.60
## Digitalis purpurea Echinacea purpurea
## 3.89 8.30
## Eupatorium cannabinum Filipendula ulmaria
## 4.00 3.00
## Fragaria vesca Galium aparine
## 0.63 2.50
## Geranium robertianum Geum rivale
## 3.00 3.20
## Lathyrus vernus Lupinus perennis
## 6.40 2.15
## Mentha longifolia Nepeta cataria
## 2.40 1.70
## Pulsatilla vulgaris Rudbeckia hirta
## 6.80 2.30
## Salvia pratensis Sanguisorba officinalis
## 2.50 4.60
## Solidago virgaurea Symphytum officinale
## 2.40 3.70
## Tanacetum vulgare Thalictrum aquilegiifolium
## 3.20 8.30
## Valeriana officinalis Verbascum thapsus
## 3.00 1.58
## Actaea racemosa Artemisia absinthium
## 6.20 5.40
## Asarum europaeum Baptisia australis
## 17.90 2.95
## Clematis vitalba Helleborus niger
## 16.50 12.30
## Hepatica nobilis Hypericum perforatum
## 24.00 1.30
## Inula helenium Origanum vulgare
## 3.20 1.63
## Ranunculus acris Scabiosa columbaria
## 11.20 2.50
## Acer rubrum Aesculus hippocastanum
## 1.14 1.44
## Betula pendula Buxus sempervirens
## 0.88 1.56
## Carya ovata Castanea sativa
## 1.30 1.75
## Corylus avellana Fagus sylvatica
## 0.88 1.14
## Fraxinus excelsior Gleditsia triacanthos
## 1.25 1.38
## Juglans regia Larix decidua
## 1.23 0.78
## Liriodendron tulipifera Malus domestica
## 2.30 1.57
## Picea abies Pinus sylvestris
## 1.11 1.48
## Populus alba Quercus robur
## 1.23 1.18
## Salix alba Tilia cordata
## 1.05 1.34
# 確認 y 和樹的端點名稱一致
tree_with_gsize$tip.label <- tree_with_gsize$tip.label[tree_with_gsize$tip.label %in% names(y)]
y
## Achillea millefolium Aconitum napellus
## 3.50 15.50
## Alchemilla vulgaris Anemone nemorosa
## 5.50 18.30
## Angelica archangelica Aquilegia vulgaris
## 7.20 6.90
## Aster amellus Campanula rotundifolia
## 8.10 2.60
## Convolvulus arvensis Delphinium elatum
## 1.43 13.60
## Digitalis purpurea Echinacea purpurea
## 3.89 8.30
## Eupatorium cannabinum Filipendula ulmaria
## 4.00 3.00
## Fragaria vesca Galium aparine
## 0.63 2.50
## Geranium robertianum Geum rivale
## 3.00 3.20
## Lathyrus vernus Lupinus perennis
## 6.40 2.15
## Mentha longifolia Nepeta cataria
## 2.40 1.70
## Pulsatilla vulgaris Rudbeckia hirta
## 6.80 2.30
## Salvia pratensis Sanguisorba officinalis
## 2.50 4.60
## Solidago virgaurea Symphytum officinale
## 2.40 3.70
## Tanacetum vulgare Thalictrum aquilegiifolium
## 3.20 8.30
## Valeriana officinalis Verbascum thapsus
## 3.00 1.58
## Actaea racemosa Artemisia absinthium
## 6.20 5.40
## Asarum europaeum Baptisia australis
## 17.90 2.95
## Clematis vitalba Helleborus niger
## 16.50 12.30
## Hepatica nobilis Hypericum perforatum
## 24.00 1.30
## Inula helenium Origanum vulgare
## 3.20 1.63
## Ranunculus acris Scabiosa columbaria
## 11.20 2.50
## Acer rubrum Aesculus hippocastanum
## 1.14 1.44
## Betula pendula Buxus sempervirens
## 0.88 1.56
## Carya ovata Castanea sativa
## 1.30 1.75
## Corylus avellana Fagus sylvatica
## 0.88 1.14
## Fraxinus excelsior Gleditsia triacanthos
## 1.25 1.38
## Juglans regia Larix decidua
## 1.23 0.78
## Liriodendron tulipifera Malus domestica
## 2.30 1.57
## Picea abies Pinus sylvestris
## 1.11 1.48
## Populus alba Quercus robur
## 1.23 1.18
## Salix alba Tilia cordata
## 1.05 1.34
#diurnal<-rbind(diurnal_BodyMass$Species,diurnal_BodyMass$BodyMassMale_kg)
#nocturnal<-rbind(remaining_bodymass$Species,remaining_bodymass$BodyMassMale_kg)
#diurnal<-as.data.frame(diurnal)
#nocturnal<-as.data.frame(nocturnal)
#BodyMass<-cbind(diurnal,nocturnal)
# 將合併後的數據框轉置
#BodyMass<- t(BodyMass)
#BodyMass<-as.data.frame(BodyMass)
y
## Achillea millefolium Aconitum napellus
## 3.50 15.50
## Alchemilla vulgaris Anemone nemorosa
## 5.50 18.30
## Angelica archangelica Aquilegia vulgaris
## 7.20 6.90
## Aster amellus Campanula rotundifolia
## 8.10 2.60
## Convolvulus arvensis Delphinium elatum
## 1.43 13.60
## Digitalis purpurea Echinacea purpurea
## 3.89 8.30
## Eupatorium cannabinum Filipendula ulmaria
## 4.00 3.00
## Fragaria vesca Galium aparine
## 0.63 2.50
## Geranium robertianum Geum rivale
## 3.00 3.20
## Lathyrus vernus Lupinus perennis
## 6.40 2.15
## Mentha longifolia Nepeta cataria
## 2.40 1.70
## Pulsatilla vulgaris Rudbeckia hirta
## 6.80 2.30
## Salvia pratensis Sanguisorba officinalis
## 2.50 4.60
## Solidago virgaurea Symphytum officinale
## 2.40 3.70
## Tanacetum vulgare Thalictrum aquilegiifolium
## 3.20 8.30
## Valeriana officinalis Verbascum thapsus
## 3.00 1.58
## Actaea racemosa Artemisia absinthium
## 6.20 5.40
## Asarum europaeum Baptisia australis
## 17.90 2.95
## Clematis vitalba Helleborus niger
## 16.50 12.30
## Hepatica nobilis Hypericum perforatum
## 24.00 1.30
## Inula helenium Origanum vulgare
## 3.20 1.63
## Ranunculus acris Scabiosa columbaria
## 11.20 2.50
## Acer rubrum Aesculus hippocastanum
## 1.14 1.44
## Betula pendula Buxus sempervirens
## 0.88 1.56
## Carya ovata Castanea sativa
## 1.30 1.75
## Corylus avellana Fagus sylvatica
## 0.88 1.14
## Fraxinus excelsior Gleditsia triacanthos
## 1.25 1.38
## Juglans regia Larix decidua
## 1.23 0.78
## Liriodendron tulipifera Malus domestica
## 2.30 1.57
## Picea abies Pinus sylvestris
## 1.11 1.48
## Populus alba Quercus robur
## 1.23 1.18
## Salix alba Tilia cordata
## 1.05 1.34
length(y)
## [1] 64
names(y)%in% tree_with_gsize$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
tree_with_gsize$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
#######################################
############ Full Analysis ##########
#######################################
#compute_betas_rrphylo(tree_with_bodymass, y)$betas.rrphylo
mfit<-modelfit(trait=y, tree=tree_with_gsize)
comb.empmodelout<-empmodel(mfit)
comb.empmodelout$parammleset
## phi phi2 theta theta2 sigsq
## AR1 0.09504789 NA NA NA 4.731867
## AR2 0.10806202 -0.08314429 NA NA 4.748257
## ARMA1 -0.60026037 NA 0.7380948 NA 4.663459
## ARMA21 -0.04633153 -0.55393033 0.7380966 NA 4.663451
## ARMA22 -0.76858571 0.16832434 0.8690300 -0.130934718 4.663469
## Weighted Estimate -0.05037934 -0.07947511 0.2019044 -0.002618694 4.722083
comb.empmodelout$fitset
## modelset paranumset loglikset aicset aicweightset rankset
## 1 AR1 2 -267.92 539.85 0.20 2
## 2 AR2 3 -265.94 537.87 0.53 1
## 3 ARMA1 3 -267.04 540.07 0.18 3
## 4 ARMA21 4 -267.04 542.07 0.07 4
## 5 ARMA22 5 -267.04 544.07 0.02 5
xtable(comb.empmodelout$parammleset,digits=3)
## % latex table generated in R 4.4.0 by xtable 1.8-4 package
## % Sun Sep 29 21:42:16 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrr}
## \hline
## & phi & phi2 & theta & theta2 & sigsq \\
## \hline
## AR1 & 0.095 & & & & 4.732 \\
## AR2 & 0.108 & -0.083 & & & 4.748 \\
## ARMA1 & -0.600 & & 0.738 & & 4.663 \\
## ARMA21 & -0.046 & -0.554 & 0.738 & & 4.663 \\
## ARMA22 & -0.769 & 0.168 & 0.869 & -0.131 & 4.663 \\
## Weighted Estimate & -0.050 & -0.079 & 0.202 & -0.003 & 4.722 \\
## \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 Sep 29 21:42:16 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rlrrrrr}
## \hline
## & modelset & paranumset & loglikset & aicset & aicweightset & rankset \\
## \hline
## 1 & AR1 & 2.000 & -267.920 & 539.850 & 0.200 & 2.000 \\
## 2 & AR2 & 3.000 & -265.940 & 537.870 & 0.530 & 1.000 \\
## 3 & ARMA1 & 3.000 & -267.040 & 540.070 & 0.180 & 3.000 \\
## 4 & ARMA21 & 4.000 & -267.040 & 542.070 & 0.070 & 4.000 \\
## 5 & ARMA22 & 5.000 & -267.040 & 544.070 & 0.020 & 5.000 \\
## \hline
## \end{tabular}
## \end{table}
## add up test.phi here
comb_betas.rrphylo<-compute_betas_rrphylo(tree_with_gsize, y)$betas.rrphylo
test.phi(param=comb.empmodelout$parammleset["AR1",][c("phi","sigsq")], betas.rrphylo=comb_betas.rrphylo,tree=tree_with_gsize)
## $z
## phi
## 1.119007
##
## $pval
## pval
## 0.2631373
##############
####AR1 test##
##############
#########Combined tree##############
comb_betas.rrphylo<-compute_betas_rrphylo(tree_with_gsize, y)$betas.rrphylo
combAr1test <- test.phi(param=comb.empmodelout$parammleset["AR1",][c("phi","sigsq")], betas.rrphylo=comb_betas.rrphylo,tree=tree_with_gsize)
#########herba tree##############
herba_betas.rrphylo<-compute_betas_rrphylo(herba_subtree, y)$betas.rrphylo
herbaAr1test <- test.phi(param=herba.empmodelout$parammleset["AR1",][c("phi","sigsq")], betas.rrphylo=herba_betas.rrphylo,tree=herba_subtree)
#########nocturnal tree##############
woody_betas.rrphylo<-compute_betas_rrphylo(woody_subtree, y)$betas.rrphylo
woodyAr1test <- test.phi(param=woody.empmodelout$parammleset["AR1",][c("phi","sigsq")], betas.rrphylo=woody_betas.rrphylo,tree=woody_subtree)
treecase<-c("combined tree", "herba tree", "woody tree")
zval<-c(combAr1test$z,herbaAr1test$z,woodyAr1test$z)
pval<-c(combAr1test$pval,herbaAr1test$pval,woodyAr1test$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 1.1190067 2.631373e-01 N
## herba tree 5.0997360 3.401275e-07 Y
## woody tree -0.1214592 9.033273e-01 N
xtable(ar1testtb,digits=3)
## % latex table generated in R 4.4.0 by xtable 1.8-4 package
## % Sun Sep 29 21:42:19 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrl}
## \hline
## & zval & pval & sig \\
## \hline
## combined tree & 1.119 & 0.263 & N \\
## herba tree & 5.100 & 0.000 & Y \\
## woody tree & -0.121 & 0.903 & N \\
## \hline
## \end{tabular}
## \end{table}
# 將日間樹、夜間樹的端點名稱轉換為縮寫形式
# herba_subtree$tip.label <- sapply(herba_subtree$tip.label, function(species) {
# parts <- strsplit(species, "_")[[1]]
# paste0(substr(parts[1], 1, 1), ".", parts[2])
# })
#
# woody_subtree$tip.label <- sapply(woody_subtree$tip.label, function(species) {
# parts <- strsplit(species, "_")[[1]]
# paste0(substr(parts[1], 1, 1), ".", parts[2])
# })
#plot(tree_with_bodymass, show.tip.label = FALSE)
#par(mfrow = c(1, 3), mar = c(1, 1, 1, 1))
#plot(diurnal_subtree,main="Diurnal 12")
#plot(tree_with_bodymass, show.tip.label = FALSE)
#plot(nocturnal_subtree_with_bodymass,main="Nocturnal 11")
par(mfrow = c(1, 3), mar = c(1, 1, 1, 1))
plot(herba_subtree, main="herba species",cex=1.10)
plot(tree_with_gsize, show.tip.label = FALSE,main="(herba+woody) species")
plot(woody_subtree, main="Woody species",cex=1.10)
###可不用
#################################
#### LRT across models ##
#################################
LL_herba <- herba.empmodelout$fitset
LL_woody <- woody.empmodelout$fitset
LL_comb <- comb.empmodelout$fitset
# Compute the test statistic
test_statistic <- 2 * ((LL_herba$loglikset + LL_woody$loglikset)-LL_comb$loglikset)
test_statistic
## [1] 847.16 830.44 843.34 840.50 842.94
#names(test_statistic)<-c(paste(LL_comb$modelset,".chi2",sep=""))
test_statistic
## [1] 847.16 830.44 843.34 840.50 842.94
# 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_herba$loglikset
## [1] 57.27 55.84 57.27 57.25 57.27
LL_woody$loglikset
## [1] 98.39 93.44 97.36 95.96 97.16
LL_comb$loglikset
## [1] -267.92 -265.94 -267.04 -267.04 -267.04
mset<-LL_comb$modelset
herbalikeset<-LL_herba$loglikset
woodylikeset<-LL_woody$loglikset
comblikeset<-LL_comb$loglikset
df2rates1rate<-data.frame(
herbalikeset=herbalikeset,
woodylikeset=woodylikeset,
comblikeset=comblikeset,
test_statistic=test_statistic,
p_value.array=p_value.array,tworates=tworates
)
rownames(df2rates1rate)<-LL_herba$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 57.27 98.39 -267.92 847.16 0 Y
## AR2 55.84 93.44 -265.94 830.44 0 Y
## ARMA1 57.27 97.36 -267.04 843.34 0 Y
## ARMA21 57.25 95.96 -267.04 840.50 0 Y
## ARMA22 57.27 97.16 -267.04 842.94 0 Y
xtable(df2rates1rate)
## % latex table generated in R 4.4.0 by xtable 1.8-4 package
## % Sun Sep 29 21:42:19 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrrl}
## \hline
## & log L\_di & log L\_no & log L\_dino & chi2 & pval & NA \\
## \hline
## AR1 & 57.27 & 98.39 & -267.92 & 847.16 & 0.00 & Y \\
## AR2 & 55.84 & 93.44 & -265.94 & 830.44 & 0.00 & Y \\
## ARMA1 & 57.27 & 97.36 & -267.04 & 843.34 & 0.00 & Y \\
## ARMA21 & 57.25 & 95.96 & -267.04 & 840.50 & 0.00 & Y \\
## ARMA22 & 57.27 & 97.16 & -267.04 & 842.94 & 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)