## Load package as R library.
library(qtlbim)

############################################################
## Don't run. Demo of how to create MCMC samples.
## data(hyper)
## hyper <- subset(hyper, chr=1:19) ## Drop X chromosome.
## hyper <- qb.genoprob(hyper, step=2) ## Genoprob with variable stepwidth.
## qbHyper <- qb.mcmc(hyper, pheno.col = 1) ## Time-consuming step.
############################################################

## Get qbHyper with MCMC samples.
data(qbHyper)

## Summaries.
summary(qbHyper)

## Plot selection.
par(ask=TRUE)

## May need to maximize plot size (for qb.coda plot).

plot(qbHyper)
## The following are included in the generic plot:
##    plot(qb.coda(qbHyper))
##    plot(qb.loci(qbHyper))
##    plot(qb.BayesFactor(qbHyper))
##    plot(qb.hpdone(qbHyper))
##    plot(qb.epistasis(qbHyper))
##    plot(qb.diag(qbHyper))

## Plot diagnostics for heritability, environmental variance.
plot(qb.diag(qbHyper, items = c("herit", "envvar")))

## 1-D scan for LPD.
one <- qb.scanone(qbHyper, chr = c(1,4,6,15), type = "LPD")
summary(one)
plot(one)
plot(out.em, chr=c(1,4,6,15), add = TRUE, col = "red", lty = 2)

## Best pattern of QTLs
best <- qb.best(qbHyper)
summary(best)
plot(best)

## Patterns close to the best.
target <- best$model[[1]]
close <- qb.close(qbHyper, target)
summary(close)
trellis.par.set(box.umbrella = list(col = "black"),
  box.rectangle = list(col = "black"), plot.symbol = list(col = "black"))
plot(close)

## Highest posterior density (HPD) pick of best chromosomes.
hpd <- qb.hpdone(qbHyper, profile = "2logBF")
summary(hpd)
plot(hpd)

## Bayes factors for comparing models.
tmp <- qb.BayesFactor(qbHyper)
summary(tmp)
plot(tmp)

## 2-D MCMC scan for LPD.
two <- qb.scantwo(qbHyper, chr = c(6,15),
  type = "2logBF")
plot(two)

## Fancier stuff: 1-D slice of 2-D surface.
plot(two, chr = 6, slice = 15, show.locus = FALSE)
plot(two, chr = 15, slice = 6, show.locus = FALSE)
two.lpd <- qb.scantwo(qbHyper, chr = c(6,15),
  type = "LPD")
plot(two.lpd, chr = 6, slice = 15, show.locus = FALSE)
plot(two.lpd, chr = 15, slice = 6, show.locus = FALSE)

## Slice both ways for detailed examination.
slice <- qb.slicetwo(qbHyper, c(6,15), c(59,19.5))
summary(slice)
plot(slice, figs = c("effects", "cellmean", "effectplot"))