Esercizio - Aggiustamento backdoor con un confondente non osservato e con un collider (slide del corso)

Testo (esempi 1 e 2 delle slide Causal Inference parte 2, slide 13 e 14; le slide danno i grafi e le formule, i numeri sono scelti qui). Variabili booleane $X,Y,Z,W$. Si vuole l'effetto causale di $X$ su $Y$.

1. Esempio 1. Grafo $Z\to X$, $Z\to W$, $W\to Y$, $X\to Y$, con $Z$ non osservato. Probabilità: $P(z)=0.4$; $P(w\mid z)=0.8$, $P(w\mid\neg z)=0.3$; $P(x\mid z)=0.8$, $P(x\mid\neg z)=0.2$; $P(y\mid x,w)=0.7$ per $(x,w)$, $0.3$ per $(x,\neg w)$, $0.5$ per $(\neg x,w)$, $0.1$ per $(\neg x,\neg w)$.
   • a. Trovare i cammini backdoor e gli insiemi validi.
   • b. Calcolare $P(y\mid do(x))$ e $P(y\mid do(\neg x))$ con il solo $W$ osservato (la formula $\sum_wP(y\mid x,w)P(w)$) usando solo le probabilità osservabili di $(X,W,Y)$.
   • c. Confrontare con $P(y\mid x)$ e con il valore vero ottenuto dal modello strutturale (chirurgia), e controllare con una simulazione.
   • d. Cosa succede se si aggiusta invece su $Z$ o su $\{Z,W\}$ (verifica del modello)?
2. Esempio 2. Grafo $Z\to X$, $Z\to W$, $Y\to W$, $X\to Y$ (collider $W$). Probabilità: $P(z)=0.5$; $P(x\mid z)=0.8$, $P(x\mid\neg z)=0.2$; $P(y\mid x)=0.4$, $P(y\mid\neg x)=0.7$; $P(w\mid z,y)=0.9$, $P(w\mid z,\neg y)=0.4$, $P(w\mid\neg z,y)=0.5$, $P(w\mid\neg z,\neg y)=0.1$. Quali insiemi si possono usare? Cosa succede aggiustando su $W$?

Teoria: Interventi e operatore do - criterio backdoor e aggiustamentoOsservare $X=x$ e forzare $X=x$ sono cose diverse: $P(y\mid x)\ne P(y\mid do(x))$ quando c'è un confondente. L'intervento $do(X=x)$ si calcola con la «chirurgia» del grafo (si tolgono le frecce entranti in X) e si traduce in probabilità osservazionali con la formula di aggiustamento $P(y\mid do(x))=\sum_z P(y\mid x,z)P(z)$, dove $Z$ sono i genitori di X, o più in generale un insieme S che soddisfa il criterio backdoor: blocca ogni cammino backdoor (cammino da X a Y che inizia con una freccia entrante in X) e non contiene discendenti di X. Effetto causale medio $ACE=P(y\mid do(x))-P(y\mid do(\neg x))$. Se il confondente non è osservabile ma c'è un mediatore Z con X→Z→Y, vale il front-door $P(y\mid do(x))=\sum_z P(z\mid x)\sum_{x'}P(y\mid x',z)P(x')$. Mediazione: effetto diretto controllato CDE con $do(x),do(z)$. Se due insiemi backdoor danno risultati diversi sui dati, il modello è sbagliato. Argomento del corso da 9 CFU.Interventi e operatore do - criterio backdoor e aggiustamento → (criterio e verifica del modello).

Soluzione

Esempio 1

a. Un cammino backdoor: $X\leftarrow Z\to W\to Y$ (fork in $Z$, poi catena). Per il criterio sono validi $\{Z\}$, $\{W\}$ e $\{Z,W\}$ (nessuno è discendente di $X$). Con $Z$ non osservato resta $\{W\}$.

b. Si usano solo i dati osservabili di $(X,W,Y)$. La distribuzione $P(x,w)$ si ricava marginalizzando $Z$: per esempio $P(x,w)=\sum_zP(z)P(w\mid z)P(x\mid z)$:

