Rate OU Model

Published

December 27, 2022

rm(list=ls())

library(ape)
library(geiger)
library(nlme)
library(phytools)

Loading required package: maps

ouraterescov<-function(regparam,modparam,y0=y0,x0=x0,tree=tree){
  a = regparam["a"]
  b = regparam["b"]
  
  alp   = modparam["alp"]
  sigma = modparam["sigma"]
  theta = modparam["theta"]
  
  Cmtx <-vcv(tree)
  round(Cmtx,3)
  #############################################################################
  #############Var[y]##########################################################
  #############################################################################
  
  covyi2yj2mtx.ou<-array(NA,dim(Cmtx))
  for(i in 1:dim(Cmtx)[1]){
    #i<-1
    t <- Cmtx[i,i]
    for(j in 1:dim(Cmtx)[2]){
      #  j<-1
      ta<- Cmtx[i,j]
      
      Ey4.ou.ta=y0^4+6*a*(y0^2*ta+a*ta^2/2+b*x0/alp*(ta-1/alp*(1-exp(-alp*ta)))+b*theta*(ta^2/2-1/alp*(ta-(1-exp(-alp*ta))/alp)))+6*b*(y0^2*x0*ta+a*x0/alp*(ta-(1-exp(-alp*ta))/alp)+a*theta*(ta^2/2-1/alp*(ta-(1-exp(-alp*ta))/alp))+b*x0^2/2/alp*(ta-(1-exp(-2*alp*ta))/2/alp)+b*theta^2/2*ta^2-2*b*theta^2/alp*(ta-(1-exp(-alp*ta))/alp)+b*theta^2/2/alp*(ta-(1-exp(-2*alp*ta))/2/alp)+2*x0*theta*b/alp*(ta-(1-exp(-alp*ta))/alp)-2*x0*theta*b/2/alp*(ta-(1-exp(-2*alp*ta))/2/alp)+ b*sigma^2*(ta^2/2/2/alp- (ta-(1-exp(-2*alp*ta))/2/alp)/((2*alp)^2)))
      
      Ey2.ou.ta =y0^2+a*ta+b*x0*1/alp*(1-exp(-alp*ta))+b*theta*(ta-1/alp*(1-exp(-alp*ta)))
      
      Vy2.ou.ta = Ey4.ou.ta - Ey2.ou.ta^2
      
      Vx.ta = sigma^2*(1-exp(-2*alp*ta))/2/alp
      
      
      Ey2x.ou.ta = y0^2*x0 + a*x0*(1-exp(-alp*ta))/alp+a*theta*(ta-(1-exp(-alp*ta))/alp)+b*x0^2*(1-exp(-2*alp*ta))/(2*alp)+b*theta^2*ta-2*b*theta^2*(1-exp(-alp*ta))/alp+b*theta^2*(1-exp(-2*alp*ta))/(2*alp)+2*x0*theta*b*((1-exp(-alp*ta))/alp - (1-exp(-2*alp*ta))/(2*alp)) + b*sigma^2*(ta/(2*alp)-(1-exp(-2*alp*ta))/((2*alp)^2))
      
      Ex.ou.ta = x0*exp(-alp*ta)+theta*(1-theta^(1-exp(-alp*ta)))
      
      covy2x.ou = Ey2x.ou.ta - Ey2.ou.ta * Ex.ou.ta
      
      covyi2yj2mtx.ou[i,j]<-Vy2.ou.ta+b^2*(1-exp(-alp*(t-ta)))^2/alp/alp*Vx.ta+2*b*(1-exp(-alp*(t-ta)))/alp*covy2x.ou
    }
  }
  
  # covyi2yj2mtx.ou
  
  #############################################################################
  ############Cov[y2,x]########################################################
  #############################################################################
  
  covyi2ximtx.ou<-array(NA,dim(Cmtx))
  for(i in 1:dim(Cmtx)[1]){
    t <- Cmtx[i,i]
    for(j in 1:dim(Cmtx)[2]){
      ta<- Cmtx[i,j]
      
      Eya2xa.ou = y0^2*x0+a*x0*(1-exp(-alp*ta))/alp+a*theta*(ta-(1-exp(-alp*ta))/alp)+b*x0^2*((1-exp(-2*alp*ta))/2/alp)+b*theta^2*ta-2*b*theta^2*(1-exp(-alp*ta))/alp+b*theta^2*(1-exp(-2*alp*ta))/2/alp+2*x0*theta*b*((1-exp(-alp*ta))/alp-(1-exp(-2*alp*ta))/2/alp)+b*sigma^2*(ta/2/alp-(1-exp(-2*alp*ta))/((2*alp)^2))
      Eya2.ou = y0^2+a*ta+b*x0*((1-exp(-alp*ta))/alp)+b*theta*(ta-(1-exp(-alp*ta))/alp)
      Exa.ou = x0*exp(-alp*ta)+theta*(1-exp(-alp*ta))
      
      covya2xa.ou = Eya2xa.ou - Eya2.ou*Exa.ou
      covyi2xi.ou = exp(-alp*(t-ta))*covya2xa.ou +  b*(1-exp(-alp*(t-ta)))/alp*exp(-alp*(t-ta))* sigma^2*(1-exp(-2*alp*ta))/2/alp
      
      
      covyi2ximtx.ou[i,j]<-exp(-alp*(t-ta))*covya2xa.ou +  b*(1-exp(-alp*(t-ta)))/alp*exp(-alp*(t-ta))* sigma^2*(1-exp(-2*alp*ta))/2/alp
    }
  }
  
  # covyi2ximtx.ou
  
  
  #############################################################################
  #############Var[x]##########################################################
  #############################################################################
  
  covxixjmtx.ou<-array(NA,dim(Cmtx))
  for(i in 1:dim(Cmtx)[1]){
    t <- Cmtx[i,i]
    for(j in 1:dim(Cmtx)[2]){
      ta<- Cmtx[i,j]
      covxixjmtx.ou[i,j]<-sigma^2*exp(-2*alp*(t-ta))*(1-exp(-2*alp*ta))/2/alp
    }
  }
  
  # covxixjmtx.ou
  
  #############################################################################
  ######################Var[r]#################################################
  #############################################################################
  
  covresmtx.ou<-covyi2yj2mtx.ou-covyi2ximtx.ou%*%solve(covxixjmtx.ou)%*%t(covyi2ximtx.ou)
  
  return(covresmtx.ou)
}

