Esercizio - Algoritmo PC sul grafo con cinque variabili (slide del corso)

Testo (slide del corso di Intelligenza Artificiale, Ing. dell'Informazione UniPD, Causal Inference parte 4, slide 9-14). Si ha un insieme di dati da un processo sconosciuto con cinque variabili $A,B,C,D,E$ e si vuole scoprirne la struttura causale $G^*$. Il grafo vero (nascosto all'algoritmo) è

1. Descrivere i quattro passi dell'algoritmo PC applicati a questo esempio, indicando gli insiemi separatori usati per togliere gli archi, quali collider si trovano e come si orientano gli altri archi.
2. Implementare PC con il test di Fisher $z$ e verificare con dati simulati da $G^*$ che si ritrova il grafo.
3. Cosa restituisce PC se il processo è una catena $A\to B\to C$? E un collider $A\to B\leftarrow C$?

Teoria: Controfattuali e scoperta causaleUn controfattuale chiede «cosa sarebbe successo a questo individuo se X fosse stato x?»: $Y_x(u)$ è il valore di Y nel modello $M_x$ (equazione di X sostituita da X=x) con le esogene fissate al valore $u$ dell'individuo. Si calcola in tre passi: abduzione (aggiornare $P(U)$ con l'evidenza, nel caso deterministico ricavare $u$ dai dati), azione (sostituire l'equazione di X con la costante), predizione (calcolare Y in $M_x$). Differenza con gli interventi: $P(Y\mid do(x))\equiv P(Y_x)$ è a livello di popolazione, ma in generale $P(Y\mid do(x),z)\ne P(Y_x\mid z)$ perché nel secondo il condizionamento è prima dell'intervento. Nei modelli lineari gaussiani i coefficienti di regressione di ogni variabile sui suoi genitori sono i coefficienti strutturali, e l'effetto causale totale è la somma dei prodotti dei coefficienti lungo i cammini non backdoor. Scoperta causale (algoritmo PC): si parte dal grafo completo non orientato, si toglie ogni arco X−Y con $X\perp!!!\perp Y\mid Z$ a $|Z|$ crescente (test di Fisher), si orientano i collider (Z mai nell'insieme separatore), poi si propaga l'orientamento; catena e fork non sono distinguibili (stessa classe di equivalenza di Markov). Assunzioni: aciclicità, sufficienza, Markov, fedeltà. Argomento del corso da 9 CFU.Controfattuali e scoperta causale → (assunzioni di PC, classi di equivalenza di Markov); d-separazione in Inferenza causale - modelli causali strutturali e scala della causalitàI dati e la distribuzione congiunta da soli non dicono quali variabili causano quali: lo mostrano il paradosso di Simpson (la medicina funziona in ogni sottogruppo ma sembra peggiorare la popolazione) e il fatto che reti bayesiane diverse (A→B e B→A) descrivono la stessa congiunta. Serve un modello causale: un modello causale strutturale (SCM) ha variabili esogene U, endogene V e funzioni $v_i=f_i(\text{genitori},u_i)$; il suo grafo è la rete causale, e una rete causale è una rete bayesiana in cui le frecce sono cause. Le tre strutture elementari decidono l'indipendenza: catena $A\to B\to C$ e fork $A\leftarrow B\to C$ si bloccano osservando B, il collider $A\to B\leftarrow C$ si blocca se B e i suoi discendenti NON sono osservati (osservarlo crea dipendenza, explaining away). d-separazione: X e Y sono separati da Z se ogni cammino è bloccato; allora $X\perp!!!\perp Y\mid Z$. La scala della causalità (Pearl): associazione $P(y\mid x)$, intervento $P(y\mid do(x))$, controfattuali. Argomento del corso da 9 CFU.Inferenza causale - modelli causali strutturali e scala della causalità →.

Soluzione

1. I quattro passi (a mano, con le indipendenze del grafo vero)

Ricordo le assunzioni: aciclicità, sufficienza (nessun confondente non osservato), Markov e fedeltà: grazie alla fedeltà le indipendenze osservate nei dati sono tutte e sole le d-separazioni di $G^*$. Quindi nel grafo vero ($A\to B\leftarrow C$, $B\to D$, $B\to E$) i test valgono:

• $A\perp\!\!\!\perp C$ non condizionato (collider $B$ non osservato), ma $A\not\perp\!\!\!\perp C\mid B$;
• $A\perp\!\!\!\perp D\mid B$, $A\perp\!\!\!\perp E\mid B$, $C\perp\!\!\!\perp D\mid B$, $C\perp\!\!\!\perp E\mid B$, $D\perp\!\!\!\perp E\mid B$ (catene e fork bloccati da $B$);
• gli archi $A-B$, $B-C$, $B-D$, $B-E$ non si rendono mai indipendenti con nessun $Z$ (sono adiacenze vere).

Passo 1: grafo completo non orientato. Dieci archi (tutte le coppie).

Passo 2: scheletro. Si prova $Z$ vuoto, poi di una variabile, ecc.

• $|Z|=0$: $A\perp\!\!\!\perp C\mid\{\}$: si toglie $A-C$ e si ricorda l'insieme separatore $\{\}$.
• $|Z|=1$: si tolgono $A-D$ e $A-E$ ($Z=\{B\}$), $C-D$ e $C-E$ ($Z=\{B\}$), $D-E$ ($Z=\{B\}$); rimangono $A-B$, $B-C$, $B-D$, $B-E$.
• $|Z|\ge2$: nessun altro arco si toglie. Lo scheletro è la «X» $A-B-C$, $B-D$, $B-E$.

Passo 3: collider. Per ogni cammino $X-Z-Y$ con $X$ e $Y$ non adiacenti, si orienta $X\to Z\leftarrow Y$ se $Z$ non era nell'insieme separatore di $X,Y$.

• $A-B-C$: $A$ e $C$ non adiacenti, insieme separatore $\{\}$, $B\notin\{\}$: collider $A\to B\leftarrow C$.
• $A-B-D$, $A-B-E$, $C-B-D$, $C-B-E$, $D-B-E$: l'insieme separatore è $\{B\}$, cioè $B$ è dentro: non sono collider; sono catene o fork, la classe di equivalenza di Markov non distingue.

Passo 4: propagazione. $A\to B-D$ con $A$ e $D$ non adiacenti: orientando $D\to B$ si creerebbe il collider $A\to B\leftarrow D$, già escluso al passo 3, quindi $B\to D$. Allo stesso modo $A\to B-E$ dà $B\to E$.

Risultato: $A\to B\leftarrow C$, $B\to D$, $B\to E$, esattamente $G^*$: qui la classe di equivalenza di Markov contiene un solo grafo, perché il collider fissa tutto il resto.

2. Implementazione con il test di Fisher $z$

