udag2wanpdag {kpcalg}

Last kPC Algorithm Step: Extend Object with Skeleton to Completed PDAG

Description

This function performs the last (generalised transitive) step in the kpc algorithm. It transforms an object of the class "pcAlgo" containing a skeleton and corresponding conditional independence information into a weakly additive noise directed acyclic graph (CPDAG). The functions first determine the v-structures in the collider step, and then performs the Generalised Transitive Step as described in Tillman et al (2009) to orient as many of the remaining edges as possible.

Usage

udag2wanpdag(gInput, suffStat, indepTest = kernelCItest, alpha = 0.2,
  verbose = FALSE, unfVect = NULL, solve.confl = FALSE,
  orientCollider = TRUE, rules = rep(TRUE, 3))

Arguments

gInput

"pcAlgo"-object containing skeleton and conditional indepedence information.

suffStat

a list of sufficient statistics, containing all necessary elements for the conditional independence decisions in the function indepTest.

indepTest

A function for testing conditional independence. It is internally called as indepTest(x,y,S,suffStat). Default is kernelCItest.

alpha

significance level (number in (0,1) for the individual conditional independence tests.

verbose

0: No output; 1: Details

unfVect

vector containing numbers that encode ambiguous triples (as returned by pc.cons.intern). This is needed in the conservative and majority rule PC algorithms.

solve.confl

if TRUE, the orientation of the v-structures and the orientation rules work with lists for candidate sets and allow bi-directed edges to resolve conflicting edge orientations. Note that therefore the resulting object is order-independent but might not be a PDAG because bi-directed edges can be present.

orientCollider

if TRUE, collider are oriented.

rules

Array of length 3 containing TRUE or FALSE for each rule. TRUE in position i means that rule i (Ri) will be applied. By default, all rules are used.gInput

Details

First we perform a collider step, that is orienting triples a-b-c as a->b<-c iff b is not in separating set of a and c. Then we orient edges a-S as a->S if b_r is independent of a set S, where b_r are the residuals of b non parametrically regressed on S and parents of b and none of the edges S_i-a can be oriented as S_i->a, that is residuals S_i_r would be independent of a.

Value

An oriented object of class "pcAlgo".

References

Tillman, R. E., Gretton, A. and Spirtes, P. (2009). Nonlinear directed acyclic structure learning with weakly additive noise model. NIPS 22, Vancouver.

Examples

## Not run: 
library(pcalg)
set.seed(4)
n <- 300
data <- NULL
x1 <- 2*(runif(n)-0.5)
x2 <- x1 + runif(n)-0.5
x3 <- x1^2 + 0.6*runif(n)
x4 <- rnorm(n)
x5 <- x3 + x4^2 + 2*runif(n)
x6 <- 10*(runif(n)-0.5)
x7 <- x6^2 + 10*runif(n)
x8 <- 2*x7^2 + rnorm(n)
x9 <- x7 + 5*runif(n)
data <- cbind(x1,x2,x3,x4,x5,x6,x7,x8,x9)
true <- matrix(0,9,9)
true[c(1),c(2,3)]<-true[c(3,4),5]<-true[c(6),c(7)]<-true[c(7),c(8)]<-true[7,9]<-1
## estimate skeleton
resU1 <- skeleton(suffStat = list(data=data, ic.method="dcc.perm", p=200),
                  indepTest = kernelCItest,
                  verbose = TRUE, alpha = 0.1, p=9)

resU2 <- skeleton(suffStat = list(data=data, ic.method="hsic.gamma",
                             sig=1, numCol = 50),
                  indepTest = kernelCItest,
                  verbose = TRUE, alpha = 0.1, p=9)

resU3 <- skeleton(suffStat = list(data=data, ic.method="hsic.perm",
                             sig=1, numCol = 50, p=200),
                  indepTest = kernelCItest,
                  verbose = TRUE, alpha = 0.1, p=9)

resU4 <- skeleton(suffStat = list(data=data, ic.method="hsic.clust",
                             p=200, sig=1, numCluster=100, numCol = 50,
                             eps = 0.1, paral = 1),
                  indepTest = kernelCItest,
                  verbose = TRUE, alpha = 0.1, p=9)

resU5 <- skeleton(suffStat = list(C = cor(data), n = n),
                  indepTest = gaussCItest,
                  verbose = TRUE, alpha = 0.1, p=9)

if (require(Rgraphviz)) {
 par(mfrow=c(2,3))
 plot(resU1,main="dpc")
 plot(resU2,main="kpc-resid-gamma")
 plot(resU3,main="kpc-resid-perm")
 plot(resU4,main="kpc-clust")
 plot(resU5,main="pc")
 plot(as(true,"graphNEL"),main="True DAG")
}

## orient edges using three different methods
resD1 <- udag2wanpdag(gInput = resU1,
                      suffStat = list(data=data, ic.method="dcc.perm", sig=1, numCol = 50, p=200),
                      indepTest = kernelCItest,
                      verbose = TRUE, alpha = 0.1)
resD2 <- udag2wanpdag(gInput = resU1,
                      suffStat = list(data=data, ic.method="hsic.gamma", sig=1, numCol = 50),
                      indepTest = kernelCItest,
                      verbose = TRUE, alpha = 0.1)
resD3 <- udag2wanpdag(gInput = resU1,
                      suffStat = list(data=data, ic.method="hsic.perm", sig=1, numCol = 50, p=200),
                      indepTest = kernelCItest,
                      verbose = TRUE, alpha = 0.1)
resD4 <- udag2pdagRelaxed(gInput = resU1, verbose = T)
if (require(Rgraphviz)) {
 par(mfrow=c(2,3))
 plot(resD1,main="dpc")
 plot(resD2,main="kpc-resid-gamma")
 plot(resD3,main="kpc-resid-perm")
 plot(resD4,main="pc")
 plot(as(true,"graphNEL"),main="True DAG")
}

## End(Not run)