regcoef.ou<-function(regparam,modparam,tree=tree,trait.x=trait.x,trait.y=trait.y,x0=x0,y0=y0){
  Dmtx<- cbind(rep(1,length(trait.x)),trait.x)
  Vr.ou<-ouraterescov(regparam,modparam,y0=y0,x0=x0,tree=tree)
  var.beta.hat.ou<-solve(t(Dmtx)%*%solve(Vr.ou)%*%Dmtx)
  beta.hat.ou <- var.beta.hat.ou%*%t(Dmtx)%*%solve(Vr.ou)%*%trait.y
  return(beta.hat.ou)
}




library(ape)
set.seed(2312)
n<-30
tree<-rcoal(n)
plot(tree)

trait.x<-matrix(rnorm(n),ncol = 1)
rownames(trait.x)<-tree$tip.label
trait.x

           [,1]
t23  0.72089130
t11 -1.36576152
t6  -0.27172013
t26  0.75459575
t15  0.94564537
t14 -2.49680212
t4  -1.67544359
t5  -0.15740153
t9  -0.28355979
t24 -0.11475313
t25  0.20864819
t12 -1.32398540
t1   0.36320511
t18 -1.50188435
t7   0.19872009
t3   0.77535322
t20 -0.91582873
t27  0.32208681
t16 -0.71889675
t21 -0.18015587
t2   0.14392351
t17 -0.33950272
t10 -0.98996927
t29  0.95623534
t22  1.38585891
t28 -1.28059886
t19  0.25374240
t30 -0.09573423
t8   0.80210522
t13 -0.23200392

#trait.y<-matrix(rnorm(n),ncol = 1)
trait.y<-3+5*trait.x + rnorm(n)
rownames(trait.y)<-tree$tip.label
trait.y

           [,1]
t23   6.4760109
t11  -4.6712806
t6   -0.4752605
t26   6.6063276
t15   8.6659486
t14 -10.8822390
t4   -4.7117791
t5    2.9891925
t9    1.7824433
t24   4.1215504
t25   4.1538349
t12  -2.9398159
t1    5.4245610
t18  -1.3284929
t7    3.9857832
t3    7.2440000
t20  -0.9640329
t27   4.3141687
t16  -1.0787281
t21   3.9318172
t2    3.4194394
t17   0.6585404
t10  -0.7839282
t29   7.5185322
t22  10.0384923
t28  -5.0653358
t19   3.2882584
t30   2.2639685
t8    6.2387700
t13   0.7484814

