Implementation.R
In CNPBayes: Bayesian mixture models for copy number polymorphisms

## ----lib, message=FALSE----------------------------------------------------
library(CNPBayes)
library(ggplot2)

## ----McmcParams------------------------------------------------------------
mp <- McmcParams(iter=1000, burnin=100, thin=1, nStarts=4)

## ----hyperparameters-------------------------------------------------------
  Hyperparameters()
  HyperparametersMultiBatch()

## ----hplist----------------------------------------------------------------
hp.list <- hpList()
names(hp.list)
hp.list[["SB"]] ## hyperparameters for SB model
hp.list <- hpList(k=3) ## each object will have hyperparameters for a 3-component mixture model
k(hp.list[["SB"]])

## ----graph-marginal, echo=FALSE--------------------------------------------
library(grid)
bayesianGrob <- function(x, r=unit(0.25, "inches")){
  tg <- textGrob(x)
  cg <- circleGrob(r=r)
  boxedText <- gTree(children=gList(tg, cg))
}

grobY <- bayesianGrob(expression(y[i]))
grobThetas <- bayesianGrob(expression(theta[h]))
grobSigma2 <- bayesianGrob(expression(sigma[h]^2))
grobH <- bayesianGrob(expression(z[i]))
grobNu0 <- bayesianGrob(expression(nu[0]))
grobSigma02 <- bayesianGrob(expression(sigma[0]^2))
grobPi <- bayesianGrob(expression(pi[h]))
grobMu <- bayesianGrob(expression(mu[h]))
grobTau2 <- bayesianGrob(expression(tau[h]^2))

h <- unit(0.25, "inches")
e <- unit(0.05, "inches")
d <- unit(0.025, "inches")

ar <- arrow(ends="last", length=unit(0.075, "inches"), type="closed")
grid.newpage()
y.x <- 0.5; y.y <- 0.1
h.x <- 0.1; h.y <- 0.3
theta.x <- 0.5; theta.y <- 0.3
sigma2.x <- 0.7; sigma2.y <- 0.3
pi.x <- 0.1; pi.y <- 0.5
mu.x <- 0.45; mu.y <- 0.5
tau2.x <- 0.55; tau2.y <- 0.5
nu0.x <- 0.65; nu0.y <- 0.5
s20.x <- 0.75; s20.y <- 0.5
grid.draw(editGrob(grobY, vp=viewport(y.x, y.y), gp=gpar(fill="gray")))
grid.draw(editGrob(grobY, vp=viewport(y.x, y.y), gp=gpar(fill="transparent")))
grid.draw(editGrob(grobH, vp=viewport(h.x, h.y)))

grid.draw(editGrob(grobThetas, vp=viewport(theta.x, theta.y)))
grid.draw(editGrob(grobSigma2, vp=viewport(sigma2.x, sigma2.y)))

## theta -> y
grid.move.to(unit(theta.x, "npc"), unit(theta.y, "npc") - h)
grid.line.to(unit(y.x, "npc"), unit(y.y, "npc") + h, arrow=ar,
             gp=gpar(fill="black"))
## sigma2 -> y
grid.move.to(unit(sigma2.x, "npc") - e, unit(sigma2.y, "npc") -h)
grid.line.to(unit(y.x, "npc") + h, unit(y.y, "npc") + h, arrow=ar,
             gp=gpar(fill="black"))

## h -> theta
grid.move.to(unit(h.x, "npc") + h, unit(h.y, "npc") - h)
grid.line.to(unit(y.x, "npc") - h, unit(y.y, "npc") + h,  arrow=ar,
             gp=gpar(fill="black"))

##pi
grid.draw(editGrob(grobPi, vp=viewport(pi.x, pi.y)))
## pi -> h
grid.move.to(x=unit(pi.x, "npc"), y=unit(pi.y, "npc") - h)
grid.line.to(x=unit(h.x, "npc"),
             y=unit(h.y, "npc")+h, arrow=ar,
             gp=gpar(fill="black"))


