Some examples of do-calculus

In [1]:

from IPython.display import display, Math

import pyagrum as gum
import pyagrum.lib.notebook as gnb

S. Tikka and J. Karvanen, 2016 [CRAN]

In [2]:

bn = gum.fastBN("w->x->z->y;w->z")

bn.cpt("w")[:] = [0.7, 0.3]

bn.cpt("x")[:] = [[0.4, 0.6], [0.3, 0.7]]

bn.cpt("z")[{"w": 0, "x": 0}] = [0.2, 0.8]
bn.cpt("z")[{"w": 0, "x": 1}] = [0.1, 0.9]
bn.cpt("z")[{"w": 1, "x": 0}] = [0.9, 0.1]
bn.cpt("z")[{"w": 1, "x": 1}] = [0.5, 0.5]

bn.cpt("y")[:] = [[0.1, 0.9], [0.8, 0.2]]
bn

Out[2]:

In [3]:

d = gum.CausalModel(bn, [("lat1", ["x", "y"])])

gnb.showCausalImpact(d, on="y", doing="x")

Causal Model

$$ \begin{equation*}P\left(y \mid \text{do}(x)\right) = \sum_{w,z}{P\left(z\mid w,x\right) \cdot P\left(w\right) \cdot \left(\sum_{x'}{P\left(y\mid z\right) \cdot P\left(x'\mid w\right)}\right)}\end{equation*} $$
Explanation : Identified via do-calculus (ID/IDC).

	y
x	0	1
0	0.5130	0.4870
1	0.6460	0.3540

Impact

Since we have the formula, let us compute by hand this intervention :

In [4]:

((bn.cpt("y") * bn.cpt("x")).sumOut(["x"]) * bn.cpt("w") * bn.cpt("z")).sumOut(["z", "w"]).putFirst("y")

Out[4]:

	y
x	0	1
0	0.5130	0.4870
1	0.6460	0.3540

In [5]:

bn = gum.fastBN("Z1->X->Z2->Y")
d = gum.CausalModel(bn, [("L1", ["Z1", "X"]), ("L2", ["Z1", "Z2"]), ("L3", ["Z1", "Y"]), ("L4", ["Y", "X"])], True)
gnb.showCausalImpact(d, "Y", "X", values={"X": 1})

Causal Model

Effect not identifiable
Explanation

No result
Impact

Front door

In [6]:

modele4 = gum.BayesNet()
modele4.add(gum.LabelizedVariable("Smoking"))
modele4.add(gum.LabelizedVariable("Cancer"))
modele4.add(gum.LabelizedVariable("Tar"))

modele4.addArc(0, 2)
modele4.addArc(2, 1)
modele4.addArc(0, 1)

# Smoking
modele4.cpt(0)[:] = [0.5, 0.5]

# Tar
modele4.cpt(2)[{"Smoking": 0}] = [0.4, 0.6]
modele4.cpt(2)[{"Smoking": 1}] = [0.3, 0.6]

# Cancer
modele4.cpt(1)[{"Smoking": 0, "Tar": 0}] = [0.1, 0.9]  # No Drug, Male -> healed in 0.8 of cases
modele4.cpt(1)[{"Smoking": 0, "Tar": 1}] = [0.15, 0.85]  # No Drug, Female -> healed in 0.4 of cases
modele4.cpt(1)[{"Smoking": 1, "Tar": 0}] = [0.2, 0.8]  # Drug, Male -> healed 0.7 of cases
modele4.cpt(1)[{"Smoking": 1, "Tar": 1}] = [0.25, 0.75]

d4 = gum.CausalModel(modele4, [("Genotype", ["Smoking", "Cancer"])], False)
gnb.showCausalModel(d4)

../_images/notebooks_64-Causality_DoCalculusExamples_10_0.svg

In [7]:

c = gum.CausalImpact(d4, on="Cancer", doing="Smoking")
if c.isIdentified():
  display(Math(c.toLatex()))

