R/mcmc.popsize.R

   1 ## mcmc.popsize.R (2010-11-16)
   2
   3 ##   Run reversible jump MCMC to sample demographic histories
   4
   5 ## Copyright 2004-2010 Rainer Opgen-Rhein and Korbinian Strimmer
   6
   7 ## Portions of this function are adapted from rjMCMC code by
   8 ## Karl W Broman (see http://www.biostat.wisc.edu/~kbroman/)
   9
  10 ## This file is part of the R-package `ape'.
  11 ## See the file ../COPYING for licensing issues.
  12
  13 # public function
  14
  15 # run rjMCMC chain
  16 mcmc.popsize <-
  17   function(tree, nstep,
  18     thinning=1, burn.in=0, progress.bar=TRUE,
  19     method.prior.changepoints=c("hierarchical", "fixed.lambda"),
  20     max.nodes=30,
  21     lambda=0.5,    # "fixed.lambda" method.prior.changepoints
  22     gamma.shape=0.5, gamma.scale=2,  # gamma distribution from which lambda is drawn (for "hierarchical" method)
  23     method.prior.heights=c("skyline", "constant", "custom"),
  24     prior.height.mean,
  25     prior.height.var
  26     )
  27 {
  28   method.prior.changepoints <- match.arg(method.prior.changepoints)
  29   method.prior.heights <- match.arg(method.prior.heights)
  30
  31
  32   #Calculate skylineplot, coalescent intervals and estimated population sizes
  33
  34   if(attr(tree, "class")=="phylo")
  35     {ci <- coalescent.intervals(tree)
  36      sk1 <- skyline(ci)}
  37   else if (attr(tree, "class")=="coalescentIntervals")
  38     {ci<-tree
  39     sk1<-skyline(ci)}
  40   else
  41   stop("tree must be an object of class phylo or coalescentIntervals")
  42
  43    #consider possibility of more than one lineage
  44    ci$lineages<-ci$lineages[sk1$interval.length>0]
  45    ci$interval.length<-ci$interval.length[sk1$interval.length>0]
  46    data<-sk1$time<-sk1$time[sk1$interval.length>0]
  47    sk1$population.size<-sk1$population.size[sk1$interval.length>0]
  48    sk1$interval.length<-sk1$interval.length[sk1$interval.length>0]
  49
  50   # constant prior for heights
  51
  52   if (method.prior.heights=="constant"){
  53     prior.height.mean<-function(position){
  54       return(mean(sk1$population.size))
  55     }
  56     prior.height.var<-function(position){
  57       return((mean(sk1$population.size))^2)
  58    }
  59   }
  60
  61   # skyline plot prior for heights
  62
  63   if (method.prior.heights=="skyline"){
  64
  65     TIME<-sk1$time
  66     numb.interv<-10
  67     prior.change.times<-abs((0:numb.interv)*max(TIME)/numb.interv)
  68     prior.height.mean.all<-prior.height.var.all<-vector(length=numb.interv)
  69     for(p.int in 1:(numb.interv))
  70     {
  71       left<-p.int
  72       right<-p.int+1
  73       sample.pop<-sk1$population.size[sk1$time>=prior.change.times[left]&sk1$time<=prior.change.times[right]]
  74       while(length(sample.pop)<10){
  75         if(left>1){left<-left-1}
  76         if(right<length(prior.change.times)){right<-right+1}
  77         sample.pop<-sk1$population.size[sk1$time>=prior.change.times[left]&sk1$time<=prior.change.times[right]]
  78       }
  79       prior.height.mean.all[p.int]<-sum(sample.pop)/length(sample.pop)
  80       prior.height.var.all[p.int]<-sum((sample.pop-prior.height.mean.all[p.int])^2)/(length(sample.pop)-1)
  81     }
  82
  83     prior.height.mean<-function(position)
  84     {
  85       j<-sum(prior.change.times<=position)
  86       if(j>=length(prior.height.mean.all)){j<-length(prior.height.mean.all)}
  87       prior.mean<-prior.height.mean.all[j]
  88       prior.mean
  89     }
  90
  91     prior.height.var<-function(position)
  92     {
  93       j<-sum(prior.change.times<=position)
  94       if(j>=length(prior.height.var.all)){j<-length(prior.height.var.all)}
  95       prior.var<-prior.height.var.all[j]
  96       prior.var
  97     }
  98   }
  99
 100   if(method.prior.heights=="custom"){
 101     if(missing(prior.height.mean)||missing(prior.height.var)){
 102          stop("custom priors not specified")}
 103   }
 104
 105   #set prior
 106   prior<-vector(length=4)
 107   prior[4]<-max.nodes
 108
 109   # set initial position of markov chain and likelihood
 110   pos<-c(0,max(data))
 111   h<-c(rep(mean(sk1$population.size), 2))
 112
 113   b.lin<-choose(ci$lineages, 2)
 114   loglik<<-loglik.pop
 115
 116   #set lists for data
 117   count.it<-floor((nstep-burn.in)/thinning)
 118   save.pos <- save.h <- vector("list",count.it)
 119   save.loglik <- 1:count.it
 120   save.steptype <- 1:count.it
 121   save.accept <- 1:count.it
 122
 123   # calculate jump probabilities for given lambda of the prior
 124   if(method.prior.changepoints=="fixed.lambda")
 125   {
 126     prior[1]<-lambda
 127     jump.prob <- matrix(ncol=4,nrow=prior[4]+1)
 128     p <- dpois(0:prior[4],prior[1])/ppois(prior[4]+1,prior[1])
 129     bk <- c(p[-1]/p[-length(p)],0)
 130     bk[bk > 1] <- 1
 131     dk <- c(0,p[-length(p)]/p[-1])
 132     dk[dk > 1] <- 1
 133     mx <- max(bk+dk)
 134     bk <- bk/mx*0.9
 135     dk <- dk/mx*0.9
 136     bk[is.na(bk)]<-0     # added
 137     dk[is.na(dk)]<-0     # added
 138     jump.prob[,3] <- bk
 139     jump.prob[,4] <- dk
 140     jump.prob[1,2] <- 0
 141     jump.prob[1,1] <- 1-bk[1]-dk[1]
 142     jump.prob[-1,1] <- jump.prob[-1,2] <-
 143     (1-jump.prob[-1,3]-jump.prob[-1,4])/2
 144   }
 145
 146
 147   # calculate starting loglik
 148   curloglik <- loglik(data,pos,h,b.lin,sk1,ci)
 149
 150   count.i<-1
 151
 152   #set progress bar
 153   if(progress.bar==TRUE)
 154   {
 155     X11(width=3, height=0.7)
 156     par(mar=c(0.5,0.5,2,0.5))
 157     plot(x=c(0,0),y=c(0,1), type="l", xlim=c(0,1), ylim=c(0,1),
 158     main="rjMCMC in progress", ylab="", xlab="", xaxs="i", yaxs="i", xaxt="n", yaxt="n")
 159   }
 160
 161  #BEGIN CALCULATION
 162
 163   for(i in (1:nstep + 1))
 164   {
 165     #progress bar
 166     if(progress.bar==TRUE)
 167     {
 168       if(i %% 100 == 0)
 169       {
 170         z<-i/nstep
 171         zt<-(i-100)/(nstep)
 172         polygon(c(zt,zt,z,z), c(1,0,0,1), col="black")
 173       }
 174     }
 175
 176   # calculate jump probabilities without given lamda
 177   if(method.prior.changepoints=="hierarchical"){
 178     prior[1]<-rgamma(1,shape=gamma.shape,scale=gamma.scale)
 179     jump.prob <- matrix(ncol=4,nrow=prior[4]+1)
 180     p <- dpois(0:prior[4],prior[1])/ppois(prior[4]+1,prior[1])
 181     bk <- c(p[-1]/p[-length(p)],0)
 182     bk[bk > 1] <- 1
 183     dk <- c(0,p[-length(p)]/p[-1])
 184     dk[dk > 1] <- 1
 185     mx <- max(bk+dk)
 186     bk <- bk/mx*0.9
 187     dk <- dk/mx*0.9
 188     bk[is.na(bk)]<-0   # added
 189     dk[is.na(dk)]<-0   # added
 190     jump.prob[,3] <- bk
 191     jump.prob[,4] <- dk
 192     jump.prob[1,2] <- 0
 193     jump.prob[1,1] <- 1-bk[1]-dk[1]
 194     jump.prob[-1,1] <- jump.prob[-1,2] <-
 195     (1-jump.prob[-1,3]-jump.prob[-1,4])/2
 196   }
 197
 198     # determine what type of jump to make
 199     wh <- sample(1:4,1,prob=jump.prob[length(h)-1,])
 200
 201     if (i %% thinning == 0& i>burn.in) {save.steptype[[count.i]] <- wh}
 202
 203     if(wh==1) {
 204       step <- ht.move(data,pos,h,curloglik,prior, b.lin, sk1, ci, prior.height.mean, prior.height.var)
 205       h <- step[[1]]
 206       curloglik <- step[[2]]
 207       if(i%%thinning==0 & i>burn.in){
 208          save.pos[[count.i]]<-pos
 209          save.h[[count.i]]<-h
 210          save.loglik[[count.i]]<-step[[2]]
 211          save.accept[[count.i]]<-step[[3]]
 212          }
 213     }
 214     else if(wh==2) {
 215       step <- pos.move(data,pos,h,curloglik, b.lin,sk1,ci)
 216       pos <- step[[1]]
 217       curloglik <- step[[2]]
 218       if(i%%thinning==0 & i>burn.in){
 219           save.pos[[count.i]]<-pos
 220           save.h[[count.i]]<-h
 221           save.loglik[[count.i]]<-step[[2]]
 222           save.accept[[count.i]]<-step[[3]]
 223           }
 224     }
 225     else if(wh==3) {
 226       step <- birth.step(data,pos,h,curloglik,prior,jump.prob, b.lin, sk1, ci, prior.height.mean, prior.height.var)
 227       pos <- step[[1]]
 228       h <- step[[2]]
 229       curloglik <- step[[3]]
 230       if(i%%thinning==0 & i>burn.in){
 231          save.pos[[count.i]]<-pos
 232          save.h[[count.i]]<-h
 233          save.loglik[[count.i]]<-step[[3]]
 234          save.accept[[count.i]]<-step[[4]]
 235          }
 236     }
 237     else {
 238       step <- death.step(data,pos,h,curloglik,prior,jump.prob, b.lin, sk1, ci, prior.height.mean, prior.height.var)
 239       pos <- step[[1]]
 240       h <- step[[2]]
 241       curloglik <- step[[3]]
 242       if(i%%thinning==0 & i>burn.in){
 243          save.pos[[count.i]]<-pos
 244          save.h[[count.i]]<-h
 245          save.loglik[[count.i]]<-step[[3]]
 246          save.accept[[count.i]]<-step[[4]]
 247          }
 248     }
 249     if (i %% thinning == 0& i>burn.in) {count.i<-count.i+1}
 250   }
 251
 252   if(progress.bar==TRUE) dev.off()
 253
 254   list(pos=save.pos,h=save.h,loglik=save.loglik,
 255        steptype=save.steptype,accept=save.accept)
 256 }
 257
 258 # private functions
 259
 260 ht.move <-
 261 function(data,pos,h,curloglik,prior, b.lin,sk1,ci, prior.height.mean, prior.height.var)
 262 {
 263 #  print("ht.move")
 264   j <- sample(1:length(h),1)
 265
 266   prior.mean<-prior.height.mean(pos[j])
 267   prior.var<-prior.height.var(pos[j])
 268
 269   prior[3]<-prior.mean/prior.var
 270   prior[2]<-(prior.mean^2)/prior.var
 271
 272   newh <- h
 273   newh[j] <- h[j]*exp(runif(1,-0.5,0.5))
 274
 275   newloglik <- loglik(data,pos,newh,b.lin,sk1,ci)
 276   lr <- newloglik - curloglik
 277
 278   ratio <- exp(lr + prior[2]*(log(newh[j])-log(h[j])) - prior[3]*(newh[j]-h[j]))
 279
 280   if(runif(1,0,1) < ratio)
 281     return(list(newh,newloglik,1))
 282   else
 283     return(list(h,curloglik,0))
 284 }
 285
 286 pos.move <-
 287 function(data,pos,h,curloglik, b.lin,sk1,ci)
 288 {
 289 #  print("pos.move")
 290   if(length(pos)==3) j <- 2
 291   else j <- sample(2:(length(pos)-1),1)
 292   newpos <- pos
 293   left <- pos[j-1]
 294   right <- pos[j+1]
 295   newpos[j] <- runif(1,left,right)
 296
 297   newloglik <- loglik(data,newpos,h, b.lin,sk1,ci)
 298   lr <-  newloglik - curloglik
 299
 300   ratio <- exp(lr) * (right-newpos[j])*(newpos[j]-left)/
 301     (right-pos[j])/(pos[j]-left)
 302
 303   if(runif(1,0,1) < ratio)
 304     return(list(newpos,newloglik,1))
 305   else
 306     return(list(pos,curloglik,0))
 307 }
 308
 309 birth.step <-
 310 function(data,pos,h,curloglik,prior,jump.prob, b.lin, sk1, ci, prior.height.mean, prior.height.var)
 311 {
 312 #  print("birth")
 313   newpos <- runif(1,0,pos[length(pos)])
 314   j <- sum(pos < newpos)
 315
 316   left <- pos[j]
 317   right <- pos[j+1]
 318
 319   prior.mean<-prior.height.mean(pos[j])
 320   prior.var<-prior.height.var(pos[j])
 321   prior[3]<-prior.mean/prior.var
 322   prior[2]<-(prior.mean^2)/prior.var
 323
 324   u <- runif(1,-0.5,0.5)
 325   oldh<-(((newpos-left)/(right-left))*(h[j+1]-h[j])+h[j])
 326   newheight<-oldh*(1+u)
 327
 328   # ratio
 329   # recall that prior = (lambda, alpha, beta, maxk)
 330   k <- length(pos) - 2
 331   L <- max(pos)
 332
 333   prior.logratio <- log(prior[1]) - log(k+1) +  log((2*k+3)*(2*k+2)) - 2*log(L) +
 334     log(newpos-left) + log(right-newpos) - log(right-left) +
 335        prior[2]*log(prior[3]) - lgamma(prior[2]) +
 336         (prior[2]-1) * log(newheight) +
 337           prior[3]*(newheight)
 338
 339    proposal.ratio <- jump.prob[k+2,4]*L/jump.prob[k+1,3]/(k+1)
 340   jacobian <- (((newpos-left)/(right-left))*(h[j+1]-h[j]))+h[j]
 341
 342   # form new parameters
 343   newpos <- sort(c(pos,newpos))
 344   newh <- c(h[1:j], newheight, h[(j+1):length(h)])
 345
 346   newloglik <- loglik(data,newpos,newh, b.lin,sk1,ci)
 347
 348   lr <- newloglik - curloglik
 349
 350   ratio <- exp(lr + prior.logratio) * proposal.ratio * jacobian
 351
 352   if(runif(1,0,1) < ratio)
 353     return(list(newpos,newh,newloglik,1))
 354   else
 355     return(list(pos,h,curloglik,0))
 356 }
 357
 358 death.step <-
 359 function(data,pos,h,curloglik,prior,jump.prob, b.lin,sk1,ci, prior.height.mean, prior.height.var)
 360 {
 361 #  print("death")
 362   # position to drop
 363   if(length(pos)==3) j <- 2
 364   else j <- sample(2:(length(pos)-1),1)
 365
 366   left <- pos[j-1]
 367   right <- pos[j+1]
 368
 369   prior.mean<-prior.height.mean(pos[j])
 370   prior.var<-prior.height.var(pos[j])
 371   prior[3]<-prior.mean/prior.var
 372   prior[2]<-(prior.mean^2)/prior.var
 373
 374   # get new height
 375   h.left <- h[j-1]
 376   h.right <- h[j+1]
 377   newheight <- (((pos[j]-left)/(right-left))*(h.right-h.left)+h.left)
 378
 379   # ratio
 380   # recall that prior = (lambda, alpha, beta, maxk)
 381   k <- length(pos) - 3
 382   L <- max(pos)
 383
 384   prior.logratio <- log(k+1) - log(prior[1]) -  log(2*(k+1)*(2*k+3)) + 2*log(L) -
 385     log(pos[j]-left) - log(right-pos[j]) + log(right-left) -
 386       prior[2]*log(prior[3]) + lgamma(prior[2]) -
 387         (prior[2]-1) * log(newheight) -
 388           prior[3]*(newheight)
 389   proposal.ratio <- (k+1)*jump.prob[k+1,3]/jump.prob[k+2,4]/L
 390   jacobian <- ((pos[j]-left)/(right-left))*(h[j+1]-h[j-1])+h[j-1]
 391
 392   # form new parameters
 393   newpos <- pos[-j]
 394
 395   newh <- h[-j]
 396
 397   newloglik <- loglik(data,newpos,newh, b.lin,sk1,ci)
 398
 399   lr <- newloglik - curloglik
 400
 401   ratio <- exp(lr + prior.logratio) * proposal.ratio * (jacobian^(-1))
 402
 403   if(runif(1,0,1) < ratio)
 404     return(list(newpos,newh,newloglik,1))
 405   else
 406     return(list(pos,h,curloglik,0))
 407
 408
 409 }
 410
 411 # calculate the log likelihood for a set of data
 412 loglik.pop <-
 413 function(time=sk1$time, pos=c(0,max(sk1$time)), h=mean(sk1$population.size),b=b.lin,sk1,ci){
 414   data.time<-c(0,time)
 415
 416   leftside<-0
 417   i<-1
 418   h1<-c(h, h[length(h)])
 419   pos1<-c(pos, pos[length(pos)])
 420   while(i<length(time)){
 421     left.pos<-sum(data.time[i+1]>=pos)
 422     right.pos<-left.pos+1
 423 h.mix<-(((data.time[i+1]-pos[left.pos])/(pos[right.pos]-pos[left.pos]))*(h[right.pos]-h[left.pos]))+h[left.pos]
 424      leftside<-leftside+log(b[i]/h.mix)
 425     i<-i+1
 426   }
 427
 428   rightside<-0
 429   time1<-c(0,time)
 430   time.count<-1
 431
 432   # heigths of jumps
 433   jumps<-sort(c(time1, pos))
 434   h.jumps<-jumps
 435   while(time.count<=length(jumps)){
 436     left.pos<-sum(jumps[time.count]>=pos)
 437     right.pos<-left.pos+1
 438      h.jumps[time.count]<-(((jumps[time.count]-pos[left.pos])/(pos[right.pos]-pos[left.pos]))*(h[right.pos]-h[left.pos]))+h[left.pos]
 439     if(is.na(h.jumps[time.count])){h.jumps[time.count]<-h[left.pos]}
 440     time.count<-time.count+1
 441   }
 442
 443   # Vektor for lineages
 444   i<-1
 445   lineages.jumps<-jumps
 446    while(i<=length(jumps)){
 447      lineages.jumps[i]<-sum(jumps[i]>=time)
 448     if(lineages.jumps[i]==0){lineages.jumps[i]<-1}
 449      i<-i+1
 450    }
 451   lineage<-ci$lineages[lineages.jumps]
 452   b1<-choose(lineage, 2)
 453
 454   #Integral
 455   a<-(h.jumps[-1]-h.jumps[-length(h.jumps)])/(jumps[-1]-jumps[-length(jumps)])
 456   c<-h.jumps[-1]-jumps[-1]*a
 457   area<-(1/a)*log(a*jumps[-1]+c)-(1/a)*log(a*jumps[-length(jumps)]+c)
 458   stepfunction<-(jumps[-1]-jumps[-length(jumps)])/h.jumps[-1]
 459   area[is.na(area)]<-stepfunction[is.na(area)]
 460
 461   rightside<-sum(area*b1[-1])
 462
 463   loglik<-leftside-rightside
 464   loglik
 465 }