R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.1  (2005-06-20), ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for a HTML browser interface to help.
Type 'q()' to quit R.

> ### * <HEADER>
> ###
> attach(NULL, name = "CheckExEnv")
> assign(".CheckExEnv", as.environment(2), pos = length(search())) # base
> ## add some hooks to label plot pages for base and grid graphics
> setHook("plot.new", ".newplot.hook")
> setHook("persp", ".newplot.hook")
> setHook("grid.newpage", ".gridplot.hook")
> 
> assign("cleanEx",
+        function(env = .GlobalEnv) {
+ 	   rm(list = ls(envir = env, all.names = TRUE), envir = env)
+            RNGkind("default", "default")
+ 	   set.seed(1)
+    	   options(warn = 1)
+ 	   delayedAssign("T", stop("T used instead of TRUE"),
+ 		  assign.env = .CheckExEnv)
+ 	   delayedAssign("F", stop("F used instead of FALSE"),
+ 		  assign.env = .CheckExEnv)
+ 	   sch <- search()
+ 	   newitems <- sch[! sch %in% .oldSearch]
+ 	   for(item in rev(newitems))
+                eval(substitute(detach(item), list(item=item)))
+ 	   missitems <- .oldSearch[! .oldSearch %in% sch]
+ 	   if(length(missitems))
+ 	       warning("items ", paste(missitems, collapse=", "),
+ 		       " have been removed from the search path")
+        },
+        env = .CheckExEnv)
> assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now
> assign("ptime", proc.time(), env = .CheckExEnv)
> grDevices::postscript("ggm-Examples.ps")
> assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv)
> options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly"))
> options(warn = 1)    
> library('ggm')
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "DAG"
> 
> ### * DAG
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: DAG
> ### Title: Directed acyclic graphs (DAGs)
> ### Aliases: DAG
> ### Keywords: graphs models multivariate
> 
> ### ** Examples
> 
> ## A Markov chain
> DAG(y ~ x, x ~ z, z ~ u)
  y x z u
y 0 0 0 0
x 1 0 0 0
z 0 1 0 0
u 0 0 1 0
> 
> ## Another DAG
> DAG(y ~ x + z + u, x ~ u, z ~ u)
  y x z u
y 0 0 0 0
x 1 0 0 0
z 1 0 0 0
u 1 1 1 0
> 
> ## A DAG with an isolated node
> DAG(v ~ v, y ~ x + z, z ~ w + u)
  v y x z w u
v 0 0 0 0 0 0
y 0 0 0 0 0 0
x 0 1 0 0 0 0
z 0 1 0 0 0 0
w 0 0 0 1 0 0
u 0 0 0 1 0 0
> 
> ## There can be repetitions
> DAG(y ~ x + u + v, y ~ z, u ~ v + z)
  y x u v z
y 0 0 0 0 0
x 1 0 0 0 0
u 1 0 0 0 0
v 1 0 1 0 0
z 1 0 1 0 0
> 
> ## Interactions are ignored
> DAG(y ~ x*z + z*v, x ~ z)
  y x z v
y 0 0 0 0
x 1 0 0 0
z 1 1 0 0
v 1 0 0 0
> 
> ## A cyclic graph returns an error!
> ## Not run: DAG(y ~ x, x ~ z, z ~ y)
> 
> ## The order can be changed
> DAG(y ~ z, y ~ x + u + v,  u ~ v + z)
  y z x u v
y 0 0 0 0 0
z 1 0 0 1 0
x 1 0 0 0 0
u 1 0 0 0 0
v 1 0 0 1 0
> 
> ## If you want to order the nodes (topological sort of the DAG)
> DAG(y ~ z, y ~ x + u + v,  u ~ v + z, order=TRUE)
  v x z u y
v 0 0 0 1 1
x 0 0 0 0 1
z 0 0 0 1 1
u 0 0 0 0 1
y 0 0 0 0 0
> 
> 
> 
> cleanEx(); ..nameEx <- "In"
> 
> ### * In
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: In
> ### Title: Indicator matrix
> ### Aliases: In
> ### Keywords: array algebra graphs multivariate
> 
> ### ** Examples
> 
> ## A simple way to find the overall induced coincentration graph
> ## The DAG on p. 198 of Cox & Wermuth (1996)
> amat <- DAG(y1 ~ y2 + y3, y3 ~ y5, y4 ~ y5)
> A <- edgeMatrix(amat)
> In(crossprod(A))
   y1 y2 y3 y5 y4
y1  1  1  1  0  0
y2  1  1  1  0  0
y3  1  1  1  1  0
y5  0  0  1  1  1
y4  0  0  0  1  1
> 
> 
> 
> cleanEx(); ..nameEx <- "InducedGraphs"
> 
> ### * InducedGraphs
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: InducedGraphs
> ### Title: Graphs induced by marginalization or conditioning
> ### Aliases: inducedCovGraph inducedConGraph inducedRegGraph
> ###   inducedChainGraph inducedDAG InducedGraphs
> ### Keywords: graphs models multivariate
> 
> ### ** Examples
> 
> ## Define a DAG
> dag <- DAG(a ~ x, c ~ b+d, d~ x)
> dag
  a x c b d
a 0 0 0 0 0
x 1 0 0 0 1
c 0 0 0 0 0
b 0 0 1 0 0
d 0 0 1 0 0
> ## Induced covariance graph of a, b, d given the empty set.
> inducedCovGraph(dag, sel=c("a", "b", "d"), cond=NULL)
  a b d
a 0 0 1
b 0 0 0
d 1 0 0
> 
> ## Induced concentration graph of a, b, c given x
> inducedConGraph(dag, sel=c("a", "b", "c"), cond="x")
  a b c
a 0 0 0
b 0 0 1
c 0 1 0
> 
> ## Overall covariance graph
> inducedCovGraph(dag)
  a x c b d
a 0 1 1 0 1
x 1 0 1 0 1
c 1 1 0 1 1
b 0 0 1 0 0
d 1 1 1 0 0
> 
> ## Overall concentration graph
> inducedConGraph(dag)
  a x c b d
a 0 1 0 0 0
x 1 0 0 0 1
c 0 0 0 1 1
b 0 0 1 0 1
d 0 1 1 1 0
> 
> ## Induced covariance graph of x, b, d given c, x.
> inducedCovGraph(dag, sel=c("a", "b", "d"), cond=c("c", "x"))
  a b d
a 0 0 0
b 0 0 1
d 0 1 0
> 
> ## Induced concentration graph of a, x, c given d, b.
> inducedConGraph(dag, sel=c("a", "x", "c"), cond=c("d", "b"))
  a x c
a 0 1 0
x 1 0 0
c 0 0 0
> 
> ## The DAG on p. 198 of Cox & Wermuth (1996)
> dag <- DAG(y1~ y2 + y3, y3 ~ y5, y4 ~ y5)
> 
> ## Cf. figure 8.7 p. 203 in Cox & Wermuth (1996)
> inducedCovGraph(dag, sel=c("y2", "y3", "y4", "y5"), cond="y1")
   y2 y3 y4 y5
y2  0  1  1  1
y3  1  0  1  1
y4  1  1  0  1
y5  1  1  1  0
> inducedCovGraph(dag, sel=c("y1", "y2", "y4", "y5"), cond="y3")
   y1 y2 y4 y5
y1  0  1  0  0
y2  1  0  0  0
y4  0  0  0  1
y5  0  0  1  0
> inducedCovGraph(dag, sel=c("y1", "y2", "y3", "y4"), cond="y5")
   y1 y2 y3 y4
y1  0  1  1  0
y2  1  0  0  0
y3  1  0  0  0
y4  0  0  0  0
> 
> ## Cf. figure 8.8 p. 203 in Cox & Wermuth (1996)
> inducedConGraph(dag, sel=c("y2", "y3", "y4", "y5"), cond="y1")
   y2 y3 y4 y5
y2  0  1  0  0
y3  1  0  0  1
y4  0  0  0  1
y5  0  1  1  0
> inducedConGraph(dag, sel=c("y1", "y2", "y4", "y5"), cond="y3")
   y1 y2 y4 y5
