bic: Pre-compute BIC values for qtlnet sampling.

In qtlnet: Causal Inference of QTL Networks

Description

Pre-compute BIC values for qtlnet sampling to speed up MCMC sampling.

Usage

bic.qtlnet(cross, pheno.col, threshold, addcov = NULL, intcov = NULL,
  max.parents = 3, parents, verbose = TRUE, ...)
bic.join(cross, pheno.col, ..., max.parents = 3)
data(Pscdbp.bic)

Arguments

cross	Object of class cross. See read.cross.
pheno.col	Phenotype identifiers from cross object. May be numeric, logical or character.
threshold	Scalar or list of thresholds, one per each node.
addcov	Additive covariates for each phenotype (NULL if not used). If entered as scalar or vector (same format as pheno.col), then the same addcov is used for all phenotypes. Altenatively, may be a list of additive covariate identifiers.
intcov	Interactive covariates, entered in the same manner as addcov.
max.parents	Maximum number of parents per node. This reduces the complexity of graphs and shortens run time. Probably best to consider values of 3-5.
parents	List containing all possible parents up to max.parents in size. May be a subset
verbose	Print iteration and number of models fit.
...	Additional arguments passed to internal routines. In the case of bic.join, these are a list of objects produced by bic.qtlnet (see example below).

Details

The most expensive part of calculations is running scanone on each phenotype with parent phenotypes as covariates. One strategy is to pre-compute the BIC contributions using a cluster and save them for later use.

We divide the job into three steps: 1) determine parents and divide into reasonable sized groups; 2) compute BIC scores using scanone on a grid of computers; 3) compute multiple MCMC runs on a grid of computers. See the example for details.

Value

Matrix with columns:

code	Binary code as decimal for the parents of a phenotype node, excluding the phenotype. Value is between 0 (no parents) and 2 ^ (length(pheno.col) - 1).
pheno.col	Phenotype column in reduced set, in range 1:length(pheno.col).
bic	BIC score for phenotype conditional on parents (and covariates).

Author(s)

Brian S. Yandell and Elias Chaibub Neto

References

Chaibub Neto E, Keller MP, Attie AD, Yandell BS (2010) Causal Graphical Models in Systems Genetics: a unified framework for joint inference of causal network and genetic architecture for correlated phenotypes. Ann Appl Statist 4: 320-339. http://dx.doi.org/10.1214/09-AOAS288

See Also

mcmc.qtlnet, parents.qtlnet

Examples

pheno.col <- 1:13
max.parents <- 12
size.qtlnet(pheno.col, max.parents)

## Not run: 
## Compute all phenotype/parent combinations.
## This shows how to break up into many smaller jobs.

#########################################
## STEP 1: Preparation. Fast. Needed in steps 2 and 3.

pheno.col <- 1:13
max.parents <- 12
threshold <- 3.83

## Load cross object. Here we use internal object.
data(Pscdbp)
## or: load("Pscdbp.RData")
cross <- Pscdbp
## or: cross <- read.cross("Pscdbp.csv", "csv")

## Break up into groups to run on several machines.
## ~53 groups of ~1000, for a total of 53248 scanone runs.
parents <- parents.qtlnet(pheno.col, max.parents)
groups <- group.qtlnet(parents = parents, group.size = 1000)

## Save all relevant objects for later steps.
save(cross, pheno.col, max.parents, threshold, parents, groups,
  file = "Step1.RData", compress = TRUE)

#########################################
## STEP 2: Compute BIC scores. Parallelize.

## NB: Configuration of parallelization determined using Step 1 results.
## Load Step 1 computations.
load("Step1.RData")

## Parallelize this:
for(i in seq(nrow(groups)))
{
  ## Pre-compute BIC scores for selected parents.
  bic <- bic.qtlnet(cross, pheno.col, threshold,
    max.parents = max.parents,
    parents = parents[seq(groups[i,1], groups[i,2])])

  save(bic, file = paste("bic", i, ".RData", sep = ""), compress = TRUE)
}

#########################################
## STEP 3: Sample Markov chain (MCMC). Parallelize.

## NB: n.runs sets the number of parallel runs.
n.runs <- 100

## Load Step 1 computations.
load("Step1.RData")

## Read in saved BIC scores and combine into one object.
bic.group <- list()
for(i in seq(nrow(groups)))
{
  load(paste("bic", i, ".RData", sep = ""))
  bic.group[[i]] <- bic
}
saved.scores <- bic.join(cross, pheno.col, bic.group)


## Parallelize this:
for(i in seq(n.runs))
{
  ## Run MCMC with randomized initial network.
  mcmc <- mcmc.qtlnet(cross, pheno.col, threshold = threshold,
    max.parents = max.parents, saved.scores = saved.scores, init.edges = NULL)

  save(mcmc, file = paste("mcmc", i, ".RData", sep = ""), compress = TRUE)
}

#########################################
## STEP 4: Combine results for post-processing.

## NB: n.runs needed from Step 3.
n.runs <- 100

## Combine outputs together.
outs.qtlnet <- list()
for(i in seq(n.runs))
{
  load(paste("mcmc", i, ".RData", sep = ""))
  outs.qtlnet[[i]] <- mcmc
}
out.qtlnet <- c.qtlnet(outs.qtlnet)
summary(out.qtlnet)
print(out.qtlnet)

## End of parallel example.
#########################################

## End(Not run)

dim(Pscdbp.bic)

qtlnet documentation built on April 14, 2020, 6:24 p.m.