Da cui $P(w)=0.292+0.208=0.5$, $P(\neg w)=0.5$ e $P(x)=0.44$. Poiché $Y$ dipende da $X$ e $W$ soltanto, $P(y\mid x,w)$ nei dati coincide con la tabella: $0.7,0.3,0.5,0.1$. Dunque $$P(y\mid do(x))=P(y\mid x,w)P(w)+P(y\mid x,\neg w)P(\neg w)=0.7\cdot0.5+0.3\cdot0.5=0.5,$$ $$P(y\mid do(\neg x))=0.5\cdot0.5+0.1\cdot0.5=0.3,\qquad ACE=0.5-0.3=0.2.$$

c. Il valore osservazionale è diverso: $P(w\mid x)=0.292/0.44=0.6636$ e $P(w\mid\neg x)=0.208/0.56=0.3714$, quindi $$P(y\mid x)=0.3\cdot0.3364+0.7\cdot0.6636=0.5655,\qquad P(y\mid\neg x)=0.1\cdot0.6286+0.5\cdot0.3714=0.2486,$$ con un «effetto» osservato di $0.317$ contro quello causale $0.2$: chi ha $X=x$ ha più spesso $W=w$ (effetto di $Z$) e $W$ aiuta $Y$, un'influenza spuria. Il valore vero si ottiene dalla chirurgia (si toglie la freccia $Z\to X$ e si fissa $X$): $P(y\mid do(x))=\sum_{z,w}P(z)P(w\mid z)P(y\mid x,w)=0.5$ e $0.3$, uguali all'aggiustamento.

d. Anche con $\{Z\}$ e con $\{Z,W\}$ la formula dà $0.5$ e $0.3$: modello coerente (più insiemi validi, stesso risultato). Se con dati reali i risultati fossero diversi, il modello causale sarebbe sbagliato.

Esempio 2

Cammino backdoor: $X\leftarrow Z\to W\leftarrow Y$, con collider $W$. Insiemi validi: $\{\}$ (il collider $W$ blocca) e $\{Z\}$ (blocca il fork $X\leftarrow Z\to W$). Non sono validi $\{W\}$ e $\{Z,W\}$: $W$ è un discendente di $X$ ($X\to Y\to W$), quindi viola la seconda condizione del criterio, e in più condizionare sul collider apre il cammino. Con $\{\}$: togliere la freccia $Z\to X$ non influenza $Y$, quindi $P(y\mid do(x))=P(y\mid x)$: $0.7$ e $0.4$ (si legge dalle tabelle: $Y$ dipende solo da $X$).

Se invece si aggiusta su $W$, errore: $\sum_wP(y\mid x,w)P(w)$ dà $0.7139$ per $\neg x$ e $0.382$ per $x$ (invece di $0.7$ e $0.4$): condizionando su $W$ si apre il cammino $X\leftarrow Z\to W\leftarrow Y$ e si introduce una dipendenza spuria tra $X$ e $Y$. Aggiustando su $Z$ si ottengono $0.7$ e $0.4$.

Codice

import itertools
import numpy as np

# ---------- Esempio 1: Z -> X, Z -> W, W -> Y, X -> Y   (Z non osservato)
pz = {1: 0.4, 0: 0.6}                                        # P(Z = z)
pw = {0: 0.3, 1: 0.8}                                        # P(W = 1 | z)
px = {0: 0.2, 1: 0.8}                                        # P(X = 1 | z)
py = {(0, 0): 0.1, (0, 1): 0.5, (1, 0): 0.3, (1, 1): 0.7}    # P(Y = 1 | x, w)
b = lambda v, p: p if v else 1 - p                           # P(V = v) quando P(V = 1) = p

# congiunta completa P(z, w, x, y)
J = {(z, w, x, y): pz[z] * b(w, pw[z]) * b(x, px[z]) * b(y, py[x, w])
     for z, w, x, y in itertools.product((0, 1), repeat=4)}
P = lambda **k: sum(p for (z, w, x, y), p in J.items() if all(dict(z=z, w=w, x=x, y=y)[n] == v for n, v in k.items()))

# dati osservabili: solo (X, W, Y); Z non si usa nelle formule che seguono
print("P(x,w):", {(x, w): round(P(x=x, w=w), 3) for x in (0, 1) for w in (0, 1)})
oss = {x: P(x=x, y=1) / P(x=x) for x in (0, 1)}
adj_w = {x: sum(P(x=x, w=w, y=1) / P(x=x, w=w) * P(w=w) for w in (0, 1)) for x in (0, 1)}
adj_z = {x: sum(P(x=x, z=z, y=1) / P(x=x, z=z) * P(z=z) for z in (0, 1)) for x in (0, 1)}
vero = {x: sum(pz[z] * b(w, pw[z]) * py[x, w] for z in (0, 1) for w in (0, 1)) for x in (0, 1)}   # chirurgia: manca P(x | z)
for nome, d in [("osservazionale P(y|x)", oss), ("aggiustamento su W", adj_w), ("aggiustamento su Z", adj_z), ("verita (chirurgia)", vero)]:
    print(f"{nome:24s} x=1: {d[1]:.4f}  x=0: {d[0]:.4f}  ACE = {d[1] - d[0]:.4f}")

# simulazione dell'SCM e stima dell'aggiustamento su W da un campione (n = 2 milioni)
rng = np.random.default_rng(5); N = 2_000_000
Z = rng.random(N) < 0.4
W = rng.random(N) < np.where(Z, 0.8, 0.3)
Ux, Uy = rng.random(N), rng.random(N)
def simula(x=None):
    X = (Ux < np.where(Z, 0.8, 0.2)) if x is None else np.full(N, bool(x))
    return X, Uy < np.select([~X & ~W, ~X & W, X & ~W, X & W], [0.1, 0.5, 0.3, 0.7])
print("simulazione do(x=1), do(x=0):", simula(1)[1].mean().round(4), simula(0)[1].mean().round(4))
X, Y = simula()
stima = lambda x: sum(Y[(X == x) & (W == w)].mean() * (W == w).mean() for w in (0, 1))
print("stima da campione (solo X, W, Y):", stima(True).round(4), stima(False).round(4))

# ---------- Esempio 2: Z -> X, Z -> W, Y -> W, X -> Y   (W e' un collider)
pz2 = {0: 0.5, 1: 0.5}; px2 = {0: 0.2, 1: 0.8}               # P(X = 1 | z)
py2 = {0: 0.7, 1: 0.4}                                       # P(Y = 1 | x)
pw2 = {(0, 0): 0.1, (0, 1): 0.5, (1, 0): 0.4, (1, 1): 0.9}   # P(W = 1 | z, y)
J2 = {(z, x, y, w): pz2[z] * b(x, px2[z]) * b(y, py2[x]) * b(w, pw2[z, y])
      for z, x, y, w in itertools.product((0, 1), repeat=4)}
P2 = lambda **k: sum(p for (z, x, y, w), p in J2.items() if all(dict(z=z, x=x, y=y, w=w)[n] == v for n, v in k.items()))
print("--- esempio 2")
print("P(y|x) con insieme vuoto:", [round(P2(x=x, y=1) / P2(x=x), 4) for x in (0, 1)], " verita: ", [py2[0], py2[1]])
print("aggiustando su W (errato):", [round(sum(P2(x=x, w=w, y=1) / P2(x=x, w=w) * P2(w=w) for w in (0, 1)), 4) for x in (0, 1)])
print("aggiustando su Z:         ", [round(sum(P2(x=x, z=z, y=1) / P2(x=x, z=z) * P2(z=z) for z in (0, 1)), 4) for x in (0, 1)])

Output:

P(x,w): {(0, 0): 0.352, (0, 1): 0.208, (1, 0): 0.148, (1, 1): 0.292}
osservazionale P(y|x)    x=1: 0.5655  x=0: 0.2486  ACE = 0.3169
aggiustamento su W       x=1: 0.5000  x=0: 0.3000  ACE = 0.2000
aggiustamento su Z       x=1: 0.5000  x=0: 0.3000  ACE = 0.2000
verita (chirurgia)       x=1: 0.5000  x=0: 0.3000  ACE = 0.2000
simulazione do(x=1), do(x=0): 0.5001 0.3
stima da campione (solo X, W, Y): 0.5005 0.3001
--- esempio 2
P(y|x) con insieme vuoto: [0.7, 0.4]  verita:  [0.7, 0.4]
aggiustando su W (errato): [0.7139, 0.382]
aggiustando su Z:          [0.7, 0.4]

L'ultima riga di ciascun gruppo mostra che la stima da un campione osservazionale, che usa soltanto $(X,W,Y)$ senza mai guardare $Z$, ritrova l'effetto causale: $0.5005$ e $0.3001$ contro $0.5$ e $0.3$. Nell'esempio 2 le liste sono nell'ordine $[x=0,\ x=1]$.