y1  0  1  0  0
y2  1  0  0  0
y4  0  0  0  1
y5  0  0  1  0
> inducedConGraph(dag, sel=c("y1", "y2", "y3", "y4"), cond="y5")
   y1 y2 y3 y4
y1  0  1  1  0
y2  1  0  1  0
y3  1  1  0  0
y4  0  0  0  0
> 
> ## Cf. figure 8.9 p. 204 in Cox & Wermuth (1996)
> inducedCovGraph(dag, sel=c("y2", "y3", "y4", "y5"), cond=NULL)
   y2 y3 y4 y5
y2  0  0  0  0
y3  0  0  1  1
y4  0  1  0  1
y5  0  1  1  0
> inducedCovGraph(dag, sel=c("y1", "y2", "y4", "y5"), cond=NULL)
   y1 y2 y4 y5
y1  0  1  1  1
y2  1  0  0  0
y4  1  0  0  1
y5  1  0  1  0
> inducedCovGraph(dag, sel=c("y1", "y2", "y3", "y4"), cond=NULL)
   y1 y2 y3 y4
y1  0  1  1  1
y2  1  0  0  0
y3  1  0  0  1
y4  1  0  1  0
> 
> ## Cf. figure 8.10 p. 204 in Cox & Wermuth (1996)
> inducedConGraph(dag, sel=c("y2", "y3", "y4", "y5"), cond=NULL)
   y2 y3 y4 y5
y2  0  0  0  0
y3  0  0  0  1
y4  0  0  0  1
y5  0  1  1  0
> inducedConGraph(dag, sel=c("y1", "y2", "y4", "y5"), cond=NULL)
   y1 y2 y4 y5
y1  0  1  0  1
y2  1  0  0  1
y4  0  0  0  1
y5  1  1  1  0
> inducedConGraph(dag, sel=c("y1", "y2", "y3", "y4"), cond=NULL)
   y1 y2 y3 y4
y1  0  1  1  0
y2  1  0  1  0
y3  1  1  0  1
y4  0  0  1  0
> 
> ## An induced regression graph
> dag2 = DAG(Y ~ X+U, W ~ Z+U)
> inducedRegGraph(dag2, sel="W",  cond=c("Y", "X", "Z"))
  W
Y 1
X 1
Z 1
> 
> ## An induced DAG
> inducedDAG(dag2, order=c("X","Y","Z","W"))
  X Y Z W
X 0 1 0 1
Y 0 0 0 1
Z 0 0 0 1
W 0 0 0 0
> 
> ## An induced multivariate regression graph
> inducedRegGraph(dag2, sel=c("Y", "W"), cond=c("X", "Z"))
  Y W
X 1 0
Z 0 1
> 
> ## An induced chain graph with LWF interpretation
> dag3 = DAG(X~W, W~Y, U~Y+Z)
> cc = list(c("W", "U"), c("X", "Y", "Z"))
> inducedChainGraph(dag3, cc=cc, type="LWF")
  W U X Y Z
W 0 1 1 1 0
U 1 0 0 1 1
X 0 0 0 0 0
Y 0 0 0 0 1
Z 0 0 0 1 0
> 
> ## ... with AMP interpretation
> inducedChainGraph(dag3, cc=cc, type="AMP")
  W U X Y Z
W 0 1 1 1 1
U 1 0 0 1 1
X 0 0 0 0 0
Y 0 0 0 0 1
Z 0 0 0 1 0
> 
> ## ... with multivariate regression interpretation
> cc= list(c("U"), c("Z", "Y"), c("X", "W"))
> inducedChainGraph(dag3, cc=cc, type="MRG")
  U Z Y X W
U 0 1 1 0 0
Z 0 0 1 0 0
Y 0 1 0 1 1
X 0 0 0 0 1
W 0 0 0 1 0
> 
> 
> 
> cleanEx(); ..nameEx <- "SimpleGraphOperations"
> 
> ### * SimpleGraphOperations
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Simple Graph Operations
> ### Title: Simple graph operations
> ### Aliases: bd ch pa vertices
> ### Keywords: graphs models multivariate
> 
> ### ** Examples
> 
> ## find boundary of a subset of nodes of a DAG
> G <- DAG(y ~ x+b+a, b~a, x~a)
> vertices(G)
[1] "y" "x" "b" "a"
> bd("b", G)
[1] "y" "a"
> bd(c("b", "x"), G)
[1] "y" "a"
> bd("x", G)
[1] "y" "a"
> bd(c("x","b"), G)
[1] "y" "a"
> ## find boundary of a subset of nodes of an UG
> G <- UG(~ y*x*z + z*h*v)
> bd("z", G)
[1] "y" "x" "h" "v"
> bd(c("y", "x"), G)
[1] "z"
> bd("v", G)
[1] "z" "h"
> bd(c("x","v"), G)
[1] "y" "z" "h"
> ## children of a subset of nodes of a DAG
> G <- DAG(y ~ x+b+a, b~a, x~a)
> ch("b", G)
[1] "y"
> ch(c("b", "x"), G)
[1] "y"
> ch("x", G)
[1] "y"
> ch(c("a","x"), G)
[1] "y" "b"
> ## parents of a subset of nodes of a DAG
> pa("b", G)
[1] "a"
> pa(c("b", "x"), G)
[1] "a"
> pa("x", G)
[1] "a"
> pa(c("x","b"), G)
[1] "a"
> 
> 
> 
> cleanEx(); ..nameEx <- "UG"
> 
> ### * UG
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: UG
> ### Title: Defining an undirected graph (UG)
> ### Aliases: UG
> ### Keywords: graphs models multivariate
> 
> ### ** Examples
> 
> ## X independent of Y given Z
> UG(~ X*Z + Y*Z)
  X Z Y
X 0 1 0
Z 1 0 1
Y 0 1 0
> 
> # The saturated model
> UG(~ X*Y*Z)
  X Y Z
X 0 1 1
Y 1 0 1
Z 1 1 0
> 
> ## The model without three-way interactions has the same graph
> UG(~ X*Y + Y*Z + Z*X)
  X Y Z
X 0 1 1
Y 1 0 1
Z 1 1 0
> UG(~ (X + Y + Z)^2)
  X Y Z
X 0 1 1
Y 1 0 1
Z 1 1 0
> 
> ## Butterfly model defined from the cliques
> UG(~ mec*vec*alg + alg*ana*sta)
    mec vec alg ana sta
mec   0   1   1   0   0
vec   1   0   1   0   0
alg   1   1   0   1   1
ana   0   0   1   0   1
sta   0   0   1   1   0
> 
> ## Some isolated nodes
> UG(~x*y*z + a + b) 
  x y z a b
x 0 1 1 0 0
y 1 0 1 0 0
z 1 1 0 0 0
a 0 0 0 0 0
b 0 0 0 0 0
> 
> 
> 
> cleanEx(); ..nameEx <- "adjMatrix"
> 
> ### * adjMatrix
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: adjMatrix
> ### Title: Adjacency matrix of a graph
> ### Aliases: adjMatrix
> ### Keywords: array algebra graphs multivariate
> 
> ### ** Examples
> 
> amat <- DAG(y ~ x+z, z~u+v)
> E <- edgeMatrix(amat)
> adjMatrix(E)
  y x z u v
y 0 0 0 0 0
x 1 0 0 0 0
z 1 0 0 0 0
u 0 0 1 0 0
v 0 0 1 0 0
> 
> 
> 
> cleanEx(); ..nameEx <- "allEdges"
> 
> ### * allEdges
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: allEdges
> ### Title: All edges of a graph
> ### Aliases: allEdges
> ### Keywords: graphs models multivariate
> 
> ### ** Examples
> 
> ## A UG graph
> allEdges(UG(~ y*v*k +v*k*d+y*d))
    
 1 2
 1 3
 1 4
 2 3
 2 4
 3 4
> 
> ## A DAG
> allEdges(DAG(u~h+o+p, h~o, o~p))
    
 2 1
 3 1
 3 2
 4 1
 4 3
