rm(list=ls())
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
ntaxa<-8
stree<-stree(ntaxa)
tree<-compute.brlen(stree(ntaxa,type="balanced"))
#tree<-reorder(tree,"postorder")
fastBM(tree,internal=T,sig2=1)->phen
phen[1:ntaxa]->y
trait<-round(y,2)
trait
## t1 t2 t3 t4 t5 t6 t7 t8
## -1.14 -1.43 -1.23 -1.64 -0.39 -0.95 0.77 0.42
edgelength<-tree$edge.length
#names(edgelength)<-edge_labels
#edgelength
betas.rrphylo <- compute_betas_rrphylo(tree, y)$betas.rrphylo
betas.rrphylo
## [,1]
## 9 0.1562941331
## 10 -0.1857322335
## 11 -0.0214344439
## 12 -0.0714316728
## 13 0.8109087659
## 14 -0.0054097547
## 15 0.4108641377
## 1 0.0081981086
## 2 -0.0189153306
## 3 0.0004793459
## 4 -0.0361951823
## 5 0.0241368545
## 6 -0.0268417318
## 7 0.1189811462
## 8 0.0864509226
# Rename node labels
node_labels <- paste("nd", (ntaxa+1):(2*ntaxa-1),sep="")
plot(tree, show.tip.label = FALSE , edge.width = 2, cex = 3)
title(main = "Raw Tree", cex.main = 1.5)
nodelabels(node_labels, cex = 1.5)
# Add tip labels
tiplabels(paste("y",1:ntaxa,sep=""), cex = 1.5)
# Get the number of edges
num_edges <- Nedge(tree)
# Create edge labels
edge_labels <- paste0("eg", seq_len(num_edges),sep="")
edgelabels(edge_labels, cex = 1.5)
df.trait<-data.frame(trait=t(trait))
colnames(df.trait)<-paste("y",1:length(trait),sep="")
df.trait
## y1 y2 y3 y4 y5 y6 y7 y8
## 1 -1.14 -1.43 -1.23 -1.64 -0.39 -0.95 0.77 0.42
xtable(df.trait)
## % latex table generated in R 4.3.2 by xtable 1.8-4 package
## % Tue Apr 30 08:43:11 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrrrrr}
## \hline
## & y1 & y2 & y3 & y4 & y5 & y6 & y7 & y8 \\
## \hline
## 1 & -1.14 & -1.43 & -1.23 & -1.64 & -0.39 & -0.95 & 0.77 & 0.42 \\
## \hline
## \end{tabular}
## \end{table}
df.edgelen<-data.frame(edgelength=round(t(edgelength),3))
colnames(df.edgelen)<-paste("eg",1:length(df.edgelen),sep="")
df.edgelen
## eg1 eg2 eg3 eg4 eg5 eg6 eg7 eg8 eg9 eg10 eg11 eg12 eg13
## 1 0.571 0.286 0.143 0.143 0.286 0.143 0.143 0.571 0.286 0.143 0.143 0.286 0.143
## eg14
## 1 0.143
xtable(df.edgelen,digits=3)
## % latex table generated in R 4.3.2 by xtable 1.8-4 package
## % Tue Apr 30 08:43:11 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrrrrrrrrrrr}
## \hline
## & eg1 & eg2 & eg3 & eg4 & eg5 & eg6 & eg7 & eg8 & eg9 & eg10 & eg11 & eg12 & eg13 & eg14 \\
## \hline
## 1 & 0.571 & 0.286 & 0.143 & 0.143 & 0.286 & 0.143 & 0.143 & 0.571 & 0.286 & 0.143 & 0.143 & 0.286 & 0.143 & 0.143 \\
## \hline
## \end{tabular}
## \end{table}
df.betas.rrphylo<-data.frame(betas.rrphylo=t(round(betas.rrphylo,3)))
colnames(df.betas.rrphylo)<-paste("nd",1:length(df.betas.rrphylo),sep="")
df.betas.rrphylo
## nd1 nd2 nd3 nd4 nd5 nd6 nd7 nd8 nd9 nd10 nd11 nd12
## 1 0.156 -0.186 -0.021 -0.071 0.811 -0.005 0.411 0.008 -0.019 0 -0.036 0.024
## nd13 nd14 nd15
## 1 -0.027 0.119 0.086
xtable(df.betas.rrphylo,digits=3)
## % latex table generated in R 4.3.2 by xtable 1.8-4 package
## % Tue Apr 30 08:43:11 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrrrrrrrrrrrr}
## \hline
## & nd1 & nd2 & nd3 & nd4 & nd5 & nd6 & nd7 & nd8 & nd9 & nd10 & nd11 & nd12 & nd13 & nd14 & nd15 \\
## \hline
## 1 & 0.156 & -0.186 & -0.021 & -0.071 & 0.811 & -0.005 & 0.411 & 0.008 & -0.019 & 0.000 & -0.036 & 0.024 & -0.027 & 0.119 & 0.086 \\
## \hline
## \end{tabular}
## \end{table}
#######################################
############ Formal Analysis ##########
#######################################
mfit<-modelfit(trait=y, tree=tree)
empmodelout<-empmodel(mfit)
empmodelout$parammleset
## phi phi2 theta theta2 sigsq
## AR1 0.29798390 NA NA NA 0.05639068
## AR2 0.30673460 -0.01756866 NA NA 0.06102567
## ARMA1 0.58670030 NA -0.40926798 NA 0.05434437
## ARMA21 -0.05687953 0.64357982 -0.40926795 NA 0.05434437
## ARMA22 -0.02356677 0.61028222 -0.02426522 -0.38503125 0.05434445
## Weighted Estimate 0.32674169 0.07412241 -0.13578639 -0.01155094 0.05621021
empmodelout$fitset
## modelset paranumset loglikset aicset aicweightset rankset
## 1 AR1 2 0.25 3.51 0.52 1
## 2 AR2 3 -0.25 6.49 0.12 3
## 3 ARMA1 3 0.49 5.03 0.24 2
## 4 ARMA21 4 0.49 7.03 0.09 4
## 5 ARMA22 5 0.49 9.03 0.03 5
xtable(empmodelout$parammleset,digits=3)
## % latex table generated in R 4.3.2 by xtable 1.8-4 package
## % Tue Apr 30 08:43:12 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrr}
## \hline
## & phi & phi2 & theta & theta2 & sigsq \\
## \hline
## AR1 & 0.298 & & & & 0.056 \\
## AR2 & 0.307 & -0.018 & & & 0.061 \\
## ARMA1 & 0.587 & & -0.409 & & 0.054 \\
## ARMA21 & -0.057 & 0.644 & -0.409 & & 0.054 \\
## ARMA22 & -0.024 & 0.610 & -0.024 & -0.385 & 0.054 \\
## Weighted Estimate & 0.327 & 0.074 & -0.136 & -0.012 & 0.056 \\
## \hline
## \end{tabular}
## \end{table}
xtable(empmodelout$fitset,digits=3)
## % latex table generated in R 4.3.2 by xtable 1.8-4 package
## % Tue Apr 30 08:43:12 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rlrrrrr}
## \hline
## & modelset & paranumset & loglikset & aicset & aicweightset & rankset \\
## \hline
## 1 & AR1 & 2.000 & 0.250 & 3.510 & 0.520 & 1.000 \\
## 2 & AR2 & 3.000 & -0.250 & 6.490 & 0.120 & 3.000 \\
## 3 & ARMA1 & 3.000 & 0.490 & 5.030 & 0.240 & 2.000 \\
## 4 & ARMA21 & 4.000 & 0.490 & 7.030 & 0.090 & 4.000 \\
## 5 & ARMA22 & 5.000 & 0.490 & 9.030 & 0.030 & 5.000 \\
## \hline
## \end{tabular}
## \end{table}
testrate2subtrees_output <- testrate2subtrees(tree, y)
## [1] "t10" "t11" "t12"
## [1] "5" "6" "7" "1" "2" "3" "4"
## [1] "t13" "t14" "t15"
## [1] "5" "6" "7" "1" "2" "3" "4"
tworatedf<- treeARanova(testrate2subtrees_output)
tworatedf
## numparam mlesoutphi mlesoutsigsq likelihoodsout teststat deg
## treeAR1 2 0.298 0.052 -0.227 16.111 2
## subtree1AR1 2 0.472 0.048 -0.006 NA NA
## subtree2AR1 2 0.300 1.101 7.834 NA NA
## pval
## treeAR1 3.172931e-04
## subtree1AR1 <NA>
## subtree2AR1 <NA>
xtable(tworatedf,digits=3)
## % latex table generated in R 4.3.2 by xtable 1.8-4 package
## % Tue Apr 30 08:43:12 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrrrl}
## \hline
## & numparam & mlesoutphi & mlesoutsigsq & likelihoodsout & teststat & deg & pval \\
## \hline
## treeAR1 & 2.000 & 0.298 & 0.052 & -0.227 & 16.111 & 2.000 & 3.172931e-04 \\
## subtree1AR1 & 2.000 & 0.472 & 0.048 & -0.006 & & & \\
## subtree2AR1 & 2.000 & 0.300 & 1.101 & 7.834 & & & \\
## \hline
## \end{tabular}
## \end{table}