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)