R/extract.popsize.R

In ape: Analyses of Phylogenetics and Evolution

## extract.popsize.R (2004-07-4)

##  Extract table with population size in dependence of time
##      from mcmc output generated by mcmc.popsize

## Copyright 2004 Rainer Opgen-Rhein and Korbinian Strimmer

## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.

extract.popsize<-function(mcmc.out, credible.interval=0.95, time.points=200, thinning=1, burn.in=0)
{

  # construct a matrix with the positions of the jumps
    b<-burn.in+1
    i<-1
    k<-array(dim=ceiling((length(mcmc.out$pos)-burn.in)/thinning))
    while(i<=length(k)) {
       k[i]<-length(mcmc.out$pos[[b]]);
       (i<-i+1);
       b<-b+thinning
    }
    o<-max(k)

    b<-burn.in+1
    i<-1
    pos.m<-matrix(nrow=length(k), ncol=o)
    while(i<=length(k)) {
        pos.m[i,]<-c(mcmc.out$pos[[b]], array(dim=o-length(mcmc.out$pos[[b]])));
        i<-i+1;
        b<-b+thinning
    }

  # construct a matrix with the heights of the jumps
    b<-burn.in+1
    i<-1
    h.m<-matrix(nrow=length(k), ncol=o)
    while(i<=length(k)) {
        h.m[i,]<-c(mcmc.out$h[[b]], array(dim=o-length(mcmc.out$h[[b]])));
        i<-i+1;
        b<-b+thinning
     }
  prep<-list("pos"=pos.m, "h"=h.m)

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

  step <- (max(prep$pos, na.rm=TRUE)-min(prep$pos, na.rm=TRUE))/(time.points-1)
  nr <- time.points

  p<-min(prep$pos, na.rm=TRUE)
  i<-1
  me<-matrix(nrow=nr, ncol=5)

  prep.l<-prep
  prep.l$pos<-cbind(prep$pos,prep$pos[,length(prep$pos[1,])])
  prep.l$h<-cbind(prep$h,prep$h[,length(prep$h[1,])])

  while (p<=max(prep$pos, na.rm=TRUE))
  {
    #Vector with position of heights
    l.prep<-prep$pos<=p
    l.prep[is.na(l.prep)]<-FALSE
    pos.of.h<-l.prep%*% array(data=1, dim=dim(prep$pos)[2])

    #Vector with heights
    z<-array(data=(1:dim(prep$pos)[1]), dim=dim(prep$pos)[1])
    index.left<-cbind(z,pos.of.h)
    index.right<-cbind(z, pos.of.h+1)

   mixed.heights<-((((p-prep$pos[index.left])/(prep$pos[index.right]-prep$pos[index.left]))*
                     (prep$h[index.right]-prep$h[index.left]))+prep$h[index.left])

    me[i,2]<-mean(mixed.heights)

    #library(MASS)
    #me[i,2]<-huber(mixed.heights)$mu

    me[i,3]<-median(mixed.heights)
    me[i,4]<-quantile(mixed.heights, probs=(1-credible.interval)/2, na.rm=TRUE)
    me[i,5]<-quantile(mixed.heights, probs=(1+credible.interval)/2, na.rm=TRUE)
    me[i,1]<-p
    p<-p+step
    i<-i+1
  }

  #av.jumps<-round((length(prep$pos)-sum(is.na(prep$pos)))/length(prep$pos[,1])-2,2)
  #print("average jumps")

  #print((length(prep$pos)-sum(is.na(prep$pos)))/length(prep$pos[,1])-2)

  colnames(me) <- c("time", "mean", "median", "lower CI", "upper CI")
  class(me) <- "popsize"

  return(me)
}