Examples

Quick start#

Open in Colab Download this notebook

You call every method in Andrey the same way: pass it your data, and read the result it returns. This notebook loads the protein measurements of Sachs et al. (2005), learns a graph from them with four classic methods, PC, GES, FCI, and DirectLiNGAM, reads each result, and draws each graph.

Start#

Choose a CPU runtime and run this cell first. In Colab it installs andrey-core==0.17.1 from PyPI; local kernels use the installed environment.

import os
import sys

IN_COLAB = "google.colab" in sys.modules

if IN_COLAB:
    os.environ["ANDREY_DEVICE"] = "cpu"
    !pip install --quiet "andrey-core[viz]==0.17.1" pandas
import andrey
import andrey.viz
from andrey.data import load_dataset

The data#

Sachs and colleagues (2005) measured 11 proteins and phospholipids in 853 single human immune cells. Biologists had already mapped how these molecules signal to one another, so a graph learned from the measurements can be checked against a known network of 17 links. load_dataset downloads the measurements on first use and caches them.

Each row is one cell and each column one molecule. For your own data, pass a pandas DataFrame or a NumPy array of the same shape: one row per sample, one column per variable.

sachs = load_dataset("sachs")  # 853 cells, 11 proteins, and the known network
data = sachs.data
data.head()
raf mek plc pip2 pip3 erk akt pka pkc p38 jnk
0 26.4 13.2 8.82 18.30 58.80 6.61 17.0 414.0 17.00 44.9 40.0
1 35.9 16.5 12.30 16.80 8.13 18.60 32.5 352.0 3.37 16.5 61.5
2 59.4 44.1 14.60 10.20 13.00 14.90 32.5 403.0 11.40 31.9 19.5
3 73.0 82.8 23.10 13.50 1.29 5.83 11.8 528.0 13.70 28.6 23.1
4 33.7 19.8 5.19 9.73 24.80 21.10 46.1 305.0 4.66 25.7 81.3

In math#

Each column of data is a variable XjX_j, here d=11d = 11 of them, and each row is one draw from a distribution PP over V={X1,…,Xd}V = \{X_1, \dots, X_d\}. We model each variable as a function of its parents in a directed acyclic graph (DAG) GG, its direct causes, and of its own noise:

Xj=fj(XpaG(j),Nj),j=1,…,d,X_j = f_j(X_{\mathrm{pa}_G(j)}, N_j), \qquad j = 1, \dots, d,

with the noise terms NjN_j jointly independent, so no hidden variable causes two columns. A method estimates GG, or the part of GG the data identify: PC and GES below return a class of DAGs, FCI drops the independence of the noise to allow hidden common causes, and DirectLiNGAM identifies GG itself by assuming linear fjf_j and non-Gaussian noise. The concepts page states each assumption and what each graph claims.

PC#

PC starts with every pair of molecules linked. It removes a link when a statistical test finds the two independent, either alone or given some of their other neighbors, and orients a link only when a pattern in the tests forces its direction.

pc = andrey.pc(data)
print(pc)
pc  cpdag  |  11 nodes  |  8 edges (6 undirected, 2 directed)

  raf -- mek
  plc -- pip3
  pip2 -- pip3
  erk -- akt
  erk -- pka
  akt -- pka
  p38 -> pkc
  jnk -> pkc

  indep_test: fisherz
  tie_break: rank

The first line names the method and the kind of graph. This one is a CPDAG: it stands for every causal graph that independence tests cannot tell apart. A directed edge (->) points the same way in all of them, and an undirected edge (--) can point either way.

You can read the same parts in code. pc.structure is the graph: its kind says how to read the edges, and oriented_edges() lists each edge as (source, target, type), with your column names. pc.metadata records the settings of the run.

graph = pc.structure
print("kind:", graph.kind)
print("directed:", [edge for edge in graph.oriented_edges() if edge[2] == "directed"])
print("test:", pc.metadata["indep_test"])
kind: cpdag
directed: [('p38', 'pkc', 'directed'), ('jnk', 'pkc', 'directed')]
test: fisherz

andrey.viz.draw renders the graph as an SVG, which the notebook shows. A directed edge ends in an arrowhead; an undirected edge has none and is drawn in amber.

andrey.viz.draw(pc.structure)
rafmekplcpip2pip3erkaktpkapkcp38jnk

GES#

GES, greedy equivalence search, reaches a graph a different way: it searches for the graph with the best score, instead of testing independence. You call it the same way.

ges = andrey.ges(data)
print(ges)
ges  cpdag  |  11 nodes  |  8 edges (6 undirected, 2 directed)

  raf -- mek
  plc -- pip3
  pip2 -- pip3
  erk -- akt
  erk -- pka
  akt -- pka
  p38 -> pkc
  jnk -> pkc

  score: 6.692e+04
  score_func: bic
  penalty_discount: 1
  ebic_gamma: 1
  tie_break: rank

GES returns a CPDAG too, and on these data it is the same graph as PC’s. Its metadata holds the score of the graph, where lower is better, and the settings that produced it.

print("same edges as PC:", ges.structure.oriented_edges() == pc.structure.oriented_edges())
andrey.viz.draw(ges.structure)
same edges as PC: True
rafmekplcpip2pip3erkaktpkapkcp38jnk

FCI#

PC and GES assume that every common cause of two measured variables is measured too. FCI does not: it allows hidden common causes, and says less about each edge in return.

fci = andrey.fci(data)
print(fci)
fci  pag  |  11 nodes  |  8 edges (6 circle, 2 partially_directed)

  raf o-o mek
  plc o-o pip3
  pip2 o-o pip3
  erk o-o akt
  erk o-o pka
  akt o-o pka
  p38 o-> pkc
  jnk o-> pkc

  indep_test: fisherz
  tie_break: rank

FCI returns a PAG. A circle (o) at an end of an edge means the tests do not say whether that end is an arrowhead or a tail. So p38 o-> pkc says that pkc does not cause p38, and leaves open whether p38 causes pkc or a hidden variable causes both. The drawing shows each circle as a hollow circle.

andrey.viz.draw(fci.structure)
rafmekplcpip2pip3erkaktpkapkcp38jnk

DirectLiNGAM#

DirectLiNGAM assumes that each variable is a linear function of its causes plus noise that is not Gaussian. When that holds, the data identify one graph, so it returns a DAG: every edge has a direction.

lingam = andrey.direct_lingam(data)
print(lingam)
direct_lingam  dag  |  11 nodes  |  11 edges (11 directed)

  mek -> raf  (+1.2)
  mek -> pip2  (+0.285)
  plc -> pip2  (+0.583)
  plc -> pip3  (+0.236)
  pip2 -> pip3  (+0.0952)
  erk -> akt  (+1.37)
  erk -> pka  (+1.84)
  pka -> akt  (+0.017)
  pkc -> p38  (+1.28)
  jnk -> pkc  (-0.0526)
  jnk -> p38  (+0.0582)

  order: erk -> plc -> mek -> pip2 -> pip3 -> pka -> jnk -> akt -> pkc -> p38 -> raf

  tie_break: rank

The number beside an edge is its weight, the coefficient of the cause in the linear equation for the effect, and order is the causal order the method found. These directions rest on the method’s assumption, so treat them as hypotheses to test.

andrey.viz.draw(lingam.structure)
rafmekplcpip2pip3erkaktpkapkcp38jnk

Next#

The Methods notebooks each take one family of methods further: its other methods, its settings, and the data it suits.