Rel 0.1.0

Author: goedman

docs/src/walkthrough.md

@@ -24,7 +24,7 @@ d_string = "DAG(
     statistics ~ algebra+analysis, 
     analysis ~ algebra)"
 
-dag = DAG("marks", d_string, df);
+dag = DAG("marks", d_string; df=df);
 show(dag)
 ```
 
@@ -50,7 +50,6 @@ OrderedDict{Symbol,Union{Nothing, Array{Symbol,1}, Symbol}} with 4 entries:
 
 Internally, a DAG object will always contain an OrderedDict representation of the DAG. This representation is used in all functions. In the definition of the OrderedDict, read `=>` as `~` in regression models or `<-` in causal models.
 
-
 Optional display the DAG using GraphViz:
 ```julia
 fname = ProjDir * "/marks.dot"
@@ -60,7 +59,7 @@ Sys.isapple() && run(`open -a GraphViz.app $(fname)`)
 
 The DAG pdf is [here](https://github.com/StatisticalRethinkingJulia/StructuralCausalModels.jl/blob/master/docs/src/marks.pdf).
 
-In this case a DataFrame with observed values has been provided and the related covariance matrix has been computed and stored in the DAG object:
+In this example a DataFrame with observed values has been provided and the related covariance matrix will be computed and stored in the DAG object:
 ```julia
 display(dag.s)
 
@@ -93,68 +92,27 @@ Although this initial version of StructuralCausalModels does not support latent
 
 # Directed separation
 
-Given a causal graph, `d_separation(dag, f, l, cond)` determines if the vertices in set `f` are `d-separated` from the vertices in set `l` given the conditioning set `cond`.
+Given a causal graph `dag`, `d_separation(dag, f, l, cond)` determines if the vertices in set `f` are `d-separated` from the vertices in set `l` given the conditioning set `cond`.
 
 Show several `d_separation` results for the marks model:
 ```julia
-f = [:statistics]; s = [:mechanics]; sel = vcat(f, s)
-cond = [:algebra]
+f = [:statistics]; s = [:mechanics];
 
-e = d_separation(dag, f, s, cond)
-println("d_separation($(dag.name), $f, $s, $cond) = $e")
+e = d_separation(dag, f, s; cond=:algebra);
+println("d_separation($(dag.name), $f, $s; cond=:algebra) = $e")
 
-e = d_separation(dag, f, s)
+e = d_separation(dag, f, s);
 println("d_separation($(dag.name), $f, $s) = $e")
 
-print("d_separation($(dag.name), [:statistics], [:mechanics], [:vectors]) = ")
-println(d_separation(dag, [:statistics], [:mechanics], [:vectors]))
-
-print("d_separation($(dag.name), [:statistics], [:mechanics], [:analysis, :vectors]) = ")
-println(d_separation(dag, [:statistics], [:mechanics], [:analysis, :vectors]))
-
-print("d_separation($(dag.name), [:statistics, :analysis], [:mechanics], [:algebra]) = ")
-println(d_separation(dag, [:statistics, :analysis], [:mechanics], [:algebra]))
-
-print("d_separation($(dag.name), [:statistics], [:mechanics, :vectors], [:algebra]) = ")
-println(d_separation(dag, [:statistics], [:mechanics, :vectors], [:algebra]))
-
-print("d_separation($(dag.name), [:statistics], [:mechanics, :analysis], [:algebra]) = ")
-println(d_separation(dag, [:statistics], [:mechanics, :analysis], [:algebra]))
-
-print("d_separation($(dag.name), [:analysis], [:vectors]) = ")
-println(d_separation(dag, [:analysis], [:vectors]))
-
-print("d_separation($(dag.name), [:analysis], [:vectors], [:algebra]) = ")
-println(d_separation(dag, [:analysis], [:vectors], [:algebra]))
-
-print("d_separation($(dag.name), [:vectors], [:statistics], [:algebra]) = ")
-println(d_separation(dag, [:analysis], [:vectors], [:algebra]))
-
-print("d_separation($(dag.name), [:statistics], [:algebra], [:analysis]) = ")
-println(d_separation(dag, [:statistics], [:algebra], [:analysis]))
-
-print("d_separation($(dag.name), [:statistics, :analysis], [:mechanics, :vectors]) = ")
-println(d_separation(dag, [:statistics, :analysis], [:mechanics, :vectors]))
-
-print("d_separation($(dag.name), [:statistics, :analysis], [:mechanics, :vectors], [:algebra]) = ")
-println(d_separation(dag, [:statistics, :analysis], [:mechanics, :vectors], [:algebra]))
+print("d_separation($(dag.name), [:statistics], [:mechanics, :analysis]; cond=[:algebra]) = ");
+println(d_separation(dag, [:statistics], [:mechanics, :analysis]; cond=[:algebra]))
 ```
 
 will produce:
 ```julia
-d_separation(marks, [:statistics], [:mechanics], [:algebra]) = false
+d_separation(marks, [:statistics], [:mechanics]; cond=:algebra) = false
 d_separation(marks, [:statistics], [:mechanics]) = true
-d_separation(marks, [:statistics], [:mechanics], [:vectors]) = true
-d_separation(marks, [:statistics], [:mechanics], [:analysis, :vectors]) = true
-d_separation(marks, [:statistics, :analysis], [:mechanics], [:algebra]) = false
-d_separation(marks, [:statistics], [:mechanics, :vectors], [:algebra]) = false
-d_separation(marks, [:statistics], [:mechanics, :analysis], [:algebra]) = false
-d_separation(marks, [:analysis], [:vectors]) = true
-d_separation(marks, [:analysis], [:vectors], [:algebra]) = false
-d_separation(marks, [:vectors], [:statistics], [:algebra]) = false
-d_separation(marks, [:statistics], [:algebra], [:analysis]) = false
-d_separation(marks, [:statistics, :analysis], [:mechanics, :vectors]) = true
-d_separation(marks, [:statistics, :analysis], [:mechanics, :vectors], [:algebra]) = false
+d_separation(marks, [:statistics], [:mechanics, :analysis]: cond=[:algebra]) = false
 ```
 
 # Basis set
@@ -176,19 +134,9 @@ BasisSet[
 ]
 ```
 
-# Shipley test
-
-Perform the Shipley test:
-```julia
-t = shipley_test(dag)
-display(t)
-
-(ctest = 2.816854003338401, df = 8, pv = 0.9453198036802164)
-```
-
 # Adjustment sets
 
-D_separation provides a set of conditional independencies given the causal model. The conditioning set closes (blocks) all paths. It provides ways to test the chosen causal model given observational data.
+The function  `basis_set()` provides a set of conditional independencies given the causal model. The conditioning set closes ("blocks the flow of causal info") all paths. The conditioning_set can be empty. It provides ways to test the chosen causal model given observational data.
 
 The function `adjustment_sets()` answers a related question, i.e. how to prevent confounding in multiple regression models assuming the chosen causal model is correct.
 
@@ -231,6 +179,16 @@ Adjustment sets:
  [:a, :m]
 ```
 
+# Shipley test
+
+Perform the Shipley test:
+```julia
+t = shipley_test(dag)
+display(t)
+
+(ctest = 2.816854003338401, df = 8, pv = 0.9453198036802164)
+```
+
 # Paths
 
 Adjustment_sets is based on several path manipulations and checks:

examples/dagitty/confounding_triangle.jl

@@ -3,21 +3,17 @@ using StructuralCausalModels
 ProjDir = @__DIR__
 cd(ProjDir) #do
 
-d = OrderedDict(
-  :w => :s,
-  :d => [:a, :w, :m],
-  :m => [:a, :s],
-  :a => :s
-);
-
 d_string = "dag {A -> {E Z}; B -> {D Z}; Z -> {D E}; E -> D}"
 
 dag = DAG("conf_triangles", d_string);
 show(dag)
 
+to_ggm(dag) |> display
+println()
+
 fname = ProjDir * "/conf_triangles.dot"
 to_graphviz(dag, fname)
-#Sys.isapple() && run(`open -a GraphViz.app $(fname)`)
+Sys.isapple() && run(`open -a GraphViz.app $(fname)`)
 
 println("Basis set:")
 bs = basis_set(dag)
@@ -30,13 +26,13 @@ A ⊥ B
 
 f = :A; s = :B;
 
-e = d_separation(dag, f, s, nothing)
-println("d_separation($(dag.name), $f, $s) = $e\n")
+e = d_separation(dag, f, s)
+fprintln("d_separation($(dag.name), $f, $s) = $e\n")
 
 f = [:B]; s = [:E];
 
-e = d_separation(dag, f, s, [:A, :Z])
-println("d_separation($(dag.name), $f, $s, [:A, :Z]) = $e\n")
+e = d_separation(dag, f, s; cond=[:A, :Z])
+println("d_separation($(dag.name), $f, $s; cond=[:A, :Z]) = $e\n")
 
 ap = all_paths(dag, :E, :D)
 ap |> display
@@ -46,10 +42,6 @@ bp = backdoor_paths(dag, ap, :E)
 bp |> display
 println()
 
-op = open_paths(dag, bp, Symbol[])
-op |> display
-println()
-
 adjustmentsets = adjustment_sets(dag, :E, :D)
 println("Adjustment sets:")
 adjustmentsets |> display

examples/marks/marks.jl

@@ -35,7 +35,7 @@ d_string = "DAG(
     statistics ~ algebra+analysis, 
     analysis ~ algebra)"
 
-dag = DAG("marks", d_string, df);
+dag = DAG("marks", d_string, df=df);
 show(dag)
 
 fname = ProjDir * "/marks.dot"
@@ -55,44 +55,44 @@ display("pcor_test = $pt"); println()
 
 f = [:statistics]; s = [:mechanics]; sel = vcat(f, s)
 
-e = d_separation(dag, f, s, :algebra)
-println("d_separation($(dag.name), $f, $s, :algebra) = $e")
+e = d_separation(dag, f, s; cond=:algebra)
+println("d_separation($(dag.name), $f, $s; cond=:algebra) = $e")
 
 e = d_separation(dag, f, s)
 println("d_separation($(dag.name), $f, $s) = $e")
 
-print("d_separation($(dag.name), [:statistics], [:mechanics], [:vectors]) = ")
-println(d_separation(dag, [:statistics], [:mechanics], [:vectors]))
+print("d_separation($(dag.name), [:statistics], [:mechanics]; cond=[:vectors]) = ")
+println(d_separation(dag, [:statistics], [:mechanics]; cond=[:vectors]))
 
-print("d_separation($(dag.name), [:statistics], [:mechanics], [:analysis, :vectors]) = ")
-println(d_separation(dag, [:statistics], [:mechanics], [:analysis, :vectors]))
+print("d_separation($(dag.name), [:statistics], [:mechanics]; cond=[:analysis, :vectors]) = ")
+println(d_separation(dag, [:statistics], [:mechanics]; cond=[:analysis, :vectors]))
 
-print("d_separation($(dag.name), [:statistics, :analysis], [:mechanics], [:algebra]) = ")
-println(d_separation(dag, [:statistics, :analysis], [:mechanics], [:algebra]))
+print("d_separation($(dag.name), [:statistics, :analysis], [:mechanics]; cond=[:algebra]) = ")
+println(d_separation(dag, [:statistics, :analysis], [:mechanics]; cond=[:algebra]))
 
-print("d_separation($(dag.name), [:statistics], [:mechanics, :vectors], [:algebra]) = ")
-println(d_separation(dag, [:statistics], [:mechanics, :vectors], [:algebra]))
+print("d_separation($(dag.name), [:statistics], [:mechanics, :vectors]; cond=[:algebra]) = ")
+println(d_separation(dag, [:statistics], [:mechanics, :vectors]; cond=[:algebra]))
 
-print("d_separation($(dag.name), [:statistics], [:mechanics, :analysis], [:algebra]) = ")
-println(d_separation(dag, [:statistics], [:mechanics, :analysis], [:algebra]))
+print("d_separation($(dag.name), [:statistics], [:mechanics, :analysis]; cond=[:algebra]) = ")
+println(d_separation(dag, [:statistics], [:mechanics, :analysis]; cond=[:algebra]))
 
 print("d_separation($(dag.name), [:analysis], [:vectors]) = ")
 println(d_separation(dag, [:analysis], [:vectors]))
 
-print("d_separation($(dag.name), [:analysis], [:vectors], [:algebra]) = ")
-println(d_separation(dag, [:analysis], [:vectors], [:algebra]))
+print("d_separation($(dag.name), [:analysis], [:vectors]; cond=[:algebra]) = ")
+println(d_separation(dag, [:analysis], [:vectors]; cond=[:algebra]))
 
-print("d_separation($(dag.name), [:vectors], [:statistics], [:algebra]) = ")
-println(d_separation(dag, [:analysis], [:vectors], [:algebra]))
+print("d_separation($(dag.name), [:vectors], [:statistics]; cond=[:algebra]) = ")
+println(d_separation(dag, [:analysis], [:vectors]; cond=[:algebra]))
 
-print("d_separation($(dag.name), [:statistics], [:algebra], [:analysis]) = ")
-println(d_separation(dag, [:statistics], [:algebra], [:analysis]))
+print("d_separation($(dag.name), [:statistics], [:algebra]; cond=[:analysis]) = ")
+println(d_separation(dag, [:statistics], [:algebra]; cond=[:analysis]))
 
 print("d_separation($(dag.name), [:statistics, :analysis], [:mechanics, :vectors]) = ")
 println(d_separation(dag, [:statistics, :analysis], [:mechanics, :vectors]))
 
-print("d_separation($(dag.name), [:statistics, :analysis], [:mechanics, :vectors], [:algebra]) = ")
-println(d_separation(dag, [:statistics, :analysis], [:mechanics, :vectors], [:algebra]))
+print("d_separation($(dag.name), [:statistics, :analysis], [:mechanics, :vectors]; cond=[:algebra]) = ")
+println(d_separation(dag, [:statistics, :analysis], [:mechanics, :vectors]; cond=[:algebra]))
 
 adjustmentsets = adjustment_sets(dag, :statistics, :mechanics)
 println("\nAdjustment sets:")

src/methods/d_separation.jl

@@ -18,7 +18,7 @@ d_separation(
 
 ### Keyword arguments
 ```julia
-* `cond::SymbolListOrNothing=nothing`  : Conditioning set
+* `cset::SymbolListOrNothing=nothing`  : Conditioning set
 * `debug=false`                        : Trace execution
 ```
 
@@ -45,7 +45,7 @@ d = OrderedDict(
 );
 
 dag = DAG("marks", d, df);
-d_separation(marks, [:statistics], [:mechanics], [:algebra]))
+d_separation(marks, [:statistics], [:mechanics]; cset=[:algebra]))
 ```
 ### Acknowledgements
 
@@ -62,9 +62,9 @@ The Julia translation is licenced under: MIT.
 Part of the API, exported.
 """
 function d_separation(d::DAG, first::SymbolList, second::SymbolList; 
-  cond::SymbolListOrNothing=nothing, debug=false)
-
-  e = induced_covariance_graph(d, vcat(first, second), SymbolList[]; debug=debug)
+  cset::SymbolListOrNothing=nothing, debug=false)
+  
+  e = induced_covariance_graph(d, vcat(first, second), cset; debug=debug)
   sum(e[first, second]) == 0
 
 end

src/methods/induced_covariance_graph.jl

@@ -57,45 +57,49 @@ $(SIGNATURES)
 
 Internal
 """
-function induced_covariance_graph(d::DAG, sel, cond; debug=false)
+function induced_covariance_graph(d::DAG, sel, cset; debug=false)
 
   @assert all([c in d.vars for c in sel]) "Selection nodes are not among vertices."
 
-  if isnothing(cond)
-    cond = []
-  elseif typeof(cond) == Symbol
-    cond = [cond]
-  elseif length(cond) > 0
-    @assert all([c in d.vars for c in cond]) "Conditioning nodes are not among vertices."
-    @assert !all([c in sel for c in cond]) "Conditioning nodes in selected nodes."
+  if isnothing(cset)
+    cs = []
+  elseif typeof(cset) == Symbol
+    cs = [cset]
+  elseif length(cset) > 0
+    cs = cset
   end
+  @assert all([c in d.vars for c in cs]) "Conditioning nodes are not among vertices."
+  #@assert !(length([cs in sel]) > 0) "Conditioning nodes in selected nodes."
 
-  l = setdiff(d.vars, union(sel, cond))
+  l = setdiff(d.vars, union(sel, cs))
   debug && println(l)
   l = union(l, sel)
   debug && println(l)
-  r = union(sel, cond)
+  r = union(sel, cs)
   debug && println(r)
 
   e = edge_matrix(d.a)                  # From adjacency matrix to edge matrix
   debug && println(e)
-  al = ancester_graph(e[l, l])
+  al = transpose(StructuralCausalModels.ancester_graph(e[l, l]))
   debug && println(al)
-  if length(cond) > 0
-    trl = indicator_matrix( e[cond, l] * al)
+  if length(cs) > 0
+    trl = StructuralCausalModels.indicator_matrix( transpose(e[l, cs]) * al)
   else
     trl = al - al
   end
   debug && println(trl)
-  dlr = indicator_matrix(I(length(l)) + transpose(trl) * trl)
+  dlr = StructuralCausalModels.indicator_matrix(I(length(l)) + transpose(trl) * trl)
   debug && println(dlr)
-  cl = transitive_closure(dlr)
+  cl = StructuralCausalModels.transitive_closure(dlr)
   debug && println(cl)
-  out = indicator_matrix(( al * cl * transpose(al)))
+  out = StructuralCausalModels.indicator_matrix(( al * cl * transpose(al)))
   debug && println(out)
   out = out[sel, sel]
   debug && println(out)
   debug && println(adjacency_matrix(out))
   adjacency_matrix(out)
 
 end
+
+export
+  induced_covariance_graph