## mu_h
grid.draw(editGrob(grobMu, vp=viewport(mu.x, mu.y)))
## mu_h -> theta
grid.move.to(x=unit(mu.x, "npc")+e, y=unit(mu.y, "npc") - h)
grid.line.to(x=unit(theta.x, "npc")-e, y=unit(theta.y, "npc")+h, arrow=ar,
             gp=gpar(fill="black"))


## sigma2_h
grid.draw(editGrob(grobTau2, vp=viewport(tau2.x, tau2.y)))
## sigma2_h -> theta_h
grid.move.to(x=unit(tau2.x, "npc")-e, y=unit(tau2.y, "npc") - h)
grid.line.to(x=unit(theta.x, "npc")+e, y=unit(theta.y, "npc")+h, arrow=ar,
             gp=gpar(fill="black"))

grid.draw(editGrob(grobNu0, vp=viewport(nu0.x, nu0.y)))
## nu_0 -> sigma2_0
grid.move.to(x=unit(nu0.x, "npc")+e, y=unit(nu0.y, "npc") - h)
grid.line.to(x=unit(sigma2.x, "npc")-e, y=unit(sigma2.y, "npc")+h, arrow=ar,
             gp=gpar(fill="black"))

grid.draw(editGrob(grobSigma02, vp=viewport(s20.x, s20.y)))
## nu_0 -> sigma2_0
grid.move.to(x=unit(s20.x, "npc")-e, y=unit(s20.y, "npc") - h)
grid.line.to(x=unit(sigma2.x, "npc")+e, y=unit(sigma2.y, "npc")+h, arrow=ar,
             gp=gpar(fill="black"))

## ----simulateData----------------------------------------------------------
set.seed(123)
sim.data <- simulateData(N=1000, p=rep(1/3, 3),
                         theta=c(-1, 0, 1),
                         sds=rep(0.1, 3))

## ----batch-----------------------------------------------------------------
## Create McmcParams for batch model
mp <- McmcParams(iter=600, burnin=100, thin=1, nStarts=4)
model.list <- gibbs(model="SB", dat=y(sim.data),
                    k_range=c(3, 3), mp=mp)

## ----model_selection-------------------------------------------------------
round(sapply(model.list, marginal_lik), 1)
names(model.list)[1]

## ----bic-------------------------------------------------------------------
round(sapply(model.list, bic), 1)

## ----mapComponents---------------------------------------------------------
model <- model.list[[1]]
cn.model <- CopyNumberModel(model)
mapping(cn.model)

## ----mapping---------------------------------------------------------------
ggMixture(cn.model)

## ----probCopyNumber--------------------------------------------------------
head(probCopyNumber(cn.model))

## ----mapping_manual--------------------------------------------------------
cn.model2 <- cn.model
mapping(cn.model2) <- c("1", "1", "2")
head(probCopyNumber(cn.model2))

## ----ggMixture-------------------------------------------------------------
ggMixture(cn.model2)

## ----collapseBatch, eval=FALSE---------------------------------------------
#  k <- 3
#  nbatch <- 3
#  means <- matrix(c(-1.2, -1.0, -1.0,
#                  -0.2, 0, 0,
#                  0.8, 1, 1), nbatch, k, byrow=FALSE)
#  sds <- matrix(0.1, nbatch, k)
#  N <- 500
#  sim.data <- simulateBatchData(N=N,
#                                batch=rep(letters[1:3], length.out=N),
#                                theta=means,
#                                sds=sds,
#                                p=c(1/5, 1/3, 1-1/3-1/5))

## ----MBP, eval=FALSE-------------------------------------------------------
#  gibbs(model="MBP", dat=y(sim.data), k_range=c(3, 3))

## ----fit_many, eval=FALSE--------------------------------------------------
#  model.list <- gibbs(model=c("SBP", "MBP"), k_range=c(1, 5),
#                      dat=y(sim.data),
#                      batches=batch(sim.data), mp=mp)