## use initial ols estimate
lm.ols<-lm(trait.y~trait.x)
regparam<-lm.ols$coefficients
names(regparam)<-c("a","b")
regparam

       a        b 
3.045972 4.974641

regparam.ols<-regparam
regparam.ini<-regparam

## Fit ou model for both trait and get root value 
fitou.x<-fitContinuous(tree, dat=trait.x, model="OU")

Warning in fitContinuous(tree, dat = trait.x, model = "OU"): 
Parameter estimates appear at bounds:
    alpha

fitbm.y<-fitContinuous(tree, dat=trait.y, model="BM")
y0 = fitbm.y$opt$z0 #1
x0 = fitou.x$opt$z0 

## get the initital estimate for sigma
sigma <- 1#sqrt(fitou.x$opt$sigsq)
alp   =  0.1#fitou.x$opt$alpha# 1e-5 #fitou.y$opt$alpha#1e-5
theta = 0#fitou.x$opt$z0  #0.8 #theta cannot be negative because  Ex.ou.ta has term 
modparam<-c(alp,theta,sigma)
names(modparam)<-c("alp","theta","sigma")
modparam

  alp theta sigma 
  0.1   0.0   1.0

regcoef.ou(regparam,modparam,tree=tree,trait.x=trait.x,trait.y=trait.y,x0=x0,y0=y0)

         [,1]
[1,] 2.865300
[2,] 7.109762

## set up bound for optimization search
sum.lm.ols<-summary(lm.ols)
bdd.a=sum.lm.ols$coefficients[,"Estimate"]-5*sum.lm.ols$coefficients[,"Std. Error"]
bdd.b=sum.lm.ols$coefficients[,"Estimate"]+5*sum.lm.ols$coefficients[,"Std. Error"]
bdd.a

(Intercept)     trait.x 
   2.005525    3.851830

bdd.b

(Intercept)     trait.x 
   4.086419    6.097453

### try different initial points
# testData<-cbind(trait.x,trait.y)
# colnames(testData)<-c("trait.x","trait.y")
# testData<-as.data.frame(testData)
# lm.bm<-gls(trait.y ~ trait.x, correlation = corBrownian(phy = tree),data = testData, method = "ML")
# lm.bm$coefficients
# lm.ou<-gls(trait.y ~ trait.x, correlation = corMartins(1, phy = tree),data = testData, method = "ML")
# lm.ou$coefficients
# lm.ols
###


cutci=10
CI.a=seq(bdd.a[1],bdd.a[2],length=cutci)
CI.b=seq(bdd.b[1],bdd.b[2],length=cutci)


Dmtx<-cbind(rep(1,n),trait.x)
beta.hat.ou.array<-array(NA,c(length(CI.a),length(CI.b),2))
dim(beta.hat.ou.array)

[1] 10 10  2

for(aIndex in 1:length(CI.a)){
  #  aIndex<-1
  a<-CI.a[aIndex]
  for(bIndex in 1:length(CI.b)){
    #    bIndex<-1
    b<-CI.b[bIndex]
    regparam<-c(a,b)    
    names(regparam)<-c("a","b")
    Vr.ou<-ouraterescov(regparam,modparam,y0=y0,x0=x0,tree=tree)
    var.beta.hat.ou<-solve(t(Dmtx)%*%solve(Vr.ou)%*%Dmtx)
    beta.hat.ou <- var.beta.hat.ou%*%t(Dmtx)%*%solve(Vr.ou)%*%trait.y
    beta.hat.ou.array[aIndex,bIndex,]<-beta.hat.ou
  }
}


for(aIndex in 1:length(CI.a)){
  for(bIndex in 1:length(CI.b)){
    print(beta.hat.ou.array[aIndex,bIndex,])
  }
}

[1] 2.866538 7.109417
[1] 2.867602 7.109128
[1] 2.868629 7.108862
[1] 2.869614 7.108618
[1] 2.870555 7.108394
[1] 2.871453 7.108188
[1] 2.872305 7.108001
[1] 2.873112 7.107829
[1] 2.873875 7.107671
[1] 2.874594 7.107527
[1] 2.865360 7.109752
[1] 2.866385 7.109459
[1] 2.867380 7.109187
[1] 2.868343 7.108935
[1] 2.869270 7.108703
[1] 2.870159 7.108488
[1] 2.87101 7.10829
[1] 2.871822 7.108108
[1] 2.872595 7.107940
[1] 2.873328 7.107785
[1] 2.864342 7.110054
[1] 2.865325 7.109760
[1] 2.866285 7.109486
[1] 2.867220 7.109229
[1] 2.868127 7.108991
[1] 2.869002 7.108769
[1] 2.869845 7.108564
[1] 2.870654 7.108373
[1] 2.871429 7.108196
[1] 2.872169 7.108033
[1] 2.863458 7.110326
[1] 2.864397 7.110034
[1] 2.865321 7.109760
[1] 2.866225 7.109501
[1] 2.867107 7.109259
[1] 2.867964 7.109033
[1] 2.868793 7.108822
[1] 2.869594 7.108625
[1] 2.870364 7.108441
[1] 2.871105 7.108271
[1] 2.862686 7.110570
[1] 2.863582 7.110283
[1] 2.864468 7.110010
[1] 2.865340 7.109752
[1] 2.866195 7.109509
[1] 2.867029 7.109280
[1] 2.867842 7.109065
[1] 2.868629 7.108863
[1] 2.869392 7.108674
[1] 2.870127 7.108498
[1] 2.862009 7.110790
[1] 2.862863 7.110508
[1] 2.863711 7.110239
[1] 2.864550 7.109983
[1] 2.865376 7.109740
[1] 2.866187 7.109510
[1] 2.866979 7.109293
[1] 2.867751 7.109088
[1] 2.868501 7.108896
[1] 2.869229 7.108715
[1] 2.861412 7.110987
[1] 2.862226 7.110712
[1] 2.863037 7.110448
[1] 2.863843 7.110195
[1] 2.864640 7.109954
[1] 2.865424 7.109724
[1] 2.866195 7.109506
[1] 2.866949 7.109300
[1] 2.867685 7.109105
[1] 2.868401 7.108921
[1] 2.860885 7.111163
[1] 2.861660 7.110896
[1] 2.862435 7.110638
[1] 2.863208 7.110390
[1] 2.863975 7.110151
[1] 2.864734 7.109923
[1] 2.865481 7.109706
[1] 2.866216 7.109499
[1] 2.866936 7.109303
[1] 2.867639 7.109117
[1] 2.860418 7.111321
[1] 2.861155 7.111063
[1] 2.861895 7.110811
[1] 2.862636 7.110568
[1] 2.863374 7.110334
[1] 2.864107 7.110108
[1] 2.864831 7.109893
[1] 2.865545 7.109686
[1] 2.866247 7.109490
[1] 2.866935 7.109302
[1] 2.860002 7.111463
[1] 2.860704 7.111213
[1] 2.861411 7.110969
[1] 2.862121 7.110732
[1] 2.862830 7.110502
[1] 2.863536 7.110280
[1] 2.864237 7.110066
[1] 2.864930 7.109862
[1] 2.865613 7.109665
[1] 2.866285 7.109478

################################
######Random select point#######
################################
## not convergence easily
regparam.ini

       a        b 
3.045972 4.974641

regparam.ols

       a        b 
3.045972 4.974641

eps<-1e-5
while(TRUE){
  a=runif(1,bdd.a[1],bdd.a[2])
  b=runif(1,bdd.b[1],bdd.b[2])
  regparam<-c(a,b)    
  names(regparam)<-c("a","b")
  regparam
  Vr.ou<-ouraterescov(regparam,modparam,y0=y0,x0=x0,tree=tree)
  var.beta.hat.ou<-solve(t(Dmtx)%*%solve(Vr.ou)%*%Dmtx)
  beta.hat.ou <- var.beta.hat.ou%*%t(Dmtx)%*%solve(Vr.ou)%*%trait.y
  beta.hat.ou
  if(sqrt(sum((beta.hat.ou-regparam.ini)^2))<eps){
    break
  }else{
    regparam.ini<-beta.hat.ou
  }
#  print(beta.hat.ou)
}
names(beta.hat.ou)<-c("a","b")
beta.hat.ou

         [,1]
[1,] 2.866654
[2,] 7.109385
attr(,"names")
[1] "a" "b"

################################################################