# R/gacFuns.R In pcalg: Methods for Graphical Models and Causal Inference

#### Documented in gac

```cond2node <- function(s, x,y,z,m,type) {
## INPUT: s in format for newStackEls
## OUTPUT: TRUE if pth is found; o/w FALSE

## Is there a path of length one ?
tmp <- newStackEls(s, x,y,z,m,type)
if (tmp\$suc) {
suc <- TRUE
} else {
## Is there a longer path ?
suc <- FALSE
if (length(tmp\$res) > 0) {
## nbrs to consider
i <- 0
while ( (!suc) & (i < length(tmp\$res)) ) {
i <- i+1
s2 <- tmp\$res[[i]]
suc <- cond2node(s2, x,y,z,m,type)
}
}
}
suc
}
cond2 <- function(x,y,z,m,type) {
## OUTPUT: True if cond 2 is true; o/w false
lx <- length(x)
pthFound <- rep(FALSE, lx)
for (i in 1:lx) {
s <- list(pth=x[i], ncp=FALSE)
## cond2node: TRUE if path is found
pthFound[i] <- cond2node(s,x,y,z,m,type)
}
## if no path is found, cond 2 is true
all(pthFound == FALSE)
}

defStat <- function(ln, cn, nn, m, type = "pag") {
## INPUT: Node positions ln, cn, nn; adj.mat m (for DAG/CPDAG/MAG/PAG)
## OUTPUT: TRUE if path ln-cn-nn in m is def.stat. o/w FALSE
res <- FALSE
stopifnot(type %in% c("dag", "cpdag", "mag", "pag"))

if (type %in% c("mag", "pag")) {
isCollider <- (m[ln,cn] == 2) & (m[nn,cn] == 2)
hasTail <- (m[ln,cn] == 3) | (m[nn,cn] == 3)
hasCircles <- (m[ln,cn] == 1) & (m[nn,cn] == 1)
isUnshielded <- (m[ln,nn] == 0) & (m[nn,ln] == 0)
isDefNonCollider <- ( (hasCircles & isUnshielded) | hasTail )
} else {
## DAG, CPDAG
isCollider <- (m[ln,cn]==0 & m[cn,ln]==1) & (m[cn,nn]==1 & m[nn,cn]==0)
hasTail <- (m[ln,cn]==1 & m[cn,ln]==0) | (m[cn,nn]==0 & m[nn,cn]==1)
isUndirected <- (m[ln,cn]==1 & m[cn,ln]==1) | (m[cn,nn]==1 & m[nn,cn]==1)
isUnshielded <- (m[ln,nn]==0 & m[nn,ln]==0)
isDefNonCollider <- ( (isUndirected & isUnshielded) | hasTail )
}

if (isCollider) {
res <- TRUE
} else {
if (isDefNonCollider) res <- TRUE
}

res
}
desc <- function(m, possible = FALSE, type = type) {
## computes all (possible) descendants of every node
## INPUT: adj.matrix m in MAG/PAG or DAG/CPDAG coding; if possible=TRUE, possible desc are found; o/w descendants
## OUPUT: list of vectors; vector at list position i contains node positions of (possible) descendants of node i
p <- ncol(m)
res <- vector(mode = "list", length = p)
for (i in 1:p) {
## res[[i]] <- possibleDeProper(m=m, x=i, y=NULL, possible=possible)
res[[i]] <- possDe(m = m, x = i, y = NULL,
possible = possible,
ds = FALSE, type = type)
}
res
}
forbiddenNodes <- function(m,x,y,type)
{
## INPUT: adj.matrix in DAG,CPDAG,MAG,PAG coding; sets of node positions x and y
## OUTPUT: set of node positions of nodes in the forbidden set (sorted)
n1 <- length(x)
n2 <- length(y)

possDeX <- possAnY <- c()
for(i in 1:max(n1,n2))
{
if (i <= n1) {
#find all possible descendants of a node x[i] that are on a proper path
#relative to x (exclude x[i] because descendants of x[i] that are not
#on a proper path are allowed)
## possDeX <- union(possDeX, setdiff(possibleDeProper(m,x[i],x),x[i]))
pdpTmp <- possDe(m = m, x = x[i], y = x, possible = TRUE,
ds = FALSE, type = type)
possDeX <- union(possDeX, setdiff(pdpTmp, x[i]))
}
if (i <= n2) {
#find all possible ancestors of a node y[i] that are
#on a proper path relative to x
## possAnY <- union(possAnY, possibleAnProper(m,y[i],x))
papTmp <- possAn(m = m, x = y[i], y = x, possible = TRUE,
ds = FALSE, type = type)
possAnY <- union(possAnY, papTmp)
}

}

#a set of all nodes on a proper possibly directed path from X to Y
pdp <- intersect(possDeX,possAnY)

#the forbiden node set are all possible descendants of nodes in pdp
fbnodes <- c()
if (length(pdp) > 0) {
for(j in 1:length(pdp))
{
## fbnodes <- union(fbnodes,possibleDeProper(m,pdp[j],c()))
pdpTmp2 <- possDe(m = m, x = pdp[j], y = c(), possible = TRUE,
ds = FALSE, type = type)
fbnodes <- union(fbnodes, pdpTmp2)
}
}

if (length(fbnodes) > 0) {
return(sort(fbnodes))
} else {
return(fbnodes)
}

}

##gac for pdags
##
## this function is the same as its equivalent in pcalg except it calls
## the pdag fucntions given above
##
##currrently separated dags,pdags,cpdags from mags,pags but could be merged

##however the res[3] must be calculated separately for dags,pdags,cpdags and mags,pags

## uses amat.cpdag encoding amat[i,j]=0,amat[j,i]=1 <=> i ->j
gac <- function (amat, x, y, z, type = "pag")
{
if (inherits(amat, "amat")) amat <- as(amat, "matrix")
if (!is.null(dimnames(amat))) dimnames(amat) <- NULL
if (type %in% c("dag", "cpdag", "pdag")) {
if (!isValidGraph(amat = amat, type = type)) {
message("The input graph is not a valid ",type,
". See function isValidGraph() for details.")
}
}
res <- c(isAmenable(amat, x = x, y = y, type = type),
NA, NA)
if (type %in% c("dag","pdag","cpdag")){
f <- bforbiddenNodes(m = amat, x = x, y = y)
res[2] <- (length(intersect(f, z)) == 0)
res[3] <- cond3fast(x = x, y = y, z = z, m = amat)
} else { ##the code for MAGs, PAGs is still the same
f <- forbiddenNodes(amat, x = x, y = y)
res[2] <- (length(intersect(f, z)) == 0)
res[3] <- cond2(x = x, y = y, z = z, m = amat, type = type)
}
list(gac = all(res), res = res, f = f)
}
isAmenable <- function(m,x,y, type = "pag") {
## INPUT: adj.matrix m; sets of node positions x and y; type in DAG, CPDAG,
## MAG or PAG
## OUTPUT: TRUE if m is amenabel wrt x,y; o/w FALSE
found <- FALSE ## if found == TRUE at any time, graph is not amenable wrt x,y
## DAG is always amenable
if (type %in% c("pdag", "cpdag","dag","pag","mag")) { ##changed added dag just in case, makes no difference
if (type == "dag")  ## added the if case for dags
return(!found)
i <- 0
p <- length(x)

## for all nodes in x, if amenability is still possible
while ( (i<p) & !found) {
i <- i+1
## posDesc of x[i] without going through any other x node
## posDesc <- possibleDeProper(m,x[i],x[-i])
posDesc <- possDe(m = m, x = x[i], y = x[-i], possible = TRUE,
ds = FALSE, type = type)
## potential problem for amenability only if there is a
## pdp from x[i] to y
if ( length(intersect(y, posDesc)) != 0 ) {
nb <- as.vector(which(m[x[i],]!=0 | m[,x[i]]!=0)) ## nbrs of x[i]
## potentially first node on pdp from x[i] to y; however, not yet sure
cand <- intersect(nb, posDesc)
j <- 0
## for all candidate nodes, if amenability is still possible
## (also covers case if cand is empty)
while ( (j<length(cand)) & !found ) {
j <- j+1
## check if there is a pdp from cand[j] to y without going through x[i]
## cand could already be in y
## pathOK <- ( length(intersect(y, possibleDeProper(m,cand[j],x[i]))) != 0 )
pdpTemp <- possDe(m = m, x = cand[j], y = x[i],
possible = TRUE, ds = FALSE,
type = type)
pathOK <- ( length(intersect(y, pdpTemp)) != 0 )
if (pathOK) {
isPDAG <- (type == "pdag" | type == "cpdag") ##changed
PDAGproblem <- (isPDAG & (m[x[i],cand[j]] == 1)) ##changed
PAGproblem1 <- (!isPDAG & (m[x[i], cand[j]] != 2) & (m[cand[j], x[i]] != 3)) ##changed
isDirEdge <- ((m[x[i], cand[j]] == 2) & (m[cand[j], x[i]] == 3))
PAGproblem2 <- (!isPDAG & isDirEdge & !visibleEdge(m, x[i], cand[j])) ##changed
found <- (PDAGproblem | PAGproblem1 | PAGproblem2)
} ## if pathOK
} ## while cand
} ## if path from x[i] to y
} ## while x
return(!found)
} else { ## if graph type not known
cat("Not a valid graph type! Should be written in lowercase! \n")
return(NULL)
}
}
mcon <- function(ln,cn,nn,m,z,descList) {
## INPUT: node positions of last, current and next node on path; adj.mat. m
## in DAG/CPDAGMAG/PAG format;
## set z; descList as returned from desc(m)
## OUTPUT: TRUE if ln and nn are mcon given z on path ln-cn-nn
res <- FALSE
isColliderPAG <- ( (m[ln,cn] == 2) & (m[nn,cn]==2) )
isColliderDAG <- ( (m[ln,cn]==0 & m[cn,ln]==1) & (m[cn,nn]==1 & m[nn,cn]==0) )
isCollider <- ( isColliderPAG | isColliderDAG)
if (isCollider) {
if (any(descList[[cn]] %in% z)) res <- TRUE
} else {
if ( !(cn %in% z) ) res <- TRUE
}
res
}
ncEdge <- function(cn, nn, m, type = "pag") {
stopifnot(type %in% c("dag", "cpdag", "mag", "pag"))
res <- FALSE
if (type %in% c("mag", "pag")) {
if (m[nn,cn]==2) res <- TRUE
} else { ## must be DAG or CPDAG
if ( (m[cn,nn] == 1) & (m[nn,cn] == 0)) res <- TRUE
}
res
}
newStackEls <- function(s,x,y,z,m,type) {
## INPUT: previous stack element of type list: {pth (num. vec.), ncp (bool);
## ncp is TRUE, if pth is non-causal; o/w FALSE
## OUTPUT: List 'res' of new stack elements (poss. empty); suc = TRUE, if an
## open path to y was found
suc <- FALSE
pth <- s\$pth; ncp <- s\$ncp
lp <- length(pth)
if (lp == 1) {
## Starting node on search
cn <- pth
nb <- setdiff(as.vector(which(m[cn,]!=0 | m[,cn]!=0)), x)
## nb does not contain x -> proper path
res <- vector("list", length(nb))
i <- 0
while ( (!suc) & (i < length(nb)) ) {
## Search through possible neighbors
i <- i + 1
nn <- nb[i]
ncpTmp <- ( ncp | ncEdge(cn,nn,m,type=type) )
## suc=T: Exists proper, m-con, d.s., non-c path from x to y
suc <- ( (nn %in% y) & ncpTmp )
res[[i]] <- list(pth = c(cn,nn), ncp = ncpTmp)
}
} else {
## Path has at least 2 nodes
descList <- desc(m, type = type)
ln <- pth[lp-1]; cn <- pth[lp]
nb <- setdiff(setdiff(as.vector(which(m[cn,]!=0 | m[,cn]!=0)), x), pth)
res <- list()
j <- 0
jj <- 0
while( (!suc) & (j<length(nb)) ) {
j <- j+1
nn <- nb[j]
mc <- mcon(ln=ln, cn=cn, nn=nn, m=m, z=z, descList=descList)
ds <- defStat(ln=ln, cn=cn, nn=nn, m=m, type = type)
ncpTmp <- ( ncp | ncEdge(cn,nn,m,type=type) )
if (mc & ds) {
jj <- jj + 1
## proper, m-connecting & def.stat. path
if ( (nn %in% y) & ncpTmp ) suc <- TRUE
res[[jj]] <- list(pth=c(pth,nn), ncp=ncpTmp)
}
}
}
list(res=res, suc=suc)
}

## ## CPDAG
## cp <- function(m) {
##     require(Rgraphviz)
##     plot(as(t(m),"graphNEL"))
## }

## ## MAG/PAG
## mp <- function(m) {
##     colnames(m) <- rownames(m) <- as.character(1:ncol(m))
##     require(Rgraphviz)
##     plotAG(m)
## }
possibleAnProper <- function(m,x,y=NULL)
{
## INPUT
## m: Adjacency matrix (coding for MAG/PAG)
## x: Starting node (col.position in m)
## y: Set Y (col.positions in m)
## OUTPUT
## All nodes with a possibly directed path from a to x not going through y
## including x itself

p <- length(m[,1]) ## nmb of nodes
q <- v <- rep(0,p) ## queue
## q has col.pos. of unvisited nodes = queue
## v has col.pos. of visited nodes
i <- k <-  1 ## i: end of queue; k: current point in queue
q[1] <- x ## x is first node in queue
tmp <- m ## ???

while(q[k]!=0 & k<=i) {## queue is not empty & current pos is within queue
t <- q[k] ## take new node from queue
v[k] <- t ## mark t as visited (if t is in v, it is visited)
k <- k+1 ## increase current position in queue
for(j in 1:p) { ## check if j is a possible parent of t
## check if edgemark at j is circle (1) or tail (3) but not head (2) or empty (0)
if (tmp[t,j] %in% c(1,3)) { ## works for CPDAG, too
## aren't already scheduled for a visit
## are not in y
if (!(j %in% q) & !(j %in% y)) {##
i <- i+1 ## incearse size of queue by one
q[i] <- j ## add node to queue
}
}
}
}
sort(setdiff(v,c(0))) ## remove trailing zeros from initial vector
}
#finds all possible descendants of a node x in a graph
#that are on a proper path relative to nodeset Y (that is, that don't go through Y)
possibleDeProper <- function(m,x,y=NULL,possible = TRUE)
{
## INPUT: adj.mat. m in MAG/PAG or DAG/CPDAG coding; node pos x;
## set of node pos y;
## If possible == TRUE, possible Desc. are found; o/w descendents are found
## OUTPUT: Node positions of (possible) descendents (ignoring path trough y); sorted
#q denotes unvisited nodes/ nodes in queue
#v denotes visited nodes
q <- v <- rep(0,length(m[,1]))
i <- k <-  1
q[i] <- x

while(q[k]!=0 & k<=i)
{
t <- q[k]
#mark t as visited
v[k] <- t
k <- k+1
#in this for cycle
#add all nodes that have a possibly directed
#edge with node t and all parents of node t to queue
for(j in 1:length(m[1,])) {
if (possible) {
## find possible descendents
collect <- ( ((m[t,j]!=0) && (m[j,t] %in% c(1,3))) || ## PAG; CPDAG t --- j
((m[t,j]==0) && (m[j,t]==1)) ) ## CPDAG: t ---> j
} else {
## find descendents
collect <- ( ((m[j,t]==3) && (m[t,j]==2)) ||
((m[j,t]==1) && (m[t,j]==0)) )
}
if (collect) {
#and that are on a proper path
if (!(j %in% q) & !(j %in% y))
{
i <- i+1
q[i] <- j
}
}
}
}
sort(setdiff(v,c(0)))
}

##The following function forms a proper back-door graph G_{XY}^{pbd}
##given the adjacency matrix, node position vectors x and y
##
gbg <-  function(m,x,y)
{
tmp <- m
for (i in seq_along(x)) {
Desc <- bPossibleDeProper(m,x[i],x[-i])
if (length(intersect(y, Desc)) != 0) {
ch <- as.vector(which(m[x[i], ] == 0 & m[, x[i]] == 1))
cand <- intersect(ch, Desc)
for(j in seq_along(cand)) {
##pathOK <- (length(intersect(y, bpossibleDeProper(m,cand[j], x[i]))) != 0)
## what if x - y? MISSING CASE!
pathOK <- (length(intersect(y, bPossibleDeProper(m,cand[j], x))) != 0) ||
(cand[j] %in% y) # added x - y case
if (pathOK) {
tmp[cand[j],x[i]] <- 0
}
}
}
}
tmp
}

## the following function tests the separation condition for gac-pdags
## it is a pseudo check for this condition, because provided
## that amenability and forbidden set are satisfied
## this will exactly check the separation condition

## however if amenability or forb are not satisfied
## this check is not necessarily accurate

## given a pdag adjacency matrix m it checks whether node position vector z
## satisfies the separation condition relative to (x,y)

## m encoding amat.cpdag as above m[i,j]=0, m[j,i]=1 <=> i -> j
cond3fast <- function(x, y, z, m)
{
## the idea is to use the lemma which says
## that if m is amenable and z satisfies the forbidden cond
## then z satisfies separation in the pdag iff
## z satisfies separation in at least one dag

## First find one dag D in the equivalence class represented by m
oneDag <- pdag2dag(as(t(m),"graphNEL"))
dagAmat <- t(as(oneDag\$graph,"matrix"))

## We instead of checking the separation crit in m
## w can check the d-separation criterion in D_{XY}^{pbd}
gb <- gbg(dagAmat,x,y)
msep(t(gb),x,y,z)  ## uses the transposed of amat.cpdag encoding
}

## note that provided that possibleDeProper is implemented in
## pcalg as looking at only def. stat. paths
## then the already implemented forbiddenNodes can be used directly
## (with the addition of type: pdag)

## m - adjacency matrix amat.cpdag encoding m[i,j]=0, m[j,i] = 1 <=> i -> j
## x - vector of node position
## y - vector of node positions

## as it is now this function can only be used for pdags, dags, cpdags
## because the bPossibleDeProper etc. functions are only defined for these
bforbiddenNodes <- function (m, x, y)
{
n1 <- length(x)
n2 <- length(y)
possDeX <- possAnY <- c()
for (i in 1:max(n1, n2)) {
if (i <= n1)
possDeX <- union(possDeX, setdiff(bPossibleDeProper(m, x[i], x), x[i])) ##changed
if (i <= n2)
possAnY <- union(possAnY, bPossibleAnProper(m, y[i], x)) ##changed
}
pdp <- intersect(possDeX, possAnY)
fbnodes <- c()
if (length(pdp) > 0) {
for (j in 1:length(pdp)) {
fbnodes <- union(fbnodes, bPossibleDeProper(m, pdp[j], c())) ##changed
}
}
if (length(fbnodes) > 0) {
return(sort(fbnodes))
}
else {
return(fbnodes)
}
}
```

## Try the pcalg package in your browser

Any scripts or data that you put into this service are public.

pcalg documentation built on June 5, 2018, 1:05 a.m.