Nothing
#' Create an htc identification function.
#'
#' A helper function for \code{\link{htcIdentifyStep}}, creates an identifier
#' function based on its given parameters. This created identifier function will
#' identify the directed edges from 'targets' to 'node.'
#'
#' @param idFunc identification of edge coefficients often requires that other
#' edge coefficients already be identified. This argument should be a
#' function that produces all such identifications. The newly created
#' identifier function will return these identifications along with its
#' own.
#' @param sources the sources of the half-trek system.
#' @param targets the targets of the half-trek system (these should be the
#' parents of \code{node}).
#' @param node the node for which all incoming edges are to be identified
#' (the tails of which are targets).
#' @param htrSources the nodes in sources which are half-trek reachable from
#' \code{node}. All incoming edges to these sources should be identified by
#' \code{idFunc} for the newly created identification function to work.
#'
#' @return an identification function
#'
#' @references
#' Foygel, R., Draisma, J., and Drton, M. (2012) Half-trek criterion for
#' generic identifiability of linear structural equation models.
#' \emph{Ann. Statist.} 40(3): 1682-1713.
createHtcIdentifier <- function(idFunc, sources, targets, node, htrSources) {
# Necessary redundant assignments
idFunc <- idFunc
sources <- sources
targets <- targets
node <- node
htrSources <- htrSources
return(function(Sigma) {
m <- nrow(Sigma)
identifiedParams <- idFunc(Sigma)
Lambda <- identifiedParams$Lambda
SigmaMinus <- Sigma
for (source in htrSources) {
SigmaMinus[source, ] <- Sigma[source, ] - t(Lambda[, source, drop = F]) %*%
Sigma
}
if (abs(det(SigmaMinus[sources, targets, drop = F])) < 10^-10) {
stop("In identification, found near-singular system. Is the input matrix generic?")
}
Lambda[targets, node] <- solve(SigmaMinus[sources, targets, drop = F], SigmaMinus[sources,
node, drop = F])
if (!any(is.na(Lambda))) {
Omega <- t(diag(m) - Lambda) %*% Sigma %*% (diag(m) - Lambda)
return(list(Lambda = Lambda, Omega = Omega))
}
return(list(Lambda = Lambda, Omega = identifiedParams$Omega))
})
}
#' Perform one iteration of HTC identification.
#'
#' A function that does one step through all the nodes in a mixed graph
#' and tries to identify new edge coefficients using the existence of
#' half-trek systems as described in Foygel, Draisma, Drton (2012).
#'
#' @param mixedGraph a \code{\link{MixedGraph}} object representing
#' the mixed graph.
#' @param unsolvedParents a list whose ith index is a vector of all the parents
#' j of i in G which for which the edge j->i is not yet known to be
#' generically identifiable.
#' @param solvedParents the complement of \code{unsolvedParents}, a list whose
#' ith index is a vector of all parents j of i for which the edge i->j
#' is known to be generically identifiable (perhaps by other algorithms).
#' @param identifier an identification function that must produce the
#' identifications corresponding to those in solved parents. That is
#' \code{identifier} should be a function taking a single argument Sigma
#' (any generically generated covariance matrix corresponding
#' to the mixed graph) and returns a list with two named arguments
#' \describe{
#' \item{Lambda}{denote the number of nodes in \code{mixedGraph} as n. Then
#' Lambda is an nxn matrix whose i,jth entry
#' \enumerate{
#' \item equals 0 if i is not a parent of j,
#' \item equals NA if i is a parent of j but \code{identifier} cannot
#' identify it generically,
#' \item equals the (generically) unique value corresponding to the weight
#' along the edge i->j that was used to produce Sigma.
#' }}
#' \item{Omega}{just as Lambda but for the bidirected edges in the mixed
#' graph}
#' }
#' such that if j is in \code{solvedParents[[i]]} we must have that
#' Lambda[j,i] is not NA.
#'
#' @return a list with four components:
#' \describe{
#' \item{\code{identifiedEdges}}{a matrix rx2 matrix where r is the number
#' of edges that where identified by this function call and
#' \code{identifiedEdges[i,1] -> identifiedEdges[i,2]} was the ith edge
#' identified}
#' \item{\code{unsolvedParents}}{as the input argument but updated with
#' any newly identified edges}
#' \item{\code{solvedParents}}{as the input argument but updated with
#' any newly identified edges}
#' \item{\code{identifier}}{as the input argument but updated with
#' any newly identified edges}
#' }
#'
#' @export
#'
#' @references
#' Foygel, R., Draisma, J., and Drton, M. (2012) Half-trek criterion for
#' generic identifiability of linear structural equation models.
#' \emph{Ann. Statist.} 40(3): 1682-1713
htcIdentifyStep <- function(mixedGraph, unsolvedParents, solvedParents, identifier) {
identifiedEdges <- numeric(0)
allNodes <- mixedGraph$nodes()
solvedNodes <- which(sapply(unsolvedParents, FUN = function(x) {
length(x) == 0
}))
for (i in setdiff(allNodes, solvedNodes)) {
htrFromNode <- mixedGraph$htrFrom(i)
allowedNodes <- setdiff(solvedNodes, mixedGraph$siblings(i))
allowedNodes <- union(setdiff(allNodes, htrFromNode), allowedNodes)
nodeParents <- mixedGraph$parents(i)
if (length(allowedNodes) < length(nodeParents)) {
next
}
halfTrekSystemResult <- mixedGraph$getHalfTrekSystem(allowedNodes, nodeParents)
if (halfTrekSystemResult$systemExists) {
identifiedEdges <- c(identifiedEdges, as.integer(rbind(nodeParents, i)))
activeFrom <- halfTrekSystemResult$activeFrom
identifier <- createHtcIdentifier(identifier, activeFrom, nodeParents,
i, intersect(activeFrom, htrFromNode))
solvedParents[[i]] <- nodeParents
unsolvedParents[[i]] <- integer(0)
solvedNodes <- c(i, solvedNodes)
}
}
return(list(identifiedEdges = matrix(identifiedEdges, byrow = T, ncol = 2), unsolvedParents = unsolvedParents,
solvedParents = solvedParents, identifier = identifier))
}
#' Determines which edges in a mixed graph are HTC-identifiable.
#'
#' Uses the half-trek criterion of Foygel, Draisma, and Drton (2012) determine
#' which edges in a mixed graph are generically identifiable. Depending on your
#' application it faster to use the \code{\link{graphID.htcID}} function
#' instead of this one, this function has the advantage of returning additional
#' information.
#'
#' @export
#'
#' @inheritParams generalGenericID
#' @inheritParams semID
#'
#' @return see the return value of \code{\link{generalGenericID}}.
#'
#' @references
#' Foygel, R., Draisma, J., and Drton, M. (2012) Half-trek criterion for
#' generic identifiability of linear structural equation models.
#' \emph{Ann. Statist.} 40(3): 1682-1713.
#'
#' Jin Tian. 2005. Identifying direct causal effects in linear models. In
#' \emph{Proceedings of the 20th national conference on Artificial intelligence
#' - Volume 1} (AAAI'05), Anthony Cohn (Ed.), Vol. 1. AAAI Press 346-352.
htcID <- function(mixedGraph, tianDecompose = T) {
result <- generalGenericID(mixedGraph, list(htcIdentifyStep), tianDecompose = tianDecompose)
result$call <- match.call()
return(result)
}
#' Determines if a mixed graph is HTC-identifiable.
#'
#' Uses the half-trek criterion of Foygel, Draisma, and Drton (2013) to check
#' if an input mixed graph is generically identifiable.
#'
#' @export
#'
#' @inheritParams graphID
#'
#' @return The vector of HTC-identifiable nodes.
#'
#' @references
#' Foygel, R., Draisma, J., and Drton, M. (2012) Half-trek criterion for
#' generic identifiability of linear structural equation models.
#' \emph{Ann. Statist.} 40(3): 1682-1713.
graphID.htcID <- function(L, O) {
m <- nrow(L)
validateMatrices(L, O)
O <- 1 * ((O + t(O)) != 0)
# 1 & 2 = source & target 2 + {1,...,m} = L(i) for i=1,...,m 2+m + {1,...,m} =
# R(i)-in for i=1,...,m 2+2*m + {1,...,m} = R(i)-out for i=1,...,m
Cap.matrix.init <- matrix(0, 2 + 3 * m, 2 + 3 * m)
for (i in 1:m) {
# edge from L(i) to R(i)-in, and to R(j)-in for all siblings j of i
Cap.matrix.init[2 + i, 2 + m + c(i, which(O[i, ] == 1))] <- 1
# edge from R(i)-in to R(i)-out
Cap.matrix.init[2 + m + i, 2 + 2 * m + i] <- 1
# edge from R(i)-out to R(j)-in for all directed edges i->j
Cap.matrix.init[2 + 2 * m + i, 2 + m + which(L[i, ] == 1)] <- 1
}
# when testing if a set A satisfies the HTC with respect to a node i, need to add
# (1) edge from source to L(j) for all j in A and (2) edge from R(j)-out to
# target for all parents j of i
Dependence.matrix <- O + diag(m)
for (i in 1:m) {
Dependence.matrix <- (Dependence.matrix + Dependence.matrix %*% L > 0)
}
Solved.nodes <- rep(0, m)
Solved.nodes[which(colSums(L) == 0)] <- 1 # nodes with no parents
change <- 1
count <- 1
while (change == 1) {
change <- 0
for (i in which(Solved.nodes == 0)) {
A <- setdiff(c(which(Solved.nodes > 0), which(Dependence.matrix[i, ] ==
0)), c(i, which(O[i, ] == 1)))
Cap.matrix <- Cap.matrix.init
Cap.matrix[1, 2 + A] <- 1
Cap.matrix[2 + 2 * m + which(L[, i] == 1), 2] <- 1
flow <- graph.maxflow(graph.adjacency(Cap.matrix), source = 1, target = 2)$value
if (flow == sum(L[, i])) {
change <- 1
count <- count + 1
Solved.nodes[i] <- count
}
}
}
if (all(Solved.nodes == 0)) {
Solved.nodes <- NULL
} else {
Solved.nodes <- order(Solved.nodes)[(1 + sum(Solved.nodes == 0)):m]
}
return(Solved.nodes)
}
#' Check for generic infinite-to-one via the half-trek criterion.
#'
#' Checks if a mixed graph is infinite-to-one using the half-trek criterion
#' presented by Foygel, Draisma, and Drton (2012).
#'
#' @export
#'
#' @inheritParams graphID
#'
#' @return TRUE if the graph could be determined to be generically
#' non-identifiable, FALSE if this test was inconclusive.
#'
#' @references
#' Foygel, R., Draisma, J., and Drton, M. (2012) Half-trek criterion for
#' generic identifiability of linear structural equation models.
#' \emph{Ann. Statist.} 40(3): 1682-1713.
graphID.nonHtcID <- function(L, O) {
m <- nrow(L)
validateMatrices(L, O)
if (m == 1) {
# Always identifiable in this case
return(F)
}
O <- 1 * ((O + t(O)) != 0)
# 1 & 2 = source & target 2 + {1,...,N} = L{i_n,j_n} for the n-th nonsibling
# pair, n=1,...,N (2+N) + 2*m^2 = 2 copies of R_i(j) -- in copy & out copy where
# (2+N) + (i-1)*m + j = R_i(j) in copy & (2+N+m^2) + (i-1)*m + j = R_i(j) out
# copy
nonsibs <- NULL
N <- 0
for (i in 1:(m - 1)) {
for (j in (i + 1):m) {
if (O[i, j] == 0) {
N <- N + 1
nonsibs <- rbind(nonsibs, c(i, j))
}
}
}
Cap.matrix <- matrix(0, 2 * m^2 + N + 2, 2 * m^2 + N + 2)
if (N != 0) {
Cap.matrix[1, 2 + (1:N)] <- 1 # edges from source to L{i,j} for each
# nonsibling pair {i,j} = nonsibs[n,1:2] edge from L{i,j} to R_i(j), and to
# R_i(k)-in for all siblings k of node j
for (n in 1:N) {
Cap.matrix[2 + n, 2 + N + (nonsibs[n, 1] - 1) * m + c(nonsibs[n, 2],
which(O[nonsibs[n, 2], ] == 1))] <- 1
# edge from L{i,j} to R_j(i), and to R_j(i)-in for all siblings k of node i
Cap.matrix[2 + n, 2 + N + (nonsibs[n, 2] - 1) * m + c(nonsibs[n, 1],
which(O[nonsibs[n, 1], ] == 1))] <- 1
}
}
for (i in 1:m) {
# edge from R_i(j)-out to target when j is a parent of i
Cap.matrix[2 + N + m^2 + (i - 1) * m + which(L[, i] == 1), 2] <- 1
for (j in 1:m) {
# edge from R_i(j)-in to R_i(j)-out
Cap.matrix[2 + N + (i - 1) * m + j, 2 + N + m^2 + (i - 1) * m + j] <- 1
# edge from R_i(j)-out to R_i(k)-in where j->k is a directed edge
Cap.matrix[2 + N + m^2 + (i - 1) * m + j, 2 + N + (i - 1) * m + which(L[j,
] == 1)] <- 1
}
}
HTC.nonID <- graph.maxflow(igraph::graph.adjacency(Cap.matrix), source = 1, target = 2)$value
return(HTC.nonID < sum(L))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.