> 
> 
> 
> cleanEx(); ..nameEx <- "basiSet"
> 
> ### * basiSet
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: basiSet
> ### Title: Basis set of a DAG
> ### Aliases: basiSet
> ### Keywords: graphs models multivariate
> 
> ### ** Examples
> 
> ## See Shipley (2000), Figure 2, p. 213
> A <- DAG(x5~ x3+x4, x3~ x2, x4~x2, x2~ x1)
> basiSet(A)
[[1]]
[1] "x1" "x4" "x2"

[[2]]
[1] "x1" "x3" "x2"

[[3]]
[1] "x1" "x5" "x4" "x3"

[[4]]
[1] "x2" "x5" "x1" "x4" "x3"

[[5]]
[1] "x4" "x3" "x2"

> 
> 
> 
> cleanEx(); ..nameEx <- "bfs"
> 
> ### * bfs
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: bfs
> ### Title: Breadth first search
> ### Aliases: bfs
> ### Keywords: graphs models multivariate
> 
> ### ** Examples
> 
> ## Finding a spanning tree of the butterfly graph
> bfs(UG(~ a*b*o + o*u*j))
$tree
  a b o u j
a 0 1 1 0 0
b 1 0 0 0 0
o 1 0 0 1 1
u 0 0 1 0 0
j 0 0 1 0 0

$branches
    
 1 2
 1 3
 3 4
 3 5

$chords
    
 2 3
 4 5

> ## Starting from another node
> bfs(UG(~ a*b*o + o*u*j), v=3)
$tree
  a b o u j
a 0 0 1 0 0
b 0 0 1 0 0
o 1 1 0 1 1
u 0 0 1 0 0
j 0 0 1 0 0

$branches
    
 3 1
 3 2
 3 4
 3 5

$chords
    
 1 2
 4 5

> 
> 
> 
> cleanEx(); ..nameEx <- "checkIdent"
> 
> ### * checkIdent
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: checkIdent
> ### Title: Identifiability of a model with one latent variable
> ### Aliases: checkIdent
> ### Keywords: graphs models multivariate
> 
> ### ** Examples
> 
> ## See DAG in Figure 4 (a) in Stanghellini & Wermuth (2003)
> d <- DAG(y1 ~ y3, y2 ~ y3 + y5, y3 ~ y4 + y5, y4 ~ y6)
> checkIdent(d, "y3")  # Identifiable
 T1.i T1.ii  T2.i T2.ii 
FALSE  TRUE  TRUE FALSE 
> checkIdent(d, "y4")  # Not identifiable?
 T1.i T1.ii  T2.i T2.ii 
FALSE FALSE FALSE FALSE 
> 
> ## See DAG in Figure 5 (a) in Stanghellini & Wermuth (2003)
> d <- DAG(y1 ~ y5+y4, y2 ~ y5+y4, y3 ~ y5+y4)
> checkIdent(d, "y4")  # Identifiable
 T1.i T1.ii  T2.i T2.ii 
FALSE FALSE FALSE  TRUE 
> checkIdent(d, "y5")  # Identifiable
 T1.i T1.ii  T2.i T2.ii 
FALSE FALSE FALSE  TRUE 
> 
> ## A simple function to check identifiability for each node
> 
> is.ident <- function(amat){
+ ### Check suff. conditions on each node of a DAG.
+    p <- nrow(amat)
+    ## Degrees of freedom
+      df <- p*(p+1)/2 - p  - sum(amat==1) - p + 1
+    if(df <= 0)
+        warning(paste("The degrees of freedom are ", df))
+     a <- rownames(amat)
+     for(i in a) {
+       b <- checkIdent(amat, latent=i)
+       if(TRUE %in% b)
+         cat("Node", i, names(b)[!is.na(b)], "\n")
+       else
+         cat("Unknown.\n")
+     }
+   }
> 
> 
> 
> cleanEx(); ..nameEx <- "cliques"
> 
> ### * cliques
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: cliques
> ### Title: Cliques of an undirected graph
> ### Aliases: cliques
> ### Keywords: graphs models multivariate
> 
> ### ** Examples
> 
> u <- UG(~ a*b*c + c*d*e*g)
> u
  a b c d e g
a 0 1 1 0 0 0
b 1 0 1 0 0 0
c 1 1 0 1 1 1
d 0 0 1 0 1 1
e 0 0 1 1 0 1
g 0 0 1 1 1 0
> cliques(u)
[[1]]
[1] "a" "b" "c"

[[2]]
[1] "c" "d" "e" "g"

> 
> graph22 <- UG(~x1*x2*x3+x3*x4*x5*x6+x5*x7*x8*x9*x10+x9*x11*x12*x13*x14+
+ x12*x15+x15*x16*x17*x18+x17*x1*x18+x1*x5*x7*x19*x20+x6*x20*x21+x22)
> 
> cliques(graph22)
[[1]]
[1] "x1" "x2" "x3"

[[2]]
[1] "x1" "x3" "x5"

[[3]]
[1] "x1"  "x5"  "x7"  "x19" "x20"

[[4]]
[1] "x1"  "x17" "x18"

[[5]]
[1] "x3" "x4" "x5" "x6"

[[6]]
[1] "x5"  "x7"  "x8"  "x9"  "x10"

[[7]]
[1] "x5"  "x6"  "x20"

[[8]]
[1] "x6"  "x20" "x21"

[[9]]
[1] "x9"  "x11" "x12" "x13" "x14"

[[10]]
[1] "x12" "x15"

[[11]]
[1] "x15" "x16" "x17" "x18"

[[12]]
[1] "x22"

> 
> 
> 
> cleanEx(); ..nameEx <- "cmpGraph"
> 
> ### * cmpGraph
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: cmpGraph
> ### Title: The complementary graph
> ### Aliases: cmpGraph
> ### Keywords: graphs models multivariate
> 
> ### ** Examples
> 
> ## A chordless four-cycle
> four <- UG(~ a*b + b*d + d*e + e*a)
> four
  a b d e
a 0 1 0 1
b 1 0 1 0
d 0 1 0 1
e 1 0 1 0
> cmpGraph(four)
  a b d e
a 0 0 1 0
b 0 0 0 1
d 1 0 0 0
e 0 1 0 0
> 
> 
> 
> cleanEx(); ..nameEx <- "conComp"
> 
> ### * conComp
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: conComp
> ### Title: Connectivity components
> ### Aliases: conComp
> ### Keywords: graphs models multivariate
> 
> ### ** Examples
> 
> ## three connected components
> conComp(UG(~a*c+c*d+e+g*o*u))
[1] 1 1 1 2 3 3 3
> ## a connected graph
> conComp(UG(~ a*b+b*c+c*d+d*a))
[1] 1 1 1 1
> 
> 
> 
> cleanEx(); ..nameEx <- "correlations"
> 
> ### * correlations
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: correlations
> ### Title: Marginal and partial correlations
> ### Aliases: correlations
> ### Keywords: array graphs models multivariate
> 
> ### ** Examples
> 
> ## See Table 6.1 in Cox & Wermuth (1996)
> data(glucose)
> correlations(glucose)
            Y          X           Z           U            V            W
Y  1.00000000 -0.2444579  0.09354755 -0.12798823  0.087942794 -0.249825305
X -0.34075247  1.0000000 -0.33532889 -0.07963908  0.039741792  0.003602872
Z  0.14979043 -0.4911540  1.00000000  0.40922980 -0.258331690  0.275120429
U  0.02734320 -0.3247271  0.52054597  1.00000000 -0.085859607 -0.076665692
V  0.04199025  0.1417811 -0.33278493 -0.22852134  1.000000000  0.158985954
W -0.12238307 -0.1102499  0.27522852  0.09900332  0.047259359  1.000000000
A -0.31889444  0.3333154 -0.25667957 -0.19593153 -0.005196467 -0.250207430
B -0.06634013  0.0861791  0.07599028 -0.06409980 -0.221254867  0.071501311
             A           B
