######################################################################
mcmc.qtlnet <- function(cross, pheno.col, threshold,
addcov=NULL, intcov=NULL,
nSamples = 1000, thinning=1,
max.parents = 3,
M0 = NULL,
burnin = 0.1, method = "hk", random.seed = NULL,
init.edges = 0,
saved.scores = NULL,
rev.method = c("nbhd", "node.edge", "single"),
verbose = FALSE, ...)
{
## Verbose: 1 or TRUE: saved count; 2: MCMC moves; 3: plot BIC; 4: 2&3.
## Check input parameters.
rev.method <- match.arg(rev.method)
## Random number generator seed.
if(!is.null(random.seed)) {
if(!is.numeric(random.seed))
stop("random seed must be numeric")
set.seed(random.seed)
}
## Burnin must be between 0 and 1.
if(is.logical(burnin))
burnin <- ifelse(burnin, 0.1, 0)
if(burnin < 0 | burnin > 1)
stop("burnin must be between 0 and 1")
## Adjust phenotypes and covariates to be numeric.
cross <- adjust.pheno(cross, pheno.col, addcov, intcov)
pheno.col <- cross$pheno.col
pheno.names <- cross$pheno.names
addcov <- cross$addcov
intcov <- cross$intcov
cross <- cross$cross
n.pheno <- length(pheno.col)
## Initial network matrix.
if(is.null(M0))
M0 <- init.qtlnet(pheno.col, max.parents, init.edges)
if(nrow(M0) != n.pheno | ncol(M0) != n.pheno)
paste("M0 must be square matrix the size of pheno.col")
## LOD threshold by phenotype.
if(length(threshold) == 1)
threshold <- rep(threshold, n.pheno)
if(length(threshold) != n.pheno)
stop("threshold must have same length as pheno.col")
## Calculate genotype probabilities if not already done.
if (!("prob" %in% names(cross$geno[[1]]))) {
warning("First running calc.genoprob.")
cross <- qtl::calc.genoprob(cross)
}
Mav <- matrix(0, n.pheno, n.pheno)
n.burnin <- burnin * nSamples * thinning
post.bic <- rep(NA,nSamples)
post.model <- rep(NA,nSamples)
all.bic <- rep(NA,nSamples)
if(is.null(saved.scores))
rev.method <- "single"
saved.scores <- make.saved.scores(pheno.names, max.parents,
saved.scores = saved.scores,
verbose = verbose, ...)
M.old <- M0
ne.old <- nbhd.size(M.old, max.parents)[[1]]
aux.new <- score.model(M.old, saved.scores, cross, addcov, intcov,
threshold, verbose, method = method, ...)
bic.old <- aux.new$model.score
tmp <- aux.new$update.scores
codes <- dimnames(saved.scores)[[1]]
if(!is.null(tmp)) {
index <- nrow(saved.scores)
index <- match(as.character(tmp$code), codes) + index * (tmp$pheno.col - 1)
saved.scores[index] <- tmp$bic
}
model.old <- aux.new$model.name
if(verbose) {
cat("\n")
if(verbose > 2)
graphics::plot(c(1,nSamples), bic.old * c(0.5,1.2), type = "n",
xlab = "sample", ylab = "BIC")
}
k <- 0
if(thinning <= 1) {
post.bic[1] <- all.bic[1] <- bic.old
post.model[1] <- model.old
k <- 1
if(verbose > 2)
graphics::points(k, post.bic[k], cex = 0.5)
}
cont.accept <- numeric(0)
accept.fn <- function(x, move, fate = "accept") {
move <- paste(fate, move, sep = ".")
if(is.na(match(move, names(x))))
x[move] <- 1
else
x[move] <- x[move] + 1
x
}
for(i in 2:(nSamples*thinning)){
## Propose new network structure.
M.new <- if(rev.method == "node.edge")
propose.new.node.edge(M.old, max.parents,
saved.scores = saved.scores, rev.method = rev.method,
verbose = (verbose %in% c(2,4)), ...)
else
propose.new.structure(M.old, max.parents,
saved.scores = saved.scores, rev.method = rev.method,
verbose = (verbose %in% c(2,4)), ...)
rev.ratio <- M.new$rev.ratio
move <- M.new$move ## Want to keep track of moves and acceptance rates.
M.new <- M.new$M
ne.new <- nbhd.size(M.new, max.parents)[[1]]
aux.new <- score.model(M.new, saved.scores, cross, addcov, intcov,
threshold, verbose, method = method, ...)
## Typically only a few update.scores.
tmp <- aux.new$update.scores
if(!is.null(tmp)) {
index <- nrow(saved.scores)
index <- match(as.character(tmp$code), codes) + index * (tmp$pheno.col - 1)
saved.scores[index] <- tmp$bic
}
bic.new <- aux.new$model.score
model.new <- aux.new$model.name
## Accept new model?
## Always if bic.new < bic.old.
## Otherwise do Metropolis-Hastings trick.
if(is.infinite(rev.ratio))
mr <- 1
else
mr <- exp(-0.5*(bic.new - bic.old)) * (ne.old / ne.new) * rev.ratio
if(is.na(mr) | is.null(mr))
browser()
if(stats::runif(1) <= min(1,mr)){
M.old <- M.new
bic.old <- bic.new
ne.old <- ne.new
model.old <- model.new
cont.accept <- accept.fn(cont.accept, move)
}
else
cont.accept <- accept.fn(cont.accept, move, "reject")
## Accumulate M's for average.
if(i > n.burnin)
Mav <- Mav + M.old
## Bookkeeping to save sample.
if((i %% thinning) == 0){
k <- k + 1
all.bic[k] <- bic.new
post.bic[k] <- bic.old
post.model[k] <- model.old
if(verbose) {
print(c(i,k))
if(verbose > 2)
graphics::points(k, post.bic[k], cex = 0.5)
}
}
}
## Make average here.
Mav <- Mav / (nSamples * thinning - n.burnin)
cont.accept <- cont.accept[sort(names(cont.accept))]
out <- list(post.model = post.model,
post.bic = post.bic,
Mav = Mav,
freq.accept = sum(cont.accept) / (nSamples * thinning),
cont.accept = cont.accept,
saved.scores=saved.scores,
all.bic=all.bic,
cross = cross)
## Attributes of qtlnet object.
attr(out, "M0") <- M0
attr(out, "threshold") <- threshold
attr(out, "nSamples") <- nSamples
attr(out, "thinning") <- thinning
attr(out, "pheno.col") <- pheno.col
attr(out, "pheno.names") <- pheno.names
attr(out, "addcov") <- addcov
attr(out, "intcov") <- intcov
attr(out, "burnin") <- burnin
attr(out, "method") <- method
attr(out, "random.seed") <- random.seed
attr(out, "random.kind") <- RNGkind()
class(out) <- c("qtlnet","list")
out
}
######################################################################
c.qtlnet <- function(...)
{
## Combine qtlnet objects.
netlist <- list(...)
out <- netlist[[1]]
if(!inherits(out, "qtlnet")) {
netlist <- out
out <- netlist[[1]]
}
if(!inherits(out, "qtlnet"))
stop("argument must be list of qtlnet objects")
burnin <- attr(out, "burnin")
n.pheno <- length(attr(out, "pheno.col"))
if(length(netlist) > 1) {
## Need to do some checking here that qtlnet objects match up.
## Minimal for now.
## Assumes that parameters stay the same, although nSamples could change.
if(n.pheno != mean(sapply(netlist, function(x) length(attr(x, "pheno.col")))))
stop("different numbers of phenotypes not allowed")
if(attr(out, "thinning") != mean(sapply(netlist, function(x) attr(x, "thinning"))))
stop("thinning values differ")
## check threshold vector.
if(!all(attr(out, "threshold") == apply(sapply(netlist,
function(x) attr(x, "threshold")), 1, mean)))
stop("threshold values differ")
if(!all(burnin == sapply(netlist, function(x) attr(x, "burnin"))))
stop("burnin values differ")
if(is.matrix(out$Mav)) {
tmp <- out$Mav
nsheet <- 1
out$Mav <- array(0, c(nrow(tmp), ncol(tmp), length(netlist)))
}
else {
## Already is array. Want to add another sheet.
tmp <- out$Mav
dtmp <- dim(tmp)
nsheet <- dtmp[3]
dtmp[3] <- nsheet + length(netlist) - 1
out$Mav <- array(0, dtmp)
}
out$Mav[, , seq(nsheet)] <- tmp
out$cont.accept <- list(out$cont.accept)
for(i in seq(2, length(netlist))) {
## Per step summaries.
for(j in c("post.model","post.bic","all.bic","freq.accept"))
out[[j]] <- c(out[[j]], netlist[[i]][[j]])
out$cont.accept[[i]] <- netlist[[i]]$cont.accept
## Attributes.
attr(out, "nSamples") <- c(attr(out, "nSamples"),
attr(netlist[[i]], "nSamples"))
## Matrices of average network structures.
out$Mav[ , , nsheet - 1 + i] <- netlist[[i]]$Mav
}
}
out
}
######################################################################
nbhd.size <- function(M, max.parents = 3)
{
n.deletions <- sum(M)
n.additions <- 0
n.reversions <- 0
le <- ncol(M)
for(j in 1:le){
add.forbid <- forbidden.additions(M, j, max.parents)
n.additions <- n.additions + (le - 1 - length(c(add.forbid$upf, add.forbid$downf)))
rev.allow <- check.reversions(M, j, max.parents)
if(!is.null(rev.allow))
n.reversions <- n.reversions + nrow( rev.allow )
}
nbhd.size <- n.deletions + n.additions + n.reversions
list(nbhd.size=nbhd.size,
n.deletions=n.deletions,
n.additions=n.additions,
n.reversions=n.reversions)
}######################################################################
init.qtlnet <- function(pheno.col, max.parents = 3, init.edges = NULL)
{
n.pheno <- length(pheno.col)
if(is.null(init.edges)) {
n.edges <- n.pheno * (n.pheno - 1) / 2
init.edges <- sample(seq(0, n.edges), 1)
}
M <- matrix(0, n.pheno, n.pheno)
if(init.edges > 0)
for(i in seq(init.edges)) {
cause <- sample(n.pheno, 1)
if(i == 1)
effect <- sample(seq(n.pheno)[-cause], 1)
else
effect <- propose.add(M, cause, max.parents)
if(effect == 0)
break
M[cause, effect] <- 1
}
M
}
######################################################################
propose.add <- function(M, node, max.parents)
{
## Add an edge with a direction.
aux1 <- forbidden.additions(M, node, max.parents)
## Is there a down option?
aux3 <- unique(c(node, aux1$upf, aux1$downf))
aux3 <- seq(ncol(M))[-aux3]
if(length(aux3)) {
## Check if any of aux3 is at or above max.parents.
wh <- which(apply(M[, aux3, drop = FALSE], 2, sum) >= max.parents)
if(length(wh))
aux3 <- aux3[-wh]
}
if(length(aux3) == 0)
aux3 <- 0
else {
if(length(aux3) > 1)
aux3 <- sample(aux3, 1)
}
aux3
}
######################################################################
propose.new.node.edge <- function(M, max.parents = 3,
saved.scores, rev.method = "node.edge",
propose = c(node = 1, edge = 2, reverse = 10, drop = 2),
check.down = FALSE,
verbose = FALSE, ...)
{
## Ref: Grzegorczyk and Husmeier (2008) Mach Learn 71: 265-305.
## Extension is to updating individual nodes and dropping edges.
rev.ratio <- 1
## Proposal weights for types of changes.
if(any(propose <= 0))
stop("propose values must be positive")
name.propose <- names(propose)
propose <- array(propose, 4)
if(length(name.propose) == 4)
names(propose) <- name.propose
else
names(propose) <- c("node","edge","reverse","drop")
## To speed up, use M as logical matrix. Still return 0/1 matrix for now.
le.nodes <- ncol(M)
le.edges <- sum(M)
prob.node <- propose["node"] * le.nodes
prob.node <- prob.node / (prob.node + propose["edge"] * le.edges)
if(stats::runif(1) <= prob.node) {
move <- "node"
## Propose new parents for a node.
node <- sample(seq(le.nodes), 1)
if(verbose)
cat(move, "")
aux2 <- update.node(M, max.parents, saved.scores, node, verbose)
rev.ratio <- aux2$rev.ratio
M <- aux2$M
if(verbose) {
up <- node.parents(M, node)$parents
cat(node, "up:", up, "\n")
}
}
else {
## Propose change to edge.
move <- "edge"
if(verbose)
cat(move, "")
## Pick a valid edge.
edge <- which(M == 1)
if(length(edge) > 1)
edge <- sample(edge, 1)
edge <- c(row(M)[edge], col(M)[edge])
if(verbose)
cat("edge ")
## Cyclic mistake checks.
## These should not happen!
if(edge[1] == edge[2]) {
cat("edge identity:", edge, "\n")
browser()
}
aux1 <- check.downstream(M, edge[2])
if(any(aux1 == edge[1])) {
cat("edge downfall:", edge, aux1, "\n")
browser()
}
prob.reverse <- propose["reverse"]
prob.reverse <- prob.reverse / (prob.reverse + propose["drop"])
if(stats::runif(1) <= prob.reverse) {
## Propose to reverse an edge.
if(verbose) cat("reverse ")
aux2 <- rev.edge(M, max.parents, saved.scores, edge, verbose)
}
else {
## Propose to drop an edge.
if(verbose) cat("drop ")
aux2 <- drop.edge(M, max.parents, saved.scores, edge, verbose)
}
rev.ratio <- aux2$rev.ratio
M <- aux2$M
if(verbose) {
for(i in 1:2) {
up <- node.parents(M, edge[i])$parents
cat(edge[i], "up:", up, "\n")
}
}
if(check.down) {
## NOTE: This takes time and should be dropped eventually.
## Check if we somehow have a created a cycle earlier. This should not happen.
down <- check.downstream(M, edge[2])[-1]
up <- check.upstream(M, edge[2])[-1]
if(length(down) & length(up)) {
if(any(down %in% up)) {
cat("downup:", edge, "\n", sort(down), "\n", sort(up), "\n")
browser()
}
}
}
}
list(M = M, rev.ratio = rev.ratio, move = move)
}
######################################################################
propose.new.structure <- function(M, max.parents = 3,
saved.scores, rev.method = "nbhd",
verbose = FALSE, ...)
{
## Acceptance rate is <20%. Could we improve here?
## Yes. See Grzegorczyk and Husmeier (2008) Mach Learn 71: 265-305.
rev.ratio <- 1
## To speed up, use M as logical matrix. Still return 0/1 matrix for now.
le.nodes <- ncol(M)
flag <- TRUE
while(flag){
## Pick node and decide on add/delete/reverse.
## Keep doing this until successful.
node <- sample(seq(le.nodes),1)
move <- sample(c("add","delete","reverse"),1)
if(verbose)
cat(node, move, "")
switch(move,
add = {
aux3 <- propose.add(M, node, max.parents)
if(!(flag <- (aux3 == 0))) {
aux1 <- check.downstream(M, aux3)
if(any(aux1 == node)) {
## This should not happen!
cat(move, "downfall:", node, aux1, "\n")
browser()
}
M[node,aux3] <- 1
if(verbose)
cat(aux3, "\n")
}
},
reverse = {
## Reverse direction of an edge.
if(flag <- (rev.method == "single")) {
aux1 <- check.reversions(M, node, max.parents)
if(!(flag <- is.null(aux1))) {
le.rev <- nrow(aux1)
aux3 <- sample(seq(le.rev), 1)
aux3 <- aux1[aux3, ]
M[aux3[1], aux3[2]] <- 0
M[aux3[2], aux3[1]] <- 1
}
}
else { ## Reverse method of Grzegorczyk and Husmeier.
aux1 <- c(which(M[,node] == 1), le.nodes + which(M[node,] == 1))
if(!(flag <- !length(aux1))) {
if(length(aux1) > 1)
aux1 <- sample(aux1, 1)
if(aux1 > le.nodes)
aux3 <- c(node, aux1 - le.nodes)
else
aux3 <- c(aux1, node)
aux1 <- check.downstream(M, aux3[2])
if(any(aux1 == aux3[1])) {
## This should not happen!
cat(move, "downfall:", aux3[1], aux1, "\n")
browser()
}
}
if(!flag) {
if(aux3[1] == aux3[2]) {
## This should not happen!
cat(move, "identity:", aux3, "\n")
browser()
}
aux2 <- rev.edge(M, max.parents, saved.scores, aux3, verbose)
rev.ratio <- aux2$rev.ratio
M <- aux2$M
}
}
if(!flag & verbose)
cat(aux3[1], aux3[2], "\n")
},
delete = {
## Delete an existing edge through node.
aux1 <- c(which(M[, node] == 1), le.nodes + which(M[node,] == 1))
if(!(flag <- !length(aux1))) {
if(length(aux1) > 1)
aux1 <- sample(aux1, 1)
if(aux1 > le.nodes)
aux3 <- c(node, aux1 - le.nodes)
else
aux3 <- c(aux1, node)
M[aux3[1], aux3[2]] <- 0
if(verbose)
cat(aux3, "\n")
}
})
}
## Check if we somehow have a created a cycle earlier. This should not happen.
down <- check.downstream(M, aux3[2])[-1]
up <- check.upstream(M, aux3[2])[-1]
if(length(down) & length(up)) {
if(any(down %in% up)) {
cat(move, "downup:", aux3, "\n", sort(down), "\n", sort(up), "\n")
browser()
}
}
list(M = M, rev.ratio = rev.ratio, move = move)
}
######################################################################
update.node <- function(M, max.parents, saved.scores, node, verbose = FALSE)
{
## Scores for possible parent sets.
down <- check.downstream(M, node)[-1]
index <- index.parents(saved.scores, down, node)
if(length(index))
z.1 <- saved.scores[-index, node]
else
z.1 <- saved.scores[, node]
z.1 <- exp(min(z.1) - z.1)
s.1 <- sum(z.1)
if(length(z.1) == 0) ## Should not happen
browser()
## Find probability for current parents of node.
rev.ratio <- z.1[find.parent.score(M, saved.scores, -index, node)]
## Make node an orphan.
M[,node] <- 0
## Sample new parents of node.
parent <- sample(seq(length(z.1)), 1, prob = z.1)
## Proposal ratio.
rev.ratio <- z.1[parent] / rev.ratio
## Add edges to graph for new parents of node.
new.parent <- find.index.parent(M, saved.scores, -index, node, parent)
if(length(new.parent))
M[new.parent, node] <- 1
if(verbose) {
cat("node rev.ratio", rev.ratio, "\n")
if(is.na(rev.ratio))
browser()
}
list(M = M, rev.ratio = rev.ratio)
}
######################################################################
drop.edge <- function(M, max.parents, saved.scores, node.pair, verbose = FALSE)
{
## Grzegorczyk and Husmeier (2008) Mach Learn 71: 265-305.
## Call provides node.pair[1 -> 2].
## Step 1: Orphan both nodes after computing reverse proposal prob.
rev.ratio <- reverse.proposal(M, saved.scores, node.pair, verbose)
M[, node.pair] <- 0
## Step 2. Sample new parents for node 1 without node 2 as one parent.
## Exclude parents downstream of 1 or 2.
tmp <- sample.exclude.down(M, saved.scores, node.pair, verbose)
q.1 <- tmp$prob
parent <- tmp$parent
## Add edges to graph for parents of node 1.
if(length(parent))
M[parent, node.pair[1]] <- 1
## Step 3. Sample new parents for node 2 without node 1 as parent.
## Exclude parents downstream of 1 or 2.
tmp <- sample.exclude.down(M, saved.scores, rev(node.pair), verbose)
q.2 <- tmp$prob
parent <- tmp$parent
## Proposal ratio.
rev.ratio <- q.1 * q.2 / rev.ratio
## Add edges to graph for parents of node 2.
if(length(parent))
M[parent, node.pair[2]] <- 1
if(verbose) {
cat("edge rev.ratio", rev.ratio, "\n")
if(is.na(rev.ratio))
browser()
}
list(M = M, rev.ratio = rev.ratio)
}
######################################################################
sample.exclude.down <- function(M, saved.scores, node.pair, verbose = FALSE)
{
## Exclude parents downstream of 1 or 2.
down <- unique(c(check.downstream(M, node.pair[2]),
check.downstream(M, node.pair[1])[-1]))
index <- index.parents(saved.scores, down, node.pair[1])
z.1 <- saved.scores[-index, node.pair[1]]
z.1 <- exp(min(z.1) - z.1)
if(length(z.1) == 0) ## Should not happen
browser()
parent <- sample(seq(length(z.1)), 1, prob = z.1)
q.1 <- z.1[parent] / sum(z.1)
new.parent <- find.index.parent(M, saved.scores, -index, node.pair[1], parent)
list(parent = new.parent, prob = q.1)
}
######################################################################
reverse.proposal <- function(M, saved.scores, node.pair, verbose = FALSE)
{
## Parents of 2 include 1 but exclude all downstream of 2.
index <- index.parents(saved.scores, node.pair[1], node.pair[2])
down <- check.downstream(M, node.pair[2])[-1]
if(length(down))
index <- index[-index.parents(saved.scores[index,, drop = FALSE], down, node.pair[2])]
z.2 <- saved.scores[index, node.pair[2]]
z.2 <- exp(min(z.2) - z.2)
q.2 <- z.2[find.parent.score(M, saved.scores, index, node.pair[2])] / sum(z.2)
## Parents of 1 exclude nodes at or downstream of 2.
down <- unique(c(check.downstream(M, node.pair[2]),
check.downstream(M, node.pair[1])[-1]))
index <- index.parents(saved.scores, down, node.pair[1])
z.1 <- saved.scores[-index, node.pair[1]]
z.1 <- exp(min(z.1) - z.1)
q.1 <- z.1[find.parent.score(M, saved.scores, -index, node.pair[1])] / sum(z.1)
rev.ratio <- q.1 * q.2
if(verbose) {
cat("rev.ratio", rev.ratio, "\n")
if(is.na(rev.ratio))
browser()
}
rev.ratio
}
######################################################################
rev.edge <- function(M, max.parents, saved.scores, node.pair, verbose = FALSE)
{
## Grzegorczyk and Husmeier (2008) Mach Learn 71: 265-305.
## Call provides node.pair[1 -> 2].
## Step 1: Orphan both nodes after computing reverse proposal prob.
rev.ratio <- reverse.proposal(M, saved.scores, node.pair, verbose)
M[, node.pair] <- 0
## Step 2. Sample new parents for node 1 with node 2 as one parent.
## Exclude parents downstream of 1 (except 2).
index <- index.parents(saved.scores, node.pair[2], node.pair[1])
down <- check.downstream(M, node.pair[1])[-1]
if(length(down))
index <- index[-index.parents(saved.scores[index,, drop = FALSE], down, node.pair[1])]
z.1 <- saved.scores[index, node.pair[1]]
z.1 <- exp(min(z.1) - z.1)
if(length(z.1) == 0) ## Should not happen
browser()
parent <- sample(seq(length(z.1)), 1, prob = z.1)
q.1 <- z.1[parent] / sum(z.1)
new.parent <- find.index.parent(M, saved.scores, index, node.pair[1], parent)
## Add edges to graph for parents of node 1.
if(length(new.parent))
M[new.parent, node.pair[1]] <- 1
## Step 3. Sample new parents for node 2 without node 1 as parent.
## Exclude parents downstream of 1 or 2.
## Somehow this is leading to cycles!
## SOmehow the rev(node.pair) is not working.
tmp <- sample.exclude.down(M, saved.scores, rev(node.pair))
q.2 <- tmp$prob
new.parent <- tmp$parent
## Add edges to graph for parents of node 2.
if(length(new.parent))
M[new.parent, node.pair[2]] <- 1
rev.ratio <- q.1 * q.2 / rev.ratio
if(verbose) {
cat("edge rev.ratio", rev.ratio, "\n")
if(is.na(rev.ratio))
browser()
}
list(M = M, rev.ratio = rev.ratio)
}
######################################################################
find.index.parent <- function(M, saved.scores, index, node, new.parent)
{
if(length(index))
parents <- dimnames(saved.scores)[[1]][index][new.parent]
else
parents <- dimnames(saved.scores)[[1]][new.parent]
unlist(sapply(strsplit(parents, ",", fixed = TRUE),
function(x, node) {
x <- as.numeric(x)
x + (x >= node)
}, node))
}
######################################################################
index.parents <- function(saved.scores, node1, node2)
{
## Row names of saved.scores renumber around missing node.
parents <- node1 - (node2 < node1)
## Which parents include node.pair[2]?
which(sapply(strsplit(dimnames(saved.scores)[[1]], ",", fixed = TRUE),
function(x, parents) any(x %in% parents),
parents))
}
######################################################################
find.parent.score <- function(M, saved.scores, index, node)
{
parents <- node.parents(M, node)$parents
if(!is.null(parents)) {
parents <- parents - (parents > node)
if(length(index))
match(paste(parents, collapse = ","), dimnames(saved.scores)[[1]][index])
else
match(paste(parents, collapse = ","), dimnames(saved.scores)[[1]])
}
else
1
}
######################################################################
forbidden.additions <- function(M, node, max.parents = 3)
{
## upf are forbidden upstream nodes (already present downstream).
## downf are forbidden downstream nodes (already present upstream).
## Check on cycles is one step. See check.reversions() for longer cycles.
## Forbidden upstream additions
upf <- which(M[node,] == 1)
if(length(upf)==0)
upf <- NULL
## Forbidden downstream additions. More complicated.
downf <- which(M[,node] == 1)
le <- length(downf)
## If at least max.parents are causal for node, forbid any more.
if(le >= max.parents)
downf <- seq(nrow(M))[-node]
else {
if(le > 0)
downf <- check.upstream(M, downf)
else
downf <- NULL
}
list(upf=upf, downf=downf)
}
######################################################################
check.qtlnet <- function(object,
min.prob = 0.9,
correct = TRUE,
verbose = FALSE,
...)
{
pheno.names <- attr(object, "pheno.names")
n.pheno <- length(pheno.names)
forbid <- 1
while(!is.null(forbid)) {
M <- threshold.net(object, min.prob = min.prob, ...)
min.prob <- attr(M, "min.prob")
M1 <- 1 * (M > 0)
## This may not be right yet.
forbid <- NULL
for(i in seq(n.pheno)) {
downf <- check.upstream(M1, i)
wh <- which(M1[i, downf] == 1)
if(length(wh))
forbid <- rbind(forbid, cbind(downf[wh],i, M[i, downf[wh]]))
}
if(!is.null(forbid)) {
forbid <- data.frame(forbid, stringsAsFactors = TRUE)
names(forbid) <- c("cause","react","prob")
for(i in 1:2)
forbid[[i]] <- ordered(pheno.names[forbid[[i]]], pheno.names)
}
if(verbose)
print(forbid)
}
attr(M, "min.prob") <- min.prob
list(forbid = forbid, M = M)
}
######################################################################
check.upstream <- function(M, nodes)
{
## After accounting for scanone, 85% of time is nbhd.size.
## Of that, 85% (75% overall) is in check.upstream.
count.upok <- nrow(M) - length(nodes)
## check upstream to see if a cycle would be created.
is.up <- apply(M[, nodes, drop = FALSE], 1, sum)
flag <- TRUE
while(flag){
aux1 <- which(is.up > 0)
if(flag <- (length(aux1) > 0 & count.upok > 0)) {
aux1 <- aux1[!(aux1 %in% nodes)]
if(flag <- length(aux1)) {
is.up <- apply(M[, aux1, drop = FALSE], 1, sum)
nodes <- c(nodes, aux1)
count.upok <- count.upok - flag
}
}
}
nodes
}
######################################################################
check.downstream <- function(M, nodes)
{
count.downok <- nrow(M) - length(nodes)
## check downstream to see if a cycle would be created.
is.down <- apply(M[nodes,, drop = FALSE], 2, sum)
flag <- TRUE
while(flag){
aux1 <- which(is.down > 0)
if(flag <- (length(aux1) > 0 & count.downok > 0)) {
aux1 <- aux1[!(aux1 %in% nodes)]
if(flag <- length(aux1)) {
is.down <- apply(M[aux1,, drop = FALSE], 2, sum)
nodes <- c(nodes, aux1)
count.downok <- count.downok - flag
}
}
}
nodes
}
######################################################################
check.reversions <- function(M, node, max.parents = 3)
{
## Check on possible reversals of directions.
## allowed if no cycles produced.
## forbidden if cycles result.
allowed <- NULL
## Check upstream.
up <- which(M[,node] == 1)
le <- length(up)
if(le) {
if(le == 1)
allowed <- cbind(up, node)
else { ## le > 1
## Multiple colliders.
forbid.up <- rep(FALSE, le)
for(k in 1:le)
forbid.up[k] <- up[k] %in% check.upstream(M, up[-k])
if(any(forbid.up)){
if(any(!forbid.up))
allowed <- cbind(up[!forbid.up], node)
}
else
allowed <- cbind(up, node)
}
}
## Check downstream.
down <- which(M[node,] == 1)
le <- length(down)
allowed2 <- NULL
if(le) {
if(le == 1)
allowed2 <- cbind(node, down)
else { ## le > 1
## Multiple colliders.
forbid.down <- rep(FALSE, le)
for(k in 1:le)
forbid.down[k] <- down[k] %in% check.downstream(M, down[-k])
if(any(forbid.down)){
if(any(!forbid.down))
allowed2 <- cbind(node, down[!forbid.down])
}
else
allowed2 <- cbind(down, node)
}
}
allowed <- rbind(allowed, allowed2)
## Final check of max.parents.
if(!is.null(allowed)) {
wh <- which(apply(M[, allowed[,1], drop = FALSE], 2, sum) >= max.parents)
if(length(wh)) {
## Remove pairs.
allowed <- allowed[-wh,, drop = FALSE]
if(nrow(allowed) == 0)
allowed <- NULL
}
}
allowed
}
##########################################################################
legacy.qtlnet <- function(qtlnet.object, codes = TRUE)
{
## Convert old qtlnet object to current format. Idempotent.
qtlnet.object$Mav <- get.model.average(qtlnet.object)
qtlnet.object$post.net.str <- NULL
qtlnet.object$model.average <- NULL
if(codes)
dimnames(qtlnet.object$saved.scores)[[1]] <-
legacy.code(dimnames(qtlnet.object$saved.scores)[[1]])
qtlnet.object
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.