library(ape)
library(MASS)
library(phytools)
## Loading required package: maps
library(geiger)
## Registered S3 method overwritten by 'geiger':
## method from
## unique.multiPhylo ape
CompareRates.multTrait<-function(phy,x,TraitCov=T,ms.err=NULL,ms.cov=NULL){
build.chol<-function(b){
c.mat<-matrix(0,nrow = p,ncol = p)
c.mat[lower.tri(c.mat)]<-b[-1]
c.mat[p,p]<-exp(b[1])
c.mat[1,1]<-sqrt(sum((c.mat[p,])^2))
if(p>2){
for(i in 2:(p-1)){
c.mat[i,i]<-ifelse((c.mat[1,1]^2-sum((c.mat[i,])^2))>0,sqrt(c.mat[1,1]^2-sum((c.mat[i,])^2)),0)
}
}
return(c.mat)
}
x<-as.matrix(x)
N<-nrow(x)
p<-ncol(x)
C<-vcv.phylo(phy)
C<-C[rownames(x),rownames(x)]
I<-diag(1,N)
if (is.matrix(ms.err)){
ms.err<-as.matrix(ms.err[rownames(x),])}
if (is.matrix(ms.cov)){
ms.cov<-as.matrix(ms.cov[rownames(x),])}
# a.obs<-colSums(solve(C)) %*% x / sum(solve(C))
# R.obs<-t(x-one%*%a.obs)%*%solve(C)%*%(x-one%*%a.obs)/N
IIDcon<-function(trait=trait,mu=mu,sigma=sigma,C=C){
z<- (trait-mu)/sigma
eC<-eigen(C)
D.n5<-diag(1/sqrt(eC$values))
C.neg.5<-eC$vectors%*%D.n5%*%t(eC$vectors)
return(C.neg.5%*%trait)
}
one<-matrix(1,N,1)
a.obs.eb<-array(NA,p)
R.obs.eb<-array(NA,c(p,p))
IIDtrait<-array(NA,dim(x))
r.array<-array(NA,p)
for(Index in 1:p){
# ?phylosig
# ?fitContinuous
# Index<-1
r<-fitContinuous(phy=phy,dat=x[,Index],model="EB")$opt$a
r.array[Index]<-r
### HERE FOR EB
Sr<- (exp(r*C)-1)/r
assign(paste("Sr",Index,sep=""),Sr)
a.obs.eb[Index]<-colSums(solve(Sr)) %*% x[,Index] / sum(solve(Sr))
R.obs.eb[Index,Index]<-t(x[,Index]-one%*%a.obs.eb[Index])%*%solve(Sr)%*%(x[,Index]-one%*%a.obs.eb[Index])/N
IIDtrait[,Index]<-IIDcon(trait=x[,Index],mu=a.obs.eb[Index],sigma=R.obs.eb[Index,Index],C=Sr)
}
for(i in 1: (p-1)){
for( j in (i+1):p){
R.obs.eb[i,j]<-R.obs.eb[j,i]<- cov(IIDtrait[,i],IIDtrait[,j])
}
}
dim(a.obs.eb)
a.obs.eb<-matrix(a.obs.eb,nrow=1,ncol=p)
a.obs.eb
dim(a.obs.eb)
R.obs.eb
D<-matrix(0,N*p,p)
for(i in 1:(N*p)){
for(j in 1:p){
if((j-1)*N < i && i<=j*N){
D[i,j]=1.0
}
}
}
y<-as.matrix(as.vector(x))
if (TraitCov==F){R.obs.eb<-diag(diag(R.obs.eb),p)}
Sr1
Sr2
r1<-r.array[1]
r2<-r.array[2]
#alpha12<-(alpha1+alpha2)/2
#alpha12<-sqrt(alpha1*alpha2)
#alpha12<-1/(1/alpha1+1/alpha2)
##################USE THIS ONE THROUOUT##############################
Sr12 <- (exp((r1+r2)/2*C)-1)/ (r1/2+r2/2)
#####################################################################
RkronSr<-rbind(cbind(R.obs.eb[1,1]*Sr1,R.obs.eb[1,2]*Sr12),
cbind(R.obs.eb[1,2]*Sr12,R.obs.eb[2,2]*Sr2))
RkronSr
LLik.obs.eb<-ifelse(is.matrix(ms.err)==TRUE,
-t(y-D%*%t(a.obs.eb))%*%ginv((RkronSr+diag(as.vector(ms.err))))%*%(y-D%*%t(a.obs.eb))/2-N*p*log(2*pi)/2-determinant((RkronSr+ diag(as.vector(ms.err))))$modulus[1]/2,
-t(y-D%*%t(a.obs.eb))%*%ginv(RkronSr)%*%(y-D%*%t(a.obs.eb))/2-N*p*log(2*pi)/2-determinant(RkronSr)$modulus[1]/2)
# sigma.mn<-mean(diag(R.obs))
sigma.mn.eb<-mean(diag(R.obs.eb))
if(is.matrix(ms.err) && is.matrix(ms.cov)){
within.spp<-cbind(ms.err,ms.cov)
rc.label<-NULL
for(i in 1:p){
rc.label<-rbind(rc.label,c(i,i))
}
for(j in 2:p){
if(i!=j&&i<j){
rc.label<-rbind(rc.label,c(i,j))
}
}
m.e<-NULL
for(i in 1:p){
temp<-NULL
for(j in 1:p){
for(k in 1:nrow(rc.label)){
if(setequal(c(i,j),rc.label[k,])==T)
{tmp<-cbind(tmp,diag(within.spp[,k]))}
}
}
m.e<-rbind(m.e,tmp)
}
}
R<-R.obs.eb
#param<-c(sigma.mn.eb,alpha1,alpha2)
#names(param)<-c("sigma","alpha1","alpha2")
lik.covF.eb<-function(param){
alpha1<-param["alpha1"]
alpha2<-param["alpha2"]
sigma<-param["sigma"]
alpha12<-(alpha1+alpha2)/2
##########################################
Sr12 <- (exp((r1+r2)/2*C)-1)/ (r1/2+r2/2)
#########################################
diag(R)<-sigma.mn.eb
RkronSr<-rbind(cbind(R[1,1]*Sr1,R[1,2]*Sr12),
cbind(R[2,1]*Sr12,R[2,2]*Sr2))
LLik<-ifelse(is.matrix(ms.err)==TRUE,
-t(y-D%*%t(a.obs.eb))%*%ginv((RkronSr+m.e))%*%(y-D%*%t(a.obs.eb))/2-N*p*log(2*pi)/2-determinant((RkronSr+ m.e))$modulus[1]/2,
-t(y-D%*%t(a.obs.eb))%*%ginv(RkronSr)%*%(y-D%*%t(a.obs.eb))/2-N*p*log(2*pi)/2-determinant(RkronSr)$modulus[1]/2)
if(LLik==-Inf){LLik<--1e+10}
return(-LLik)
}
#param<-c(alpha1,alpha2,sigma.mn.eb,0)
#names(param)<-c("alpha1","alpha2","sigma","R.offd")
lik.covT.eb<-function(param){
alpha1<-param["alpha1"]
alpha2<-param["alpha2"]
sigma<-param["sigma"]
R.offd<-param["R.offd"]
alpha12<-(alpha1+alpha2)/2
#########################################
Sr12 <- (exp((r1+r2)/2*C)-1)/ (r1/2+r2/2)
#########################################
low.chol<-build.chol(c(sigma,R.offd))
R<-low.chol%*%t(low.chol)
RkronSr<-rbind(cbind(R[1,1]*Sr1,R[1,2]*Sr12),
cbind(R[1,2]*Sr12,R[2,2]*Sr2))
LLik<-ifelse(is.matrix(ms.err)==TRUE,
-t(y-D%*%t(a.obs.eb))%*%ginv((RkronSr+m.e))%*%(y-D%*%t(a.obs.eb))/2-N*p*log(2*pi)/2-determinant((RkronSr+ m.e))$modulus[1]/2,
-t(y-D%*%t(a.obs.eb))%*%ginv(RkronSr)%*%(y-D%*%t(a.obs.eb))/2-N*p*log(2*pi)/2-determinant(RkronSr)$modulus[1]/2)
# print(LLik)
if(LLik==-Inf){LLik<--1e+10}
return(-LLik)
}
sigma.upper<-2*max(apply(x,2,sd))
p0<-c(r1,r2,sigma.mn.eb)
names(p0)<-c("r1","r2","sigma")
if(TraitCov==F){model.eb<-
optim(p0,fn=lik.covF.eb,method="L-BFGS-B",lower = c(1e-5,1e-5,1e-5),upper=c(5,5,sigma.upper))
}
R.offd<-rep(0,(p*(p-1)/2))
p0<-c(r1,r2,sigma.mn.eb,0)
names(p0)<-c("r1","r2","sigma","R.offd")
if(TraitCov==T){model1.eb<-
optim(p0,fn=lik.covT.eb,method="L-BFGS-B",lower = c(1e-5,1e-5,1e-5,1e-5),upper=c(5,5,sigma.upper,sigma.upper))
}
if(TraitCov==F){R.constr.eb<-diag(model.eb$par["sigma"],p)}
if(TraitCov==T){
chol.mat<-build.chol(model1.eb$par[c("sigma","R.offd")])
R.constr.eb<-chol.mat%*%t(chol.mat)
}
if(model1.eb$convergence==0){
message.eb<-"Optimization has converged."}else{
message.eb<-"Optim may not have converrged.
Consideer changing startt value or lower/upper limits."}
LRT.eb<-(-2*((-model1.eb$value-LLik.obs.eb)))
LRT.prob.eb<-pchisq(LRT.eb, (p-1),lower.tail = FALSE)
AIC.obs.eb<- -2*LLik.obs.eb+2*p+2*p+2*p #(2p twice: 1x for rates, 1x for anc.states,1x for h)
AIC.common.eb<--2*(-model1.eb$value)+2+2*p+2*p #(2*1:for 1 rate 2p for anc.states)
return(
list(
Robs.eb=R.obs.eb,
Rconstrained.eb=R.constr.eb,
Lobs.eb=LLik.obs.eb,
Lconstrained.eb=(-model1.eb$value),
LRTest.eb=LRT.eb,
Prob.eb=LRT.prob.eb,
AICc.obs.eb=AIC.obs.eb,
AICc.constrained.eb=AIC.common.eb,
optimmessage.eb=message.eb
)
)
}
### Sample code
phy<-rcoal(5)
plot(phy)
x<- matrix(c(rnorm(5,2,1),rnorm(5,0,0.5)),ncol=2)
rownames(x)<-phy$tip.label#LETTERS[1:N]
CompareRates.multTrait(phy=phy,x=x,TraitCov=T,ms.err=NULL,ms.cov=NULL)
## $Robs.eb
## [,1] [,2]
## [1,] 63.67079 80.77437
## [2,] 80.77437 65.89467
##
## $Rconstrained.eb
## [,1] [,2]
## [1,] 20.130945 6.489075
## [2,] 6.489075 20.130945
##
## $Lobs.eb
## [1] -26.74995
##
## $Lconstrained.eb
## [1] -32.08474
##
## $LRTest.eb
## [1] 10.66958
##
## $Prob.eb
## [1] 0.001089117
##
## $AICc.obs.eb
## [1] 65.4999
##
## $AICc.constrained.eb
## [1] 74.16949
##
## $optimmessage.eb
## [1] "Optimization has converged."
#####
tree<-read.tree("http://tonyjhwueng.info/phyrates/ple.nwk")
plot(tree)
tree$tip.label
## [1] "Plethodon_dorsalis" "Plethodon_ventralis"
## [3] "Plethodon_angusticlavius" "Plethodon_welleri"
## [5] "Plethodon_punctatus" "Plethodon_wehrlei"
## [7] "Plethodon_websteri" "Plethodon_teyahalee"
## [9] "Plethodon_cylindraceus" "Plethodon_variolatus"
## [11] "Plethodon_chlorobryonis" "Plethodon_chattahoochee"
## [13] "Plethodon_cheoah" "Plethodon_shermani"
## [15] "Plethodon_amplus" "Plethodon_meridianus"
## [17] "Plethodon_montanus" "Plethodon_albagula"
## [19] "Plethodon_sequoyah" "Plethodon_ocmulgee"
## [21] "Plethodon_savannah" "Plethodon_grobmani"
## [23] "Plethodon_kisatchie" "Plethodon_mississippi"
## [25] "Plethodon_kiamichi" "Plethodon_aureolus"
## [27] "Plethodon_glutinosus" "Plethodon_jordani"
## [29] "Plethodon_metcalfi" "Plethodon_ouachitae"
## [31] "Plethodon_fourchensis" "Plethodon_caddoensis"
## [33] "Plethodon_kentucki" "Plethodon_petraeus"
## [35] "Plethodon_yonahlossee" "Plethodon_hubrichti"
## [37] "Plethodon_nettingi" "Plethodon_richmondi"
## [39] "Plethodon_electromorphus" "Plethodon_cinereus"
## [41] "Plethodon_shenandoah" "Plethodon_hoffmani"
## [43] "Plethodon_virginia" "Plethodon_serratus"
df<-read.csv("http://tonyjhwueng.info/phyrates/Adams2012-SystBiolData.csv")
spX<-strsplit(as.character(df$X),"_")
phy<-tree
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)
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
## $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."