Y -0.288992196 -0.05320535
X  0.174071199  0.13665929
Z -0.002712233  0.08902947
U -0.119850983 -0.14027029
V -0.068351329 -0.22966551
W -0.248252892  0.03783576
A  1.000000000 -0.13733490
B -0.085404874  1.00000000
> 
> 
> 
> cleanEx(); ..nameEx <- "cycleMatrix"
> 
> ### * cycleMatrix
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: cycleMatrix
> ### Title: Fundamental cycles
> ### Aliases: cycleMatrix
> ### Keywords: graphs models multivariate
> 
> ### ** Examples
> 
> ## Three cycles
> cycleMatrix(UG(~a*b*d+d*e+e*a*f))
  1 2 1 3 1 4 1 5 2 3 3 4 4 5
1   1   1   0   0   1   0   0
2   0   1   1   0   0   1   0
3   0   0   1   1   0   0   1
> ## No cycle
>  cycleMatrix(UG(~a*b))
NULL
> ## two cycles: the first is even and the second is odd
> cm <- cycleMatrix(UG(~a*b+b*c+c*d+d*a+a*u*v))
> apply(cm, 1, sum)
1 2 
4 3 
> 
> 
> 
> cleanEx(); ..nameEx <- "dSep"
> 
> ### * dSep
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: dSep
> ### Title: d-separation
> ### Aliases: dSep
> ### Keywords: graphs models multivariate
> 
> ### ** Examples
> 
> ## Conditioning on a transition node
> dSep(DAG(y ~ x, x ~ z), first="y", second="z", cond = "x")
[1] TRUE
> ## Conditioning on a collision node (collider)
> dSep(DAG(y ~ x, y ~ z), first="x", second="z", cond = "y")
[1] FALSE
> ## Conditioning on a source node
> dSep(DAG(y ~ x, z ~ x), first="y", second="z", cond = "x")
[1] TRUE
> ## Marginal independence
> dSep(DAG(y ~ x, y ~ z), first="x", second="z", cond = NULL)
[1] TRUE
> ## The DAG defined on p.~47 of Lauritzen (1996)
> dag <- DAG(g ~ x, h ~ x+f, f ~ b, x ~ l+d, d ~ c, c ~ a, l ~ y, y ~ b)
> dSep(dag, first="a", second="b", cond=c("x", "y"))
[1] TRUE
> dSep(dag, first="a", second=c("b", "d"), cond=c("x", "y"))
[1] FALSE
> 
> 
> 
> cleanEx(); ..nameEx <- "drawGraph"
> 
> ### * drawGraph
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: drawGraph
> ### Title: Drawing a graph with a simple point and click interface.
> ### Aliases: drawGraph
> ### Keywords: graphs hplot iplot
> 
> ### ** Examples
> 
> ## A directed acyclic graph
> d <- DAG(y1 ~ y2+y6, y2 ~ y3, y3 ~ y5+y6, y4 ~ y5+y6)
> ## Not run: drawGraph(d)
> 
> ## An undirected graph
> g <- UG(~giova*anto*armo + anto*arj*sara) 
> ## Not run: drawGraph(d)
> 
> ## An ancestral graph
> ag <- makeAG(ug=UG(~y0*y1), dag=DAG(y4~y2, y2~y1), bg=UG(~y2*y3+y3*y4))
> ## Not run: drawGraph(ag)
> ## Not run: drawGraph(ag, bid=FALSE)
> 
> ## A more complex example with coordinates: the UNIX evolution
> xy <-
+ structure(c(5, 15, 23, 25, 26, 17, 8, 6, 6, 7, 39, 33, 23, 49, 
+ 19, 34, 13, 29, 50, 68, 70, 86, 89, 64, 81, 45, 64, 49, 64, 87, 
+ 65, 65, 44, 37, 64, 68, 73, 85, 83, 95, 84, 0, 7, 15, 27, 44, 
+ 37, 36, 20, 51, 65, 44, 64, 59, 73, 69, 78, 81, 90, 97, 89, 72, 
+ 85, 74, 62, 68, 59, 52, 48, 43, 50, 34, 21, 18, 5, 1, 10, 2, 
+ 11, 2, 1, 44), .Dim = c(41, 2), .Dimnames = list(NULL, c("x", 
+ "y")))
> Unix <- DAG(
+                 SystemV.3 ~ SystemV.2,
+                 SystemV.2 ~ SystemV.0,
+                 SystemV.0 ~ TS4.0,
+                 TS4.0 ~ Unix.TS3.0 + Unix.TS.PP + CB.Unix.3,
+                 PDP11.SysV ~ CB.Unix.3,
+                 CB.Unix.3 ~ CB.Unix.2,
+                 CB.Unix.2 ~ CB.Unix.1,
+                 Unix.TS.PP ~ CB.Unix.3,
+                 Unix.TS3.0 ~ Unix.TS1.0 + PWB2.0 + USG3.0 + Interdata,
+                 USG3.0 ~ USG2.0,
+                 PWB2.0 ~ Interdata + PWB1.2,
+                 USG2.0 ~ USG1.0,
+                 CB.Unix.1 ~ USG1.0,
+                 PWB1.2 ~ PWB1.0,
+                 USG1.0 ~ PWB1.0,
+                 PWB1.0 ~ FifthEd,
+                 SixthEd ~ FifthEd,
+                 LSX ~ SixthEd,
+                 MiniUnix ~ SixthEd,
+                 Interdata ~ SixthEd,
+                 Wollongong ~ SixthEd,
+                 SeventhEd ~ Interdata,
+                 BSD1 ~ SixthEd,
+                 Xenix ~ SeventhEd,
+                 V32 ~ SeventhEd,
+                 Uniplus ~ SeventhEd,
+                 BSD3 ~ V32,
+                 BSD2 ~ BSD1,
+                 BSD4 ~ BSD3,
+                 BSD4.1 ~ BSD4,
+                 EigthEd ~ SeventhEd + BSD4.1,
+                 NinethEd ~ EigthEd,
+                 Ultrix32 ~ BSD4.2,
+                 BSD4.2 ~ BSD4.1,
+                 BSD4.3 ~ BSD4.2,
+                 BSD2.8 ~ BSD4.1 + BSD2,
+                 BSD2.9 ~ BSD2.8,
+                 Ultrix11 ~ BSD2.8 + V7M + SeventhEd,
+                 V7M ~ SeventhEd
+                 )
> drawGraph(Unix, coor=xy, adjust=FALSE)
> # dev.print(file="unix.fig", device=xfig) # Edit the graph with Xfig
> 
> 
> 
> cleanEx(); ..nameEx <- "edgeMatrix"
> 
> ### * edgeMatrix
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: edgeMatrix
> ### Title: Edge matrix of a graph
> ### Aliases: edgeMatrix
> ### Keywords: array algebra graphs multivariate
> 
> ### ** Examples
> 
> amat <- DAG(y ~ x+z, z~u+v)
> amat
  y x z u v
y 0 0 0 0 0
x 1 0 0 0 0
z 1 0 0 0 0
u 0 0 1 0 0
v 0 0 1 0 0
> edgeMatrix(amat)
  y x z u v
y 1 1 1 0 0
x 0 1 0 0 0
z 0 0 1 1 1
u 0 0 0 1 0
v 0 0 0 0 1
> edgeMatrix(amat, inv=TRUE)
  y x z u v
y 1 1 1 0 0
x 0 1 0 0 0
z 0 0 1 1 1
u 0 0 0 1 0
v 0 0 0 0 1
> 
> 
> 
> cleanEx(); ..nameEx <- "essentialGraph"
> 
> ### * essentialGraph
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: essentialGraph
> ### Title: Essential graph
> ### Aliases: essentialGraph
> ### Keywords: graphs models multivariate
> 
> ### ** Examples
> 
> dag = DAG(U ~ Y+Z, Y~X, Z~X)
> essentialGraph(dag)
  U Y Z X
U 0 0 0 0
Y 1 0 0 1
Z 1 0 0 1
X 0 1 1 0
> 
> 
> 
> cleanEx(); ..nameEx <- "findPath"
> 
> ### * findPath
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: findPath
> ### Title: Finding paths
> ### Aliases: findPath
> ### Keywords: graphs
> 
> ### ** Examples
> 
> ## A (single) path on a spanning tree
> findPath(bfs(UG(~ a*b*c + b*d + d*e+ e*c))$tree, st=1, en=5)
[1] 1 3 5
> 
> 
> 
> cleanEx(); ..nameEx <- "fitAncestralGraph"
> 
> ### * fitAncestralGraph
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: fitAncestralGraph
> ### Title: Fitting of Gaussian Ancestral Graph Models
> ### Aliases: fitAncestralGraph
> ### Keywords: graphs models multivariate
> 
> ### ** Examples
> 
> ## A covariance matrix
> "S" <- structure(c(2.93, -1.7, 0.76, -0.06,
+                   -1.7, 1.64, -0.78, 0.1,
+                    0.76, -0.78, 1.66, -0.78,
+                   -0.06, 0.1, -0.78, 0.81), .Dim = c(4,4),
+                  .Dimnames = list(c("y", "x", "z", "u"), c("y", "x", "z", "u")))
> ## The following should give the same fit.   
> ## Fit an ancestral graph y -> x <-> z <- u
> fitAncestralGraph(ag1 <- makeAG(dag=DAG(x~y,z~u), bg = UG(~x*z)), S, n=100)
$Shat
          y          x          z          u
y  2.930000 -1.4344254  0.0000000  0.0000000
x -1.434425  1.3799680 -0.3430373  0.0000000
z  0.000000 -0.3430373  1.5943070 -0.7442518
u  0.000000  0.0000000 -0.7442518  0.8100000

$Lhat
     y x z    u
y 2.93 0 0 0.00
x 0.00 0 0 0.00
z 0.00 0 0 0.00
u 0.00 0 0 0.81

$Bhat
         y x z         u
y 1.000000 0 0 0.0000000
x 0.489565 1 0 0.0000000
z 0.000000 0 1 0.9188294
u 0.000000 0 0 1.0000000

$Ohat
  y          x          z u
y 0  0.0000000  0.0000000 0
x 0  0.6777235 -0.3430373 0
z 0 -0.3430373  0.9104666 0
u 0  0.0000000  0.0000000 0

$dev
[1] 21.57711

$df
[1] 3

$it
[1] 4

> 
> ## Fit an ancestral graph y <-> x <-> z <-> u
> fitAncestralGraph(ag2 <- makeAG(bg= UG(~y*x+x*z+z*u)), S, n=100)
$Shat
          y          x          z          u
y  2.930000 -1.4344255  0.0000000  0.0000000
x -1.434425  1.3799680 -0.3430373  0.0000000
z  0.000000 -0.3430373  1.5943070 -0.7442518
u  0.000000  0.0000000 -0.7442518  0.8100000

$Lhat
          y          x          z          u
y  2.930000 -1.4344255  0.0000000  0.0000000
x -1.434425  1.3799680 -0.3430373  0.0000000
z  0.000000 -0.3430373  1.5943070 -0.7442518
u  0.000000  0.0000000 -0.7442518  0.8100000

$Bhat
  y x z u
y 1 0 0 0
x 0 1 0 0
z 0 0 1 0
u 0 0 0 1

$Ohat
          y          x          z          u
y  2.930000 -1.4344255  0.0000000  0.0000000
x -1.434425  1.3799680 -0.3430373  0.0000000
z  0.000000 -0.3430373  1.5943070 -0.7442518
u  0.000000  0.0000000 -0.7442518  0.8100000

$dev
[1] 21.57711

$df
[1] 3

$it
[1] 14

> 
> ## Fit the same graph with fitCovGraph
> fitCovGraph(ag2, S, n=100)    
$Shat
          y          x          z          u
y  2.930000 -1.4344255  0.0000000  0.0000000
x -1.434425  1.3799680 -0.3430373  0.0000000
z  0.000000 -0.3430373  1.5943070 -0.7442518
u  0.000000  0.0000000 -0.7442518  0.8100000

$dev
[1] 21.57711

$df
[1] 3

$it
[1] 14

> 
> ## Another example for the mathematics marks data
> 
> data(marks)
> S <- var(marks)
> mag1 <- makeAG(bg=UG(~mechanics*vectors*algebra+algebra*analysis*statistics))
> fitAncestralGraph(mag1, S, n=88)
$Shat
           mechanics   vectors  algebra  analysis statistics
mechanics  305.68848 127.04336 53.15421   0.00000    0.00000
vectors    127.04336 172.84222 43.11507   0.00000    0.00000
algebra     53.15421  43.11507 88.39060  84.78364   92.65154
analysis     0.00000   0.00000 84.78364 220.38036  155.53553
statistics   0.00000   0.00000 92.65154 155.53553  297.75536

$Lhat
           mechanics   vectors  algebra  analysis statistics
mechanics  305.68848 127.04336 53.15421   0.00000    0.00000
vectors    127.04336 172.84222 43.11507   0.00000    0.00000
algebra     53.15421  43.11507 88.39060  84.78364   92.65154
analysis     0.00000   0.00000 84.78364 220.38036  155.53553
statistics   0.00000   0.00000 92.65154 155.53553  297.75536

$Bhat
           mechanics vectors algebra analysis statistics
mechanics          1       0       0        0          0
vectors            0       1       0        0          0
algebra            0       0       1        0          0
analysis           0       0       0        1          0
statistics         0       0       0        0          1

$Ohat
           mechanics   vectors  algebra  analysis statistics
mechanics  305.68848 127.04336 53.15421   0.00000    0.00000
vectors    127.04336 172.84222 43.11507   0.00000    0.00000
algebra     53.15421  43.11507 88.39060  84.78364   92.65154
analysis     0.00000   0.00000 84.78364 220.38036  155.53553
statistics   0.00000   0.00000 92.65154 155.53553  297.75536

$dev
[1] 31.97538

$df
[1] 4

$it
[1] 21

> 
> mag2 <- makeAG(ug=UG(~mechanics*vectors+analysis*statistics),
+                dag=DAG(algebra~mechanics+vectors+analysis+statistics))
> fitAncestralGraph(mag2, S, n=88) # Same fit as above
$Shat
           mechanics   vectors  algebra  analysis statistics
mechanics  305.68848 127.04336 53.15421   0.00000    0.00000
vectors    127.04336 172.84222 43.11507   0.00000    0.00000
algebra     53.15421  43.11507 88.39060  84.78364   92.65154
analysis     0.00000   0.00000 84.78364 220.38036  155.53553
statistics   0.00000   0.00000 92.65154 155.53553  297.75536

$Lhat
           mechanics  vectors algebra analysis statistics
mechanics   305.6885 127.0434       0   0.0000     0.0000
vectors     127.0434 172.8422       0   0.0000     0.0000
algebra       0.0000   0.0000       0   0.0000     0.0000
analysis      0.0000   0.0000       0 220.3804   155.5355
statistics    0.0000   0.0000       0 155.5355   297.7554

$Bhat
            mechanics    vectors algebra   analysis statistics
mechanics   1.0000000  0.0000000       0  0.0000000  0.0000000
vectors     0.0000000  1.0000000       0  0.0000000  0.0000000
algebra    -0.1010961 -0.1751394       1 -0.2615174 -0.1745604
analysis    0.0000000  0.0000000       0  1.0000000  0.0000000
statistics  0.0000000  0.0000000       0  0.0000000  1.0000000

$Ohat
           mechanics vectors  algebra analysis statistics
mechanics          0       0  0.00000        0          0
vectors            0       0  0.00000        0          0
algebra            0       0 37.12009        0          0
analysis           0       0  0.00000        0          0
statistics         0       0  0.00000        0          0

$dev
[1] 31.97538

$df
[1] 4

$it
[1] 2