$\displaystyle P\left(Cancer \mid \text{do}(Smoking)\right) = \sum_{Tar}{P\left(Tar\mid Smoking\right) \cdot \left(\sum_{Smoking'}{P\left(Cancer\mid Smoking',Tar\right) \cdot P\left(Smoking'\right)}\right)}$

In [8]:

c.toDict()

Out[8]:

{'op': 'sum',
 'var': 'Tar',
 'term': {'op': '*',
  'op1': {'op': 'P', 'vars': ['Tar'], 'knowing': ['Smoking']},
  'op2': {'op': 'sum',
   'var': 'Smoking',
   'term': {'op': '*',
    'op1': {'op': 'P', 'vars': ['Cancer'], 'knowing': ['Tar', 'Smoking']},
    'op2': {'op': 'P', 'vars': ['Smoking']}}}}}

In [9]:

adjj = c.eval()

print(adjj)


      ║  Cancer           │
Smokin║0        │1        │
──────║─────────│─────────│
0     ║ 0.1800  │ 0.8200  │
1     ║ 0.1650  │ 0.7350  │

In [10]:

formula, adj, exp = gum.causalImpact(d4, on="Cancer", doing="Smoking", values={"Smoking": 0})

In [11]:

display(Math(formula.toLatex()))
adj

$\displaystyle P\left(Cancer \mid \text{do}(Smoking)\right) = \sum_{Tar}{P\left(Tar\mid Smoking\right) \cdot \left(\sum_{Smoking'}{P\left(Cancer\mid Smoking',Tar\right) \cdot P\left(Smoking'\right)}\right)}$

Out[11]:

Cancer
0	1
0.1800	0.8200

## Last example from R

In [12]:

m = gum.fastBN("z2->x->z1->y;z2->z1;z2->z3->y")

m.cpt("z2")[:] = [0.5, 0.5]
m.cpt("x")[:] = [
  [0.4, 0.6],  # z2=0
  [0.4, 0.6],
]  # z2=1
m.cpt("z3")[:] = [
  [0.3, 0.7],  # z2=0
  [0.3, 0.7],
]  # z2=1
m.cpt("z1")[{"z2": 0, "x": 0}] = [0.2, 0.8]
m.cpt("z1")[{"z2": 0, "x": 1}] = [0.25, 0.75]
m.cpt("z1")[{"z2": 1, "x": 0}] = [0.1, 0.9]
m.cpt("z1")[{"z2": 1, "x": 1}] = [0.15, 0.85]

m.cpt("y")[{"z1": 0, "z3": 0}] = [0.5, 0.5]
m.cpt("y")[{"z1": 0, "z3": 1}] = [0.45, 0.55]
m.cpt("y")[{"z1": 1, "z3": 0}] = [0.4, 0.6]
m.cpt("y")[{"z1": 1, "z3": 1}] = [0.35, 0.65]


d = gum.CausalModel(m, [("X-Z2", ["x", "z2"]), ("X-Z3", ["x", "z3"]), ("X-Y", ["x", "y"]), ("Y-Z2", ["y", "z2"])], True)

gnb.showCausalModel(d)

../_images/notebooks_64-Causality_DoCalculusExamples_17_0.svg

In [13]:

try:
  formula, result, msg = gum.causalImpact(d, on={"y", "z2", "z1", "z3"}, doing={"x"})
except gum.HedgeException as h:
  print(h.message)

print(msg)
display(Math(formula.toLatex()))

Identified via do-calculus (ID/IDC).

$\displaystyle P\left(y,z1,z2,z3 \mid \text{do}(x)\right) = P\left(z1\mid x,z2\right) \cdot P\left(z3\mid z2\right) \cdot \frac {\sum_{x'}{P\left(y\mid z1,z3\right) \cdot P\left(z3\mid z2\right) \cdot P\left(x'\mid z2\right) \cdot P\left(z2\right)}}{\sum_{x',y'}{P\left(y'\mid z1,z3\right) \cdot P\left(z3\mid z2\right) \cdot P\left(x'\mid z2\right) \cdot P\left(z2\right)}} \cdot P\left(z2\right)$

In [14]:

# computation for this formula directly in pyAgrum
f1 = m.cpt("x") * m.cpt("z2") * m.cpt("z3") * m.cpt("y")
f2 = f1.sumOut(["x"])
f3 = f1.sumOut(["x", "y"])
f4 = f2 / f3
pyResult = m.cpt("z3") * m.cpt("z1") * m.cpt("z2") * f4

According to [ref], the result should be :

\[P(z1|x,z2)\cdot \left( \sum_{x'} P(y,z3|z1,z2)\cdot P(x',z2) \right)\]

In [15]:

# computation for that formula
ie = gum.LazyPropagation(m)
refResult = (ie.evidenceJointImpact(["y", "z3"], ["z1", "z2"]) * ie.evidenceJointImpact(["x", "z2"], [])).sumOut(
  ["x"]
) * m.cpt("z1")

In [16]:

# the three tensors in the same variable context
r1 = gum.Tensor(result).fillWith(pyResult)
r2 = gum.Tensor(result).fillWith(refResult)
print(
  "Maximum error between these 3 versions : {}".format(
    max((r1 - r2).abs().max(), (r1 - result).abs().max(), (r2 - result).new_abs().max())
  )
)

Maximum error between these 3 versions : 5.551115123125783e-17

Unidentifiabilty

In [17]:

m1 = gum.fastBN("z1->x->z2->y")

cdg = gum.CausalModel(
  m1, [("Z1−X", ["z1", "x"]), ("Z1-Y", ["z1", "y"]), ("Z1-Z1", ["z1", "z2"]), ("X−Y", ["x", "y"])], True
)
gnb.showCausalModel(cdg)

../_images/notebooks_64-Causality_DoCalculusExamples_24_0.svg

In [18]:

err = gnb.showCausalImpact(cdg, "y", "x", values={"x": 0})

Causal Model

Effect not identifiable
Explanation

No result
Impact

another one

In [19]:

# EXEMPLE PAGE 17 : http://ftp.cs.ucla.edu/pub/stat_ser/r350.pdf

m1 = gum.fastBN("y<-z2->z3<-z1->x<-z3->y<-x")

gnb.showBN(m1)
d = gum.CausalModel(m1)

../_images/notebooks_64-Causality_DoCalculusExamples_27_0.svg

In [20]:

display(Math(gum.CausalImpact(d, on={"z1", "z2", "z3", "y"}, doing={"x"}).toLatex()))

$\displaystyle P\left(y,z1,z2,z3 \mid \text{do}(x)\right) = P\left(z3\mid z1,z2\right) \cdot P\left(y\mid x,z2,z3\right) \cdot P\left(z2\right) \cdot P\left(z1\right)$

In [21]:

display(Math(gum.causalImpact(d, on={"y"}, doing={"x"})[0].toLatex()))

$\displaystyle P\left(y \mid \text{do}(x)\right) = \sum_{z2,z3}{P\left(y\mid x,z2,z3\right) \cdot P\left(z2,z3\right)}$

In [22]:

display(Math(gum.CausalImpact(d, on={"z1", "z3", "y"}, doing={"x", "z2"}).toLatex()))

$\displaystyle P\left(y,z1,z3 \mid \text{do}(x,z2)\right) = P\left(z3\mid z1,z2\right) \cdot P\left(y\mid x,z2,z3\right) \cdot P\left(z1\right)$

In [23]:

display(Math(gum.CausalImpact(d, on={"y"}, doing={"x", "z2"}).toLatex()))

$\displaystyle P\left(y \mid \text{do}(x,z2)\right) = \sum_{z1,z3}{P\left(z3\mid z1,z2\right) \cdot P\left(y\mid x,z2,z3\right) \cdot P\left(z1\right)}$

Other example

In [24]:

# http://www.stats.ox.ac.uk/~lienart/gml15_causalinference.html

m1 = gum.fastBN("a->p->b->y<-p;a->y")
gnb.showBN(m1)
d = gum.CausalModel(m1)

../_images/notebooks_64-Causality_DoCalculusExamples_33_0.svg

In [25]:

display(Math(gum.CausalImpact(d, on={"y"}, doing={"a", "b"}).toLatex()))

$\displaystyle P\left(y \mid \text{do}(a,b)\right) = \sum_{p}{P\left(y\mid a,b,p\right) \cdot P\left(p\mid a\right)}$

example f

In [26]:

# https://cse.sc.edu/~mgv/talks/AIM2010.ppt , example (f)

m1 = gum.fastBN("X->Z1->Z2->Y<-Z1;X->Y")
d = gum.CausalModel(m1, [("l1", ["X", "Z1"]), ("l2", ["Y", "Z1"])], True)
gnb.showCausalModel(d)

../_images/notebooks_64-Causality_DoCalculusExamples_36_0.svg

In [27]:

display(Math(gum.CausalImpact(d, on={"Y"}, doing={"X"}).toLatex()))

$\displaystyle [unidentifiable]$

Example [Pearl,2009] Causality, p66

In [28]:

bn = gum.fastBN("Z1->Z2->Z3->Y<-X->Z2;Z2->Y;Z1->X->Z3<-Z1")
gnb.showBN(bn)

../_images/notebooks_64-Causality_DoCalculusExamples_39_0.svg

In [29]:

c = gum.CausalModel(bn, [("Z0", ("X", "Z1", "Z3"))], False)
gnb.showCausalModel(c)

../_images/notebooks_64-Causality_DoCalculusExamples_40_0.svg

In [30]:

formula, impact, explanation = gum.causalImpact(c, on="Y", doing="X")
gnb.showCausalImpact(c, "Y", "X")

Causal Model

$$ \begin{equation*}P\left(Y \mid \text{do}(X)\right) = \sum_{Z1,Z2,Z3}{\left(\sum_{X'}{P\left(Z3\mid Z2\right) \cdot P\left(X'\right) \cdot P\left(Z1\right)}\right) \cdot P\left(Z2\mid X,Z1\right) \cdot P\left(Y\mid X,Z2,Z3\right)}\end{equation*} $$
Explanation : Identified via do-calculus (ID/IDC).

	X
Y	0	1
0	0.3786	0.4964
1	0.6214	0.5036

Impact

Causality without numerical computations

You may want to compute some impact without having a (discrete) quantifications of the model. In that cas, you are not intersted in causalImpact which will include some computations of a numerical (discrete) impact.

In [31]:

f = gum.CausalImpact(c, on="Y", doing="X")
print("- Version latex")
display(Math(f.toLatex()))
print("- Version Abstract Syntax Tree (in a dict)")
f.toDict()

- Version latex

$\displaystyle P\left(Y \mid \text{do}(X)\right) = \sum_{Z1,Z2,Z3}{\left(\sum_{X'}{P\left(Z3\mid Z2\right) \cdot P\left(X'\right) \cdot P\left(Z1\right)}\right) \cdot P\left(Z2\mid X,Z1\right) \cdot P\left(Y\mid X,Z2,Z3\right)}$

- Version Abstract Syntax Tree (in a dict)

Out[31]:

{'op': 'sum',
 'var': 'Z2',
 'term': {'op': 'sum',
  'var': 'Z3',
  'term': {'op': 'sum',
   'var': 'Z1',
   'term': {'op': '*',
    'op1': {'op': 'sum',
     'var': 'X',
     'term': {'op': '*',
      'op1': {'op': 'P', 'vars': ['Z3'], 'knowing': ['Z2']},
      'op2': {'op': '*',
       'op1': {'op': 'P', 'vars': ['X'], 'knowing': []},
       'op2': {'op': 'P', 'vars': ['Z1']}}}},
    'op2': {'op': '*',
     'op1': {'op': 'P', 'vars': ['Z2'], 'knowing': ['Z1', 'X']},
     'op2': {'op': 'P', 'vars': ['Y'], 'knowing': ['Z3', 'Z2', 'X']}}}}}}

In [32]:

f.print_ast()

┌────────────┐
│ P(Y|do(X)) │
└────────────┘
      └─ sum
         ├─ var: Z2
         └─ sum
            ├─ var: Z3
            └─ sum
               ├─ var: Z1
               └─ *
                  ├─ sum
                  │  ├─ var: X
                  │  └─ *
                  │     ├─ P
                  │     │  ├─ vars: ['Z3']
                  │     │  └─ knowing: ['Z2']
                  │     └─ *
                  │        ├─ P
                  │        │  ├─ vars: ['X']
                  │        │  └─ knowing: []
                  │        └─ P
                  │           └─ vars: ['Z1']
                  └─ *
                     ├─ P
                     │  ├─ vars: ['Z2']
                     │  └─ knowing: ['Z1', 'X']
                     └─ P
                        ├─ vars: ['Y']
                        └─ knowing: ['Z3', 'Z2', 'X']

In [33]:

m1 = gum.CausalModel(gum.fastBN("y<-z2->z3<-z1->x<-z3->y<-x"))
f1 = gum.CausalImpact(m1, on={"z1", "y"}, doing={"x", "z2"}, knowing={"z3"})
f1.print_ast()

┌─────────────────────┐
│ P(y,z1|do(x,z2),z3) │
└─────────────────────┘
          └─ /
             ├─ *
             │  ├─ P
             │  │  ├─ vars: ['z3']
             │  │  └─ knowing: ['z1', 'z2']
             │  └─ *
             │     ├─ P
             │     │  ├─ vars: ['y']
             │     │  └─ knowing: ['z3', 'z2', 'x']
             │     └─ P
             │        └─ vars: ['z1']
             └─ sum
                ├─ var: z1
                └─ *
                   ├─ P
                   │  ├─ vars: ['z3']
                   │  └─ knowing: ['z1', 'z2']
                   └─ P
                      └─ vars: ['z1']

In [ ]: