## ----lib-----------------------------------------------------------------
library(CNPBayes)
library(GenomicRanges)

## ----find_cnps-----------------------------------------------------------
se <- readRDS(system.file("extdata", "simulated_se.rds", package="CNPBayes"))
grl <- readRDS(system.file("extdata", "grl_deletions.rds", package="CNPBayes"))

## ----summary, message=FALSE----------------------------------------------
cnv.region <- consensusCNP(grl, max.width=5e6)
i <- subjectHits(findOverlaps(cnv.region, rowRanges(se)))
med.summary <- matrixStats::colMedians(assays(se)[["cn"]][i, ], na.rm=TRUE)

## ----model_construction--------------------------------------------------
model.marginal <- MarginalModel(data=med.summary)
model.batch <- BatchModel(data=med.summary, 
                          batch=c(rep(1, 12), rep(2, 23)))

## ----collapse------------------------------------------------------------
model.batch <- BatchModel(data=med.summary,
                          batch=collapseBatch(model.batch))

## ----posteriorSimulation-------------------------------------------------
set.seed(1337)
mlist.marginal <- posteriorSimulation(model.marginal, k=1:4)
mlist.batch <- posteriorSimulation(model.batch, k=1:4)

## ----likelihood----------------------------------------------------------
# marginal likelihood of each model
marginal.lik <- marginalLikelihood(mlist.marginal)
batch.lik <- marginalLikelihood(mlist.batch)

# bayes factor for comparing top two models
logBayesFactor(marginal.lik)
logBayesFactor(batch.lik)

## ----map-----------------------------------------------------------------
selected.model <- mlist.marginal[[which.max(marginal.lik)]]
map(selected.model)