Il test di indipendenza condizionata per variabili gaussiane usa la correlazione parziale $r$ di $X$ e $Y$ dato $Z$ (dall'inversa della matrice di correlazione) e la statistica $$z=\tfrac12\ln\frac{1+r}{1-r}\sqrt{n-|Z|-3},\qquad p=2\big(1-\Phi(|z|)\big);$$ si accetta l'indipendenza se $p>\alpha$ (qui $\alpha=0.01$), e $z$ è circa normale standard sotto l'ipotesi di indipendenza. I dati sono simulati da $G^*$ con un SCM lineare gaussiano (coefficienti scelti qui: $B=0.9A+0.8C+0.7\varepsilon$, $D=0.9B+0.7\varepsilon$, $E=0.8B+0.7\varepsilon$, $n=5000$). SciPy serve solo per $\Phi$.

import itertools
import numpy as np
from scipy.stats import norm

def indipendenti(C, n, i, j, S, alpha=0.01):
    """test di Fisher z: i e j indipendenti dato S? (C = matrice di correlazione dei dati)"""
    idx = [i, j] + list(S)
    P = np.linalg.inv(C[np.ix_(idx, idx)])
    r = -P[0, 1] / np.sqrt(P[0, 0] * P[1, 1])                          # correlazione parziale
    z = 0.5 * np.log((1 + r) / (1 - r)) * np.sqrt(n - len(S) - 3)
    p = 2 * (1 - norm.cdf(abs(z)))
    return p > alpha, r, p

def pc(dati, nomi, alpha=0.01, stampa=True):
    n, k = dati.shape
    C = np.corrcoef(dati.T)
    adj = {i: set(range(k)) - {i} for i in range(k)}                    # 1. grafo completo non orientato
    sep = {}                                                            # insiemi separatori
    L = 0
    while any(len(adj[i]) - 1 >= L for i in range(k)):                  # 2. scheletro, |Z| = L crescente
        for i, j in itertools.combinations(range(k), 2):
            if j not in adj[i]: continue
            for S in itertools.combinations(sorted(adj[i] - {j}), L):
                ind, r, p = indipendenti(C, n, i, j, S, alpha)
                if ind:
                    adj[i].discard(j); adj[j].discard(i); sep[i, j] = sep[j, i] = set(S)
                    if stampa: print(f"L={L}: {nomi[i]} _|_ {nomi[j]} | {{{','.join(nomi[s] for s in S)}}}  (r = {r:+.3f}, p = {p:.2f}) -> arco tolto")
                    break
        L += 1
    orient = set()                                                      # coppie (x, z) = arco x -> z
    for z in range(k):                                                  # 3. collider
        for x, y in itertools.combinations(sorted(adj[z]), 2):
            if y not in adj[x] and z not in sep[x, y]: orient |= {(x, z), (y, z)}
    cambiato = True
    while cambiato:                                                     # 4. propagazione: x -> z - y, x e y non adiacenti
        cambiato = False
        for x, z in list(orient):
            for y in adj[z]:
                if (z, y) not in orient and (y, z) not in orient and y != x and y not in adj[x]:
                    orient.add((z, y)); cambiato = True
    non_orient = {tuple(sorted((i, j))) for i in adj for j in adj[i] if (i, j) not in orient and (j, i) not in orient}
    return orient, non_orient

def descrivi(orient, non_orient, nomi):
    return sorted(f"{nomi[x]}->{nomi[z]}" for x, z in orient) + sorted(f"{nomi[i]}-{nomi[j]}" for i, j in non_orient)

# dati da G*: A -> B <- C, B -> D, B -> E  (SCM lineare gaussiano; coefficienti scelti qui)
rng = np.random.default_rng(3); N = 5000
A, C = rng.normal(size=N), rng.normal(size=N)
B = 0.9 * A + 0.8 * C + 0.7 * rng.normal(size=N)
D = 0.9 * B + 0.7 * rng.normal(size=N)
E = 0.8 * B + 0.7 * rng.normal(size=N)
nomi = list("ABCDE")
o, n_o = pc(np.column_stack([A, B, C, D, E]), nomi)
print("risultato:", descrivi(o, n_o, nomi))

# controllo: una catena A -> B -> C non e' orientabile; un collider si
A = rng.normal(size=N); B = A + 0.7 * rng.normal(size=N); C = B + 0.7 * rng.normal(size=N)
print("catena   :", descrivi(*pc(np.column_stack([A, B, C]), list("ABC"), stampa=False), list("ABC")))
A, B = rng.normal(size=N), rng.normal(size=N); C = A + B + 0.7 * rng.normal(size=N)
print("collider :", descrivi(*pc(np.column_stack([A, B, C]), list("ABC"), stampa=False), list("ABC")))

Output:

L=0: A _|_ C | {}  (r = +0.015, p = 0.30) -> arco tolto
L=1: A _|_ D | {B}  (r = -0.007, p = 0.61) -> arco tolto
L=1: A _|_ E | {B}  (r = +0.001, p = 0.97) -> arco tolto
L=1: C _|_ D | {B}  (r = +0.005, p = 0.73) -> arco tolto
L=1: C _|_ E | {B}  (r = +0.013, p = 0.37) -> arco tolto
L=1: D _|_ E | {B}  (r = -0.005, p = 0.75) -> arco tolto
risultato: ['A->B', 'B->D', 'B->E', 'C->B']
catena   : ['A-B', 'B-C']
collider : ['A->C', 'B->C']

Le sei rimozioni coincidono con quelle a mano (e con le slide). Le coppie non rimosse hanno correlazioni parziali grandi con $p\approx0$ (per esempio $A\perp\!\!\!\perp B\mid D$ è rifiutata, $r=0.39$, perché $A$ e $B$ sono adiacenti).

3. Catena e collider

• Catena $A\to B\to C$. Il test dà $A\perp\!\!\!\perp C\mid B$: si toglie $A-C$ con insieme separatore $\{B\}$, che contiene $B$: nessun collider, nessuna propagazione. Risultato: ['A-B', 'B-C'], due archi non orientati: la catena, il fork opposto $A\leftarrow B\to C$ e $A\leftarrow B\leftarrow C$ sono nella stessa classe di equivalenza di Markov e dai soli dati non si distinguono.
• Collider $A\to B\leftarrow C$. Il test dà $A\perp\!\!\!\perp C\mid\{\}$: insieme separatore vuoto, $B$ non c'è: collider, orientamento completo. Risultato: entrambi gli archi sono orientati verso il nodo comune. Nel codice di controllo le due cause sono generate come $A$ e $B$ e l'effetto comune è la variabile chiamata C, quindi l'uscita è ['A->C', 'B->C'].