R/marginal.posterior.profile.R

In mcveanlab/TreeWASDir: Functions for TreeWAS on -/0/+ states

#' Function to get marginal posterior on -/0/+ profile
#'
#' @param pars A list object as returned by the function set.parameters
#' @param data.sub A 3-dimension array (N. SNPs x N. phenotypes x 3) with state likelihoods.
#' @return pp. A 2-dimension array (N. Phenotypes x 3) with posterior probabilties.
#' 

marginal.posterior.profile <- function(
    pars = NULL,
    data.sub = NULL
) {

    ## Get forward and G matrices
    tmp <- calc.llk.tree(pars, data.sub, returnForwardMatrix=T, returnGMatrix=T);
    f <- tmp$f;
    g <- tmp$g;
    
    ## Build parents list and reverse F traversal order
    parents <- rep(0, pars$n.phenos);
    for (i in pars$t.path) parents[pars$ontology[[i]]] <- i;
    ord <- c(rev(pars$t.path), pars$terminals);

    ## Construct backward matrix
    b <- array(0, dim(f));
    b[ord[1],] <- log(pars$stat.dist);
    for (i in ord[-1]) {
        r.i <- b[parents[i],]+f[parents[i],]-g[i,];
        mx <- max(r.i);
        tmp <- mx+log(sum(exp(r.i-mx)));
        tmp <- tmp+pars$lp.switch+log(pars$stat.dist);
        tmp2 <- -pars$theta.tree+b[parents[i],]+f[parents[i],]-g[i,];
        tmp3 <- cbind(tmp2, tmp);
        mx <- apply(tmp3, 1, max);
        b[i,] <- log(rowSums(exp(tmp3-mx)))+mx;
    }
    
    ## Posteriors
    tmp <- b+f;
    mx <- apply(tmp, 1, max);
    pp <- exp(tmp-mx);
    pp <- pp/rowSums(pp);
    return(pp);
}