library(xtable)
library(ape)
# analyze salamander data under four models
# read data
phy<-read.tree("https://tonyjhwueng.info/phyrates/ple.nwk")
plot(phy)
df<-read.csv("https://tonyjhwueng.info/phyrates/Adams2012-SystBiolData.csv")
head(df)
## X HeadLength BodyWidth Forelimb
## 1 P_albagula 12.320000 8.160000 16.12000
## 2 P_amplus 11.120000 7.580000 16.44000
## 3 P_angusticlavius 7.033333 3.700000 8.10000
## 4 P_aureolus 9.200000 5.314286 12.52857
## 5 P_caddoensis 8.200000 3.600000 9.20000
## 6 P_chattahoochee 10.371429 6.400000 13.67143
spX<-strsplit(as.character(df$X),"_")
spname<-array(NA,length(phy$tip.label))
for(Index in 1:length(phy$tip.label)){
spname[Index]<-paste("Plethodon_", spX[[Index]][2],sep="")
}
spname
## [1] "Plethodon_albagula" "Plethodon_amplus"
## [3] "Plethodon_angusticlavius" "Plethodon_aureolus"
## [5] "Plethodon_caddoensis" "Plethodon_chattahoochee"
## [7] "Plethodon_cheoah" "Plethodon_chlorobryonis"
## [9] "Plethodon_cinereus" "Plethodon_cylindraceus"
## [11] "Plethodon_dorsalis" "Plethodon_electromorphus"
## [13] "Plethodon_fourchensis" "Plethodon_glutinosus"
## [15] "Plethodon_grobmani" "Plethodon_hoffmani"
## [17] "Plethodon_hubrichti" "Plethodon_jordani"
## [19] "Plethodon_kentucki" "Plethodon_kiamichi"
## [21] "Plethodon_kisatchie" "Plethodon_meridianus"
## [23] "Plethodon_metcalfi" "Plethodon_mississippi"
## [25] "Plethodon_montanus" "Plethodon_nettingi"
## [27] "Plethodon_ocmulgee" "Plethodon_ouachitae"
## [29] "Plethodon_petraeus" "Plethodon_punctatus"
## [31] "Plethodon_richmondi" "Plethodon_savannah"
## [33] "Plethodon_sequoyah" "Plethodon_serratus"
## [35] "Plethodon_shenandoah" "Plethodon_shermani"
## [37] "Plethodon_teyahalee" "Plethodon_variolatus"
## [39] "Plethodon_ventralis" "Plethodon_virginia"
## [41] "Plethodon_websteri" "Plethodon_wehrlei"
## [43] "Plethodon_welleri" "Plethodon_yonahlossee"
df$X<-spname
HeadLength<-df$HeadLength
names(HeadLength)<-spname
BodyWidth<-df$BodyWidth
names(BodyWidth)<-spname
HeadLength<-HeadLength[phy$tip.label]
BodyWidth<-BodyWidth[phy$tip.label]
x<-cbind(HeadLength,BodyWidth)
# BM
source("https://tonyjhwueng.info/phyrates/ComparingBM.r")
## Loading required package: maps
result.bm<-CompareRates.multTrait(phy=phy,x=x,TraitCov=T,ms.err=NULL,ms.cov=NULL)
# PMM
source("https://tonyjhwueng.info/phyrates/ComparingPMM.r")
result.pmm<-CompareRates.multTrait(phy=phy,x=x,TraitCov=T,ms.err=NULL,ms.cov=NULL)
# OU
source("https://tonyjhwueng.info/phyrates/ComparingOU.r")
result.ou<-CompareRates.multTrait(phy=phy,x=x,TraitCov=T,ms.err=NULL,ms.cov=NULL)
# EB
source("https://tonyjhwueng.info/phyrates/ComparingEB.r")
result.eb<-CompareRates.multTrait(phy=phy,x=x,TraitCov=T,ms.err=NULL,ms.cov=NULL)
## Warning in fitContinuous(phy = phy, dat = x[, Index], model = "EB"):
## Parameter estimates appear at bounds:
## a
## Warning in fitContinuous(phy = phy, dat = x[, Index], model = "EB"):
## Parameter estimates appear at bounds:
## a
print(result.bm)
## $Robs
## HeadLength BodyWidth
## HeadLength 0.2282753 0.1252281
## BodyWidth 0.1252281 0.1192467
##
## $Rconstrained
## [,1] [,2]
## [1,] 1.0154436 0.1252281
## [2,] 0.1252281 1.0154436
##
## $Lobs
## [1] -126.8484
##
## $Lconstrained
## [1] -188.316
##
## $LRTest
## [1] 122.9352
##
## $Prob
## [1] 1.440774e-28
##
## $AICc.obs
## [1] 261.6968
##
## $AICc.constrained
## [1] 382.632
##
## $optimmessage
## [1] "Optimization has converged."
print(result.pmm)
## $Robs.pmm
## [,1] [,2]
## [1,] 0.2222583 0.1252084
## [2,] 0.1252084 0.1251008
##
## $Rconstrained.pmm
## [,1] [,2]
## [1,] 1.0166341 0.1300416
## [2,] 0.1300416 1.0166341
##
## $Lobs.pmm
## [1] -126.824
##
## $Lconstrained.pmm
## [1] -187.8177
##
## $LRTest.pmm
## [1] 121.9874
##
## $Prob.pmm
## [1] 2.323066e-28
##
## $AICc.obs.pmm
## [1] 265.648
##
## $AICc.constrained.pmm
## [1] 385.6354
##
## $optimmessage.pmm
## [1] "Optimization has converged."
print(result.ou)
## $Robs.ou
## [,1] [,2]
## [1,] 0.2343296 0.1322438
## [2,] 0.1322438 0.1302819
##
## $Rconstrained.ou
## [,1] [,2]
## [1,] 1.0243149 0.1577517
## [2,] 0.1577517 1.0243149
##
## $Lobs.ou
## [1] -126.3931
##
## $Lconstrained.ou
## [1] -185.9488
##
## $LRTest.ou
## [1] 119.1114
##
## $Prob.ou
## [1] 9.901163e-28
##
## $AICc.obs.ou
## [1] 264.7862
##
## $AICc.constrained.ou
## [1] 381.8976
##
## $optimmessage.ou
## [1] "Optimization has converged."
print(result.eb)
## $Robs.eb
## [,1] [,2]
## [1,] 0.2282800 0.1238533
## [2,] 0.1238533 0.1192491
##
## $Rconstrained.eb
## [,1] [,2]
## [1,] 1.0154639 0.1252306
## [2,] 0.1252306 1.0154639
##
## $Lobs.eb
## [1] -126.8613
##
## $Lconstrained.eb
## [1] -188.316
##
## $LRTest.eb
## [1] 122.9094
##
## $Prob.eb
## [1] 1.459645e-28
##
## $AICc.obs.eb
## [1] 265.7226
##
## $AICc.constrained.eb
## [1] 386.632
##
## $optimmessage.eb
## [1] "Optimization has converged."
####### BM ########
##OBS
bm.obs.table<-array(NA,c(2,8))
rownames(bm.obs.table)<-c("HL","BW")
colnames(bm.obs.table)<-
c("","","Lobs","Lconstrained","LRTest",
"Prob","AICc.obs","AICc.constrained")
bm.obs.table[1:2,1:2]<-round(result.bm$Robs,3)
bm.obs.table[1,"Lobs"]<-round(result.bm$Lobs,3)
bm.obs.table
## Lobs Lconstrained LRTest Prob AICc.obs AICc.constrained
## HL 0.228 0.125 -126.848 NA NA NA NA NA
## BW 0.125 0.119 NA NA NA NA NA NA
## COMMON
bm.com.table<-array(NA,c(2,8))
rownames(bm.com.table)<-c("HL","BW")
colnames(bm.com.table)<-
c("","","Lobs","Lconstrained","LRTest",
"Prob","AICc.obs","AICc.constrained")
bm.com.table[1:2,1:2]<-round(result.bm$Rconstrained,3)
bm.com.table[1,"Lconstrained"]<-round(result.bm$Lconstrained,3)
bm.com.table[1,"LRTest"]<-round(result.bm$LRTest,3)
bm.com.table[1,"Prob"]<-round(result.bm$Prob,3)
bm.com.table[1,"AICc.obs"]<-round(result.bm$AICc.obs,3)
bm.com.table[1,"AICc.constrained"]<-round(result.bm$AICc.constrained,3)
bm.com.table
## Lobs Lconstrained LRTest Prob AICc.obs AICc.constrained
## HL 1.015 0.125 NA -188.316 122.935 0 261.697 382.632
## BW 0.125 1.015 NA NA NA NA NA NA
xtable(rbind(bm.obs.table,array(NA,c(1,8)),bm.com.table))
## % latex table generated in R 4.0.2 by xtable 1.8-4 package
## % Wed Dec 23 22:40:08 2020
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrrrrr}
## \hline
## & V1 & V2 & Lobs & Lconstrained & LRTest & Prob & AICc.obs & AICc.constrained \\
## \hline
## HL & 0.23 & 0.12 & -126.85 & & & & & \\
## BW & 0.12 & 0.12 & & & & & & \\
## X & & & & & & & & \\
## HL.1 & 1.01 & 0.12 & & -188.32 & 122.94 & 0.00 & 261.70 & 382.63 \\
## BW.1 & 0.12 & 1.01 & & & & & & \\
## \hline
## \end{tabular}
## \end{table}
####### PMM ########
##OBS
pmm.obs.table<-array(NA,c(2,8))
rownames(pmm.obs.table)<-c("HL","BW")
colnames(pmm.obs.table)<-
c("","","Lobs","Lconstrained","LRTest",
"Prob","AICc.obs","AICc.constrained")
pmm.obs.table[1:2,1:2]<-round(result.pmm$Robs,3)
pmm.obs.table[1,"Lobs"]<-round(result.pmm$Lobs,3)
pmm.obs.table
## Lobs Lconstrained LRTest Prob AICc.obs AICc.constrained
## HL 0.222 0.125 -126.824 NA NA NA NA NA
## BW 0.125 0.125 NA NA NA NA NA NA
## COMMON
pmm.com.table<-array(NA,c(2,8))
rownames(pmm.com.table)<-c("HL","BW")
colnames(pmm.com.table)<-
c("","","Lobs","Lconstrained","LRTest",
"Prob","AICc.obs","AICc.constrained")
pmm.com.table[1:2,1:2]<-round(result.pmm$Rconstrained,3)
pmm.com.table[1,"Lconstrained"]<-round(result.pmm$Lconstrained,3)
pmm.com.table[1,"LRTest"]<-round(result.pmm$LRTest,3)
pmm.com.table[1,"Prob"]<-round(result.pmm$Prob,3)
pmm.com.table[1,"AICc.obs"]<-round(result.pmm$AICc.obs,3)
pmm.com.table[1,"AICc.constrained"]<-round(result.pmm$AICc.constrained,3)
pmm.com.table
## Lobs Lconstrained LRTest Prob AICc.obs AICc.constrained
## HL 1.017 0.130 NA -187.818 121.987 0 265.648 385.635
## BW 0.130 1.017 NA NA NA NA NA NA
xtable(rbind(pmm.obs.table,pmm.com.table))
## % latex table generated in R 4.0.2 by xtable 1.8-4 package
## % Wed Dec 23 22:40:08 2020
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrrrrr}
## \hline
## & V1 & V2 & Lobs & Lconstrained & LRTest & Prob & AICc.obs & AICc.constrained \\
## \hline
## HL & 0.22 & 0.12 & -126.82 & & & & & \\
## BW & 0.12 & 0.12 & & & & & & \\
## HL.1 & 1.02 & 0.13 & & -187.82 & 121.99 & 0.00 & 265.65 & 385.63 \\
## BW.1 & 0.13 & 1.02 & & & & & & \\
## \hline
## \end{tabular}
## \end{table}
####### OU ########
##OBS
ou.obs.table<-array(NA,c(2,8))
rownames(ou.obs.table)<-c("HL","BW")
colnames(ou.obs.table)<-
c("","","Lobs","Lconstrained","LRTest",
"Prob","AICc.obs","AICc.constrained")
ou.obs.table[1:2,1:2]<-round(result.ou$Robs,3)
ou.obs.table[1,"Lobs"]<-round(result.ou$Lobs,3)
ou.obs.table
## Lobs Lconstrained LRTest Prob AICc.obs AICc.constrained
## HL 0.234 0.132 -126.393 NA NA NA NA NA
## BW 0.132 0.130 NA NA NA NA NA NA
## COMMON
ou.com.table<-array(NA,c(2,8))
rownames(ou.com.table)<-c("HL","BW")
colnames(ou.com.table)<-
c("","","Lobs","Lconstrained","LRTest",
"Prob","AICc.obs","AICc.constrained")
ou.com.table[1:2,1:2]<-round(result.ou$Rconstrained,3)
ou.com.table[1,"Lconstrained"]<-round(result.ou$Lconstrained,3)
ou.com.table[1,"LRTest"]<-round(result.ou$LRTest,3)
ou.com.table[1,"Prob"]<-round(result.ou$Prob,3)
ou.com.table[1,"AICc.obs"]<-round(result.ou$AICc.obs,3)
ou.com.table[1,"AICc.constrained"]<-round(result.ou$AICc.constrained,3)
ou.com.table
## Lobs Lconstrained LRTest Prob AICc.obs AICc.constrained
## HL 1.024 0.158 NA -185.949 119.111 0 264.786 381.898
## BW 0.158 1.024 NA NA NA NA NA NA
xtable(rbind(ou.obs.table,ou.com.table))
## % latex table generated in R 4.0.2 by xtable 1.8-4 package
## % Wed Dec 23 22:40:08 2020
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrrrrr}
## \hline
## & V1 & V2 & Lobs & Lconstrained & LRTest & Prob & AICc.obs & AICc.constrained \\
## \hline
## HL & 0.23 & 0.13 & -126.39 & & & & & \\
## BW & 0.13 & 0.13 & & & & & & \\
## HL.1 & 1.02 & 0.16 & & -185.95 & 119.11 & 0.00 & 264.79 & 381.90 \\
## BW.1 & 0.16 & 1.02 & & & & & & \\
## \hline
## \end{tabular}
## \end{table}
####### EB ########
##OBS
eb.obs.table<-array(NA,c(2,8))
rownames(eb.obs.table)<-c("HL","BW")
colnames(eb.obs.table)<-
c("","","Lobs","Lconstrained","LRTest",
"Prob","AICc.obs","AICc.constrained")
eb.obs.table[1:2,1:2]<-round(result.eb$Robs,3)
eb.obs.table[1,"Lobs"]<-round(result.eb$Lobs,3)
eb.obs.table
## Lobs Lconstrained LRTest Prob AICc.obs AICc.constrained
## HL 0.228 0.124 -126.861 NA NA NA NA NA
## BW 0.124 0.119 NA NA NA NA NA NA
## COMMON
eb.com.table<-array(NA,c(2,8))
rownames(eb.com.table)<-c("HL","BW")
colnames(eb.com.table)<-
c("","","Lobs","Lconstrained","LRTest",
"Prob","AICc.obs","AICc.constrained")
eb.com.table[1:2,1:2]<-round(result.eb$Rconstrained,3)
eb.com.table[1,"Lconstrained"]<-round(result.eb$Lconstrained,3)
eb.com.table[1,"LRTest"]<-round(result.eb$LRTest,3)
eb.com.table[1,"Prob"]<-round(result.eb$Prob,3)
eb.com.table[1,"AICc.obs"]<-round(result.eb$AICc.obs,3)
eb.com.table[1,"AICc.constrained"]<-round(result.eb$AICc.constrained,3)
eb.com.table
## Lobs Lconstrained LRTest Prob AICc.obs AICc.constrained
## HL 1.015 0.125 NA -188.316 122.909 0 265.723 386.632
## BW 0.125 1.015 NA NA NA NA NA NA
xtable(rbind(eb.obs.table,eb.com.table))
## % latex table generated in R 4.0.2 by xtable 1.8-4 package
## % Wed Dec 23 22:40:08 2020
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrrrrr}
## \hline
## & V1 & V2 & Lobs & Lconstrained & LRTest & Prob & AICc.obs & AICc.constrained \\
## \hline
## HL & 0.23 & 0.12 & -126.86 & & & & & \\
## BW & 0.12 & 0.12 & & & & & & \\
## HL.1 & 1.01 & 0.12 & & -188.32 & 122.91 & 0.00 & 265.72 & 386.63 \\
## BW.1 & 0.12 & 1.01 & & & & & & \\
## \hline
## \end{tabular}
## \end{table}