> 
> 
> 
> cleanEx(); ..nameEx <- "fitConGraph"
> 
> ### * fitConGraph
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: fitConGraph
> ### Title: Fitting of Gaussian concentration graph models
> ### Aliases: fitConGraph
> ### Keywords: graphs models multivariate
> 
> ### ** Examples
> 
> ## A model for the sample covariance matrix of the
> ## mathematics marks (Whittaker, 1990)
> data(marks)
> S <- cov(marks) * 87 / 88
> ## A butterfly concentration graph
> fitConGraph(UG(~ mechanics*vectors*algebra + algebra*analysis*statistics), S , n = 88)
$Shat
           mechanics   vectors   algebra  analysis statistics
mechanics  302.21475 125.59969 100.31599  99.62942   108.3001
vectors    125.59969 170.87810  84.18957  83.61337    90.8902
algebra    100.31599  84.18957 111.60318 110.83936   120.4857
analysis    99.62942  83.61337 110.83936 217.87603   153.7681
statistics 108.30013  90.89021 120.48567 153.76808   294.3718

$dev
[1] 0.900883

$df
[1] 4

$it
[1] 2

> 
> 
> 
> cleanEx(); ..nameEx <- "fitCovGraph"
> 
> ### * fitCovGraph
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: fitCovGraph
> ### Title: Fitting of Gaussian covariance graph models
> ### Aliases: fitCovGraph
> ### Keywords: graphs models multivariate
> 
> ### ** Examples
> 
> ## Correlations among four strategies to cope with stress for 
> ## 72 students. Cox & Wermuth (1996), p. 73.
> ##  Y = cognitive avoidance
> ##  X = vigilance
> ##  V = blunting
> ##  U = monitoring
> 
> R <- matrix(c(
+    1.00, -0.20,  0.46,  0.01,
+   -0.20,  1.00,  0.00,  0.47,
+    0.46,  0.00,  1.00, -0.15,
+    0.01,  0.47, -0.15,  1.00), 4, 4)
> nam <- c("Y", "X", "V", "U") 
> dimnames(R) <- list(nam, nam)
> 
> ## A chordless 4-cycle covariance graph
> gr <- UG(~ Y*X + X*U + U*V + V*Y)
> fitCovGraph(gr, R, n=72)
$Shat
           Y          X          V          U
Y  0.9995028 -0.2038079  0.4610985  0.0000000
X -0.2038079  1.0020514  0.0000000  0.4717660
V  0.4610985  0.0000000  1.0015117 -0.1536671
U  0.0000000  0.4717660 -0.1536671  0.9996165

$dev
[1] 0.007675671

$df
[1] 2

$it
[1] 7

> fitCovGraph(gr, R, n=72, alg="dual")
$Shat
           Y          X          V          U
Y  0.9994039 -0.2037920  0.4610611  0.0000000
X -0.2037920  1.0020263  0.0000000  0.4717259
V  0.4610611  0.0000000  1.0014904 -0.1536580
U  0.0000000  0.4717259 -0.1536580  0.9995170

$dev
[1] 0.007676413

$df
[1] 2

$it
[1] 6

> 
> 
> 
> cleanEx(); ..nameEx <- "fitDag"
> 
> ### * fitDag
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: fitDag
> ### Title: Fitting of Gaussian DAG models
> ### Aliases: fitDag
> ### Keywords: graphs models multivariate
> 
> ### ** Examples
> 
> dag <- DAG(y ~ x+u, x ~ z, z ~ u)
> "S" <- structure(c(2.93, -1.7, 0.76, -0.06,
+                    -1.7, 1.64, -0.78, 0.1,
+                     0.76, -0.78, 1.66, -0.78,
+                     -0.06, 0.1, -0.78, 0.81), .Dim = c(4,4),
+          .Dimnames = list(c("y", "x", "z", "u"), c("y", "x", "z", "u")))
> fitDag(dag, S, 200)
$Shat
           y         x          z          u
y  2.8998982 -1.685527  0.7687591 -0.3371388
x -1.6855265  1.640000 -0.7800000  0.3665060
z  0.7687591 -0.780000  1.6600000 -0.7800000
u -0.3371388  0.366506 -0.7800000  0.8100000

$Ahat
  y        x         z           u
y 1 1.039897 0.0000000 -0.05430825
x 0 1.000000 0.4698795  0.00000000
z 0 0.000000 1.0000000  0.96296296
u 0 0.000000 0.0000000  1.00000000

$Dhat
        y         x         z         u 
1.1654339 1.2734940 0.9088889 0.8100000 

$dev
[1] 26.90038

$df
[1] 2

> 
> 
> 
> cleanEx(); ..nameEx <- "fitDagLatent"
> 
> ### * fitDagLatent
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: fitDagLatent
> ### Title: Fitting Gaussian DAG models with one latent variable
> ### Aliases: fitDagLatent
> ### Keywords: graphs models multivariate
> 
> ### ** Examples
> 
> ## data from Joreskog and Goldberger (1975)
> V <- matrix(c(1,     0.36,   0.21,  0.10,  0.156, 0.158,
+               0.36,  1,      0.265, 0.284, 0.192, 0.324,
+               0.210, 0.265,  1,     0.176, 0.136, 0.226,
+               0.1,   0.284,  0.176, 1,     0.304, 0.305, 
+               0.156, 0.192,  0.136, 0.304, 1,     0.344,
+               0.158, 0.324,  0.226, 0.305, 0.344, 1),     6,6)
> nod <- c("y1", "y2", "y3", "x1", "x2", "x3")
> dimnames(V) <- list(nod,nod)
> dag <- DAG(y1 ~ z, y2 ~ z, y3 ~ z, z ~ x1 + x2 + x3, x1~x2+x3, x2~x3) 
> fitDagLatent(dag, V, n=530, latent="z", seed=4564)
...............................................................................
..............................................................
$Shat
          y1        y2        y3        x1        x2        x3         z
y1 1.0000000 0.3424515 0.1870304 0.1692573 0.1315980 0.2038965 0.4656999
y2 0.3424515 1.0000000 0.2953241 0.2672602 0.2077955 0.3219560 0.7353480
y3 0.1870304 0.2953241 1.0000000 0.1459645 0.1134878 0.1758367 0.4016114
x1 0.1692573 0.2672602 0.1459645 1.0000000 0.3040000 0.3050000 0.3634472
x2 0.1315980 0.2077955 0.1134878 0.3040000 1.0000000 0.3440000 0.2825812
x3 0.2038965 0.3219560 0.1758367 0.3050000 0.3440000 1.0000000 0.4378280
z  0.4656999 0.7353480 0.4016114 0.3634472 0.2825812 0.4378280 1.0000000

$Ahat
   y1 y2 y3         x1          x2         x3          z
y1  1  0  0  0.0000000  0.00000000  0.0000000 -0.4656999
y2  0  1  0  0.0000000  0.00000000  0.0000000 -0.7353480
y3  0  0  1  0.0000000  0.00000000  0.0000000 -0.4016114
x1  0  0  0  1.0000000 -0.22580030 -0.2273247  0.0000000
x2  0  0  0  0.0000000  1.00000000 -0.3440000  0.0000000
x3  0  0  0  0.0000000  0.00000000  1.0000000  0.0000000
z   0  0  0 -0.2321414 -0.09726341 -0.3335663  1.0000000

$Dhat
       y1        y2        y3        x1        x2        x3         z 
0.7831236 0.4592633 0.8387083 0.8620227 0.8816640 1.0000000 0.7420994 

$dev
[1] 12.52184

$df
[1] 6

$it
[1] 142

> fitDagLatent(dag, V, n=530, latent="z", norm=2, seed=145)
...............................................................................
........................................................................
$Shat
           y1         y2         y3         x1         x2         x3          z
y1  1.0000000  0.3424514  0.1870305  0.1692574  0.1315981  0.2038965 -0.5405989
y2  0.3424514  1.0000000  0.2953241  0.2672601  0.2077955  0.3219559 -0.8536145
y3  0.1870305  0.2953241  1.0000000  0.1459646  0.1134879  0.1758368 -0.4662031
x1  0.1692574  0.2672601  0.1459646  1.0000000  0.3040000  0.3050000 -0.4219009
x2  0.1315981  0.2077955  0.1134879  0.3040000  1.0000000  0.3440000 -0.3280291
x3  0.2038965  0.3219559  0.1758368  0.3050000  0.3440000  1.0000000 -0.5082445
z  -0.5405989 -0.8536145 -0.4662031 -0.4219009 -0.3280291 -0.5082445  1.3475286

$Ahat
   y1 y2 y3       x1         x2         x3         z
y1  1  0  0 0.000000  0.0000000  0.0000000 0.4011781
y2  0  1  0 0.000000  0.0000000  0.0000000 0.6334667
y3  0  0  1 0.000000  0.0000000  0.0000000 0.3459690
x1  0  0  0 1.000000 -0.2258003 -0.2273247 0.0000000
x2  0  0  0 0.000000  1.0000000 -0.3440000 0.0000000
x3  0  0  0 0.000000  0.0000000  1.0000000 0.0000000
z   0  0  0 0.269477  0.1129064  0.3872142 1.0000000

$Dhat
       y1        y2        y3        x1        x2        x3         z 
0.7831236 0.4592637 0.8387082 0.8620227 0.8816640 1.0000000 1.0000000 

$dev
[1] 12.52184

$df
[1] 6

$it
[1] 152

> 
> 
> 
> cleanEx(); ..nameEx <- "fundCycles"
> 
> ### * fundCycles
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: fundCycles
> ### Title: Fundamental cycles
> ### Aliases: fundCycles
> ### Keywords: graphs models multivariate
> 
> ### ** Examples
> 
> ## Three fundamental cycles
> fundCycles(UG(~a*b*d + d*e + e*a*f))
[[1]]
 [,1] [,2]
    2    1
    1    3
    3    2

[[2]]
 [,1] [,2]
    3    1
    1    4
    4    3

[[3]]
 [,1] [,2]
    4    1
    1    5
    5    4

> 
> 
> 
> cleanEx(); ..nameEx <- "glucose"
> 
> ### * glucose
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: glucose
> ### Title: Glucose control
> ### Aliases: glucose
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(glucose)
> ## See Cox & Wermuth (1996), Figure 6.3 p. 140
> coplot(Y ~ W | A, data=glucose)
> 
> 
> 
> cleanEx(); ..nameEx <- "isAcyclic"
> 
> ### * isAcyclic
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: isAcyclic
> ### Title: Graph queries
> ### Aliases: isAcyclic
> ### Keywords: graphs models multivariate
> 
> ### ** Examples
> 
> ## A cyclic graph
> d <- matrix(0,3,3)
> rownames(d) <- colnames(d) <- c("x", "y", "z")
> d["x","y"] <- d["y", "z"] <- d["z", "x"] <- 1
> ## Test if the graph is acyclic
> isAcyclic(d)
[1] FALSE
> 
> 
> 
> cleanEx(); ..nameEx <- "isGident"
> 
> ### * isGident
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: isGident
> ### Title: G-identifiability of an UG
> ### Aliases: isGident
> ### Keywords: graphs models multivariate
> 
> ### ** Examples
> 
> ## A not G-identifiable UG
> G1 <- UG(~ a*b + u*v)
> isGident(G1)
[1] FALSE
> ## G-identifiable UG
> G2 <- UG(~ a + b + u*v)
> isGident(G2)
[1] TRUE
> ## G-identifiable UG
> G3 <- cmpGraph(UG(~a*b*c+x*y*z))
> isGident(G3)
[1] TRUE
> 
> 
> 
> cleanEx(); ..nameEx <- "makeAG"
> 
> ### * makeAG
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: makeAG
> ### Title: Ancestral Graphs
> ### Aliases: makeAG
> ### Keywords: graphs models multivariate
> 
> ### ** Examples
> 
> ## Examples from Richardson and Spirtes (2002)
> ## Not run: a1 <- makeAG(dag=DAG(a~b, b~d, d~c), bg=UG(~a*c))  # Not an AG. (a2) p.969
> a2 <- makeAG(dag=DAG(b ~ a, d~c), bg=UG(~a*c+c*b+b*d))           # Fig. 3 (b1) p.969
>  a3 <- makeAG(ug = UG(~ a*c), dag=DAG(b ~ a, d~c), bg=UG(~ b*d)) # Fig. 3 (b2) p.969
>  a5 <- makeAG(bg=UG(~alpha*beta+gamma*delta), dag=DAG(alpha~gamma,
+ delta~beta))  # Fig. 6 p. 973
> ## Another Example
>  a4 <- makeAG(ug=UG(~y0*y1), dag=DAG(y4~y2, y2~y1), bg=UG(~y2*y3+y3*y4))
> 
> 
> 
> cleanEx(); ..nameEx <- "marks"
> 
> ### * marks
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: marks
> ### Title: Mathematics marks
> ### Aliases: marks
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(marks)
> pairs(marks)
> 
> 
> 
> cleanEx(); ..nameEx <- "parcor"
> 
> ### * parcor
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: parcor
> ### Title: Partial correlations
> ### Aliases: parcor
> ### Keywords: array graphs models multivariate
> 
> ### ** Examples
> 
> ### Partial correlations for the mathematics marks data
> data(marks)
> S <- var(marks)
> parcor(S)
               mechanics    vectors   algebra      analysis statistics
mechanics   1.0000000000 0.32846282 0.2292442 -0.0007121818 0.02550897
vectors     0.3284628169 1.00000000 0.2816015  0.0778303605 0.01996942
algebra     0.2292441885 0.28160149 1.0000000  0.4317713938 0.35673607
analysis   -0.0007121818 0.07783036 0.4317714  1.0000000000 0.25277656
statistics  0.0255089719 0.01996942 0.3567361  0.2527765574 1.00000000
> 
> 
> 
> cleanEx(); ..nameEx <- "pcor"
> 
> ### * pcor
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: pcor
> ### Title: Partial correlation
> ### Aliases: pcor
> ### Keywords: models multivariate
> 
> ### ** Examples
> 
> data(marks)
> ## The correlation between vectors and algebra given analysis and statistics
>  pcor(c("vectors", "algebra", "analysis", "statistics"), var(marks))
[1] 0.3882030
> ## The same
> pcor(c(2,3,4,5), var(marks))
[1] 0.3882030
> ## The correlation between vectors and algebra given statistics
>  pcor(c("vectors", "algebra", "statistics"), var(marks))
[1] 0.4753595
> ## The marginal correlation between analysis and statistics 
> pcor(c("analysis","statistics"), var(marks))
[1] 0.6071743
> 
> 
> 
> cleanEx(); ..nameEx <- "pcor.test"
> 
> ### * pcor.test
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: pcor.test
> ### Title: Test for zero partial association
> ### Aliases: pcor.test
> ### Keywords: htest multivariate
> 
> ### ** Examples
> 
> ## Are 2,3 independent given 1?
> data(marks)
> pcor.test(pcor(c(2,3,1), var(marks)), 1, n=88)
$tval
[1] 4.528232

$df
[1] 85

$pvalue
[1] 1.923857e-05

> 
> 
> 
> cleanEx(); ..nameEx <- "rcorr"
> 
> ### * rcorr
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: rcorr
> ### Title: Random correlation matrix
> ### Aliases: rcorr
> ### Keywords: distribution multivariate
> 
> ### ** Examples
> 
> ## A random correlation matrix of order 3
> rcorr(3)
           [,1]        [,2]         [,3]
[1,]  1.0000000 0.521541162 -0.217068990
[2,]  0.5215412 1.000000000  0.001206383
[3,] -0.2170690 0.001206383  1.000000000
> ## A random correlation matrix of order 5
> rcorr(5)
           [,1]       [,2]         [,3]         [,4]        [,5]
[1,] 1.00000000  0.1684914  0.139394244  0.191157141  0.02119760
[2,] 0.16849139  1.0000000  0.590519025 -0.071782494 -0.75399064
[3,] 0.13939424  0.5905190  1.000000000 -0.005817047 -0.84219944
[4,] 0.19115714 -0.0717825 -0.005817047  1.000000000  0.45304494
[5,] 0.02119760 -0.7539906 -0.842199438  0.453044941  1.00000000
> 
> 
> 
> cleanEx(); ..nameEx <- "rnormDag"
> 
> ### * rnormDag
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: rnormDag
> ### Title: Random sample from a decomposable Gaussian model
> ### Aliases: rnormDag
> ### Keywords: distribution multivariate
> 
> ### ** Examples
> 
> ## Generate a sample of 100 observation from a multivariate normal
> ## The matrix of the path coefficients 
> A <- matrix(
+ c(1, -2, -3,  0, 0,  0,  0,
+   0,  1,  0, -4, 0,  0,  0,
+   0,  0,  1,  2, 0,  0,  0,
+   0,  0,  0,  1, 1, -5,  0,
+   0,  0,  0,  0, 1,  0,  3,
+   0,  0,  0,  0, 0,  1, -4,
+   0,  0,  0,  0, 0,  0,  1), 7, 7, byrow=TRUE)
> D <- rep(1, 7)
> X <- rnormDag(100, A, D)
> 
> ## The true covariance matrix
> solve(A) %*% diag(D) %*% t(solve(A))
      [,1]  [,2]  [,3]  [,4] [,5] [,6] [,7]
[1,]  2238  4450 -2221  1112 -140  194   46
[2,]  4450  8897 -4448  2224 -280  388   92
[3,] -2221 -4448  2225 -1112  140 -194  -46
[4,]  1112  2224 -1112   556  -70   97   23
[5,]  -140  -280   140   -70   10  -12   -3
[6,]   194   388  -194    97  -12   17    4
[7,]    46    92   -46    23   -3    4    1
> 
> ## Triangular decomposition of the sample covariance matrix
> triDec(cov(X))$A
     [,1]      [,2]       [,3]        [,4]        [,5]       [,6]       [,7]
[1,]    1 -1.983356 -3.0155982 -0.04543888 -0.05244152 -0.1656270 -0.7566530
[2,]    0  1.000000  0.0465952 -3.83744867 -0.04859890 -0.3316160 -0.3175740
[3,]    0  0.000000  1.0000000  1.89111947 -0.13189087  0.5756604 -0.1721831
[4,]    0  0.000000  0.0000000  1.00000000  0.97995432 -4.9036606 -0.3272093
[5,]    0  0.000000  0.0000000  0.00000000  1.00000000  0.1104357  2.4880929
[6,]    0  0.000000  0.0000000  0.00000000  0.00000000  1.0000000 -4.0068651
[7,]    0  0.000000  0.0000000  0.00000000  0.00000000  0.0000000  1.0000000
> 
> 
> 
> cleanEx(); ..nameEx <- "rsphere"
> 
> ### * rsphere
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: rsphere
> ### Title: Random vectors on a sphere
> ### Aliases: rsphere
> ### Keywords: distribution multivariate
> 
> ### ** Examples
> 
> ## 100 points on circle
> z <- rsphere(100,2)
> plot(z)
> 
> ## 100 points on a sphere
> z <- rsphere(100, 3)
> pairs(z)
> 
> 
> 
> cleanEx(); ..nameEx <- "shipley.test"
> 
> ### * shipley.test
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: shipley.test
> ### Title: Test of all independencies implied by a given DAG
> ### Aliases: shipley.test
> ### Keywords: graphs models multivariate
> 
> ### ** Examples
> 
> ## A decomposable model for the mathematics marks data
> data(marks)
> dag <- DAG(mechanics ~ vectors+algebra, vectors ~ algebra, statistics ~ algebra+analysis, analysis ~ algebra)
> shipley.test(dag, cov(marks), n=88)
$ctest
[1] 2.850461

$df
[1] 8

$pvalue
[1] 0.9433884

> 
> 
> 
> cleanEx(); ..nameEx <- "swp"
> 
> ### * swp
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: swp
> ### Title: Sweep operator
> ### Aliases: swp
> ### Keywords: array algebra models multivariate
> 
> ### ** Examples
> 
> ## A very simple example
> V <- matrix(c(10, 1, 1, 2), 2, 2)
> swp(V, 2)
     [,1] [,2]
[1,]  9.5  0.5
[2,]  0.5 -0.5
> 
> 
> 
> cleanEx(); ..nameEx <- "topSort"
> 
> ### * topSort
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: topSort
> ### Title: Topological sort
> ### Aliases: topSort topOrder
> ### Keywords: graphs models multivariate
> 
> ### ** Examples
> 
> ## A simple example
> dag <- DAG(a ~ b, c ~ a + b, d ~ c + b)
> dag
  a b c d
a 0 0 1 0
b 1 0 1 1
c 0 0 0 1
d 0 0 0 0
> topOrder(dag)
[1] 2 1 3 4
> topSort(dag)
  b a c d
b 0 1 1 1
a 0 0 1 0
c 0 0 0 1
d 0 0 0 0
> 
> 
> 
> cleanEx(); ..nameEx <- "transClos"
> 
> ### * transClos
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: transClos
> ### Title: Transitive closure of a graph
> ### Aliases: transClos
> ### Keywords: graphs models multivariate
> 
> ### ** Examples
> 
> ## Closure of a DAG
> d <- DAG(y ~ x, x ~ z)
> transClos(d)
  y x z
y 0 0 0
x 1 0 0
z 1 1 0
> 
> ## Closure of an UG
> g <- UG(~ x*y*z+z*u+u*v)
> transClos(g)
  x y z u v
x 0 1 1 1 1
y 1 0 1 1 1
z 1 1 0 1 1
u 1 1 1 0 1
v 1 1 1 1 0
> 
> 
> 
> cleanEx(); ..nameEx <- "triDec"
> 
> ### * triDec
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: triDec
> ### Title: Triangular decomposition of a covariance matrix
> ### Aliases: triDec
> ### Keywords: array algebra models multivariate
> 
> ### ** Examples
> 
> ## Triangular decomposition of a covariance matrix
> B <- matrix(c(1,  -2, 0, 1,
+               0,   1, 0, 1,
+               0,   0, 1, 0,
+               0,   0, 0, 1), 4, 4, byrow=TRUE)
> B
     [,1] [,2] [,3] [,4]
[1,]    1   -2    0    1
[2,]    0    1    0    1
[3,]    0    0    1    0
[4,]    0    0    0    1
> D <- diag(c(3, 1, 2, 1))
> S <- B %*% D %*% t(B)
> triDec(S)
$A
     [,1] [,2] [,3] [,4]
[1,]    1    2    0   -3
[2,]    0    1    0   -1
[3,]    0    0    1    0
[4,]    0    0    0    1

$B
     [,1] [,2] [,3] [,4]
[1,]    1   -2    0    1
[2,]    0    1    0    1
[3,]    0    0    1    0
[4,]    0    0    0    1

$Delta
[1] 3 1 2 1

> solve(B)
     [,1] [,2] [,3] [,4]
[1,]    1    2    0   -3
[2,]    0    1    0   -1
[3,]    0    0    1    0
[4,]    0    0    0    1
> 
> 
> 
> cleanEx(); ..nameEx <- "unmakeAG"
> 
> ### * unmakeAG
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: unmakeAG
> ### Title: Ancestral Graph components
> ### Aliases: unmakeAG
> ### Keywords: graphs models multivariate
> 
> ### ** Examples
> 
> ag <- makeAG(ug=UG(~y0*y1), dag=DAG(y4~y2, y2~y1), bg=UG(~y2*y3+y3*y4))
> unmakeAG(ag)
$dag
   y2 y3 y4 y0 y1
y2  0  0  1  0  0
y3  0  0  0  0  0
y4  0  0  0  0  0
y0  0  0  0  0  0
y1  1  0  0  0  0

$ug
   y2 y3 y4 y0 y1
y2  0  0  0  0  0
y3  0  0  0  0  0
y4  0  0  0  0  0
y0  0  0  0  0  1
y1  0  0  0  1  0

$bg
   y2 y3 y4 y0 y1
y2  0  1  0  0  0
y3  1  0  1  0  0
y4  0  1  0  0  0
y0  0  0  0  0  0
y1  0  0  0  0  0

> 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  5.77 0.07 5.91 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')