Esercizio - Effetto causale di pioggia e irrigatore sul prato bagnato (laboratorio 7)

Testo (laboratorio 7, Causal inference - Part A, esercizio 1). Si consideri la rete dello sprinkler, supposta una descrizione affidabile delle relazioni causali tra le quattro variabili: nuvoloso $C$, irrigatore $S$, pioggia $R$, erba bagnata $G$, con $C\to R$, $C\to S$, $R\to G$, $S\to G$ e le tabelle

1. Scrivere un programma Python che calcoli $P(g\mid do(r))$ visto a lezione.
2. Qual è l'effetto causale del tempo nuvoloso $c$ sul prato bagnato $g$, cioè l'espressione di $P(g\mid do(c))$?
3. Tra pioggia e irrigatore, quale ha l'effetto causale più forte, $P(g\mid do(r))$ o $P(g\mid do(s))$?
4. Calcolare l'effetto causale medio (ACE) dell'irrigatore su $g$.

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 →; il grafo è 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. $P(g\mid do(r))$

La pioggia non si può modificare fisicamente, ma si può calcolare il suo effetto causale «come se» si intervenisse, usando le probabilità osservazionali. L'unico genitore di $R$ è $C$, quindi la formula di aggiustamento è $$P(g\mid do(r))=\sum_cP(g\mid r,c)P(c)=P(g\mid r,c)P(c)+P(g\mid r,\neg c)P(\neg c).$$

$P(c)$ è nella rete. Resta $\mathbf P(G\mid r,c)=\langle P(g\mid r,c),P(\neg g\mid r,c)\rangle$, che si ricava dalla congiunta (Reti bayesiane - semantica, indipendenze e inferenza esattaUna rete bayesiana è un grafo diretto aciclico (DAG) con una variabile per nodo e una tabella delle probabilità condizionate (CPT) $P(X_i\mid\text{Parents}(X_i))$ per nodo; rappresenta la congiunta come prodotto $P(x_1,\dots,x_n)=\prod_i P(x_i\mid\text{parents}(X_i))$, con molti meno numeri della tabella completa (rete dell'allarme: 10 contro 31). Tre strutture elementari: catena $X\to Y\to Z$ e fork $Y\leftarrow X\to Z$ (estremi indipendenti dato il centro), collider $X\to Z\leftarrow Y$ (indipendenti, ma dipendenti se si osserva $Z$ o un suo discendente). Modelli causali (genitori = cause) danno reti più piccole. Inferenza esatta per enumerazione: $\mathbf P(X\mid\mathbf e)=\alpha\sum_{\mathbf y}\mathbf P(X,\mathbf e,\mathbf y)$ calcolata come somma di prodotti, in profondità e da sinistra a destra; per l'allarme $\mathbf P(B\mid j,m)=\langle0{,}284;\ 0{,}716\rangle$. Molti prodotti si ripetono: l'eliminazione di variabili li salva come fattori e li riusa.Reti bayesiane - semantica, indipendenze e inferenza esatta →): $$\mathbf P(G\mid r,c)=\frac{\mathbf P(G,r,c)}{P(r,c)}=\alpha\sum_s\mathbf P(G,r,c,s)=\alpha\sum_sP(c)P(r\mid c)P(s\mid c)\mathbf P(G\mid r,s)=\alpha'\sum_sP(s\mid c)\mathbf P(G\mid r,s).$$ Il prodotto $P(c)P(r\mid c)$ non dipende da $s$: si porta fuori dalla somma e si assorbe nella costante di normalizzazione ($\alpha'$ è una nuova costante, e il risultato sommato su $S$ è già normalizzato perché $P(s\mid c)+P(\neg s\mid c)=1$). Con i numeri: $$\mathbf P(G\mid r,c)=0.1\cdot\langle0.95,0.05\rangle+0.9\cdot\langle0.90,0.10\rangle=\langle0.095+0.81,\ 0.005+0.09\rangle=\langle0.905,0.095\rangle.$$ Per $\neg c$ cambiano i pesi $P(s\mid\neg c)=0.5$: $$\mathbf P(G\mid r,\neg c)=0.5\cdot\langle0.95,0.05\rangle+0.5\cdot\langle0.90,0.10\rangle=\langle0.925,0.075\rangle.$$ Sostituendo nella formula: $$P(g\mid do(r))=0.905\cdot0.5+0.925\cdot0.5=0.915.$$

import numpy as np

t, f = 0, 1                                   # indici di "vero" e "falso"
P_C = np.array([0.5, 0.5])                    # P(C)
P_S_C = np.array([[0.1, 0.5], [0.9, 0.5]])    # P(S | C):  riga = valore di S, colonna = valore di C
P_R_C = np.array([[0.8, 0.1], [0.2, 0.9]])    # P(R | C)
P_G_SR = np.array([[[0.95, 0.90], [0.90, 0.10]],     # P(G | S, R): indici [g][s][r]
                   [[0.05, 0.10], [0.10, 0.90]]])

# P(G | r, c) = sum_s P(s | c) P(G | r, s)
P_G_rc = P_S_C[t, t] * P_G_SR[:, t, t] + P_S_C[f, t] * P_G_SR[:, f, t]     # P(G | r, c)
P_G_rnc = P_S_C[t, f] * P_G_SR[:, t, t] + P_S_C[f, f] * P_G_SR[:, f, t]    # P(G | r, no c)
p_g_do_r = P_C[t] * P_G_rc[t] + P_C[f] * P_G_rnc[t]
print(P_G_rc, P_G_rnc, p_g_do_r)

Output:

[0.905 0.095] [0.925 0.075] 0.915

Controllo indipendente: simulando l'SCM dello sprinkler (quattro variabili esogene uniformi, funzioni a soglia) un milione di volte con $R$ forzata a vero, la frequenza di erba bagnata è $0.9153$.

2. $P(g\mid do(c))$

$C$ non ha genitori: l'intervento non toglie nessuna freccia e il grafo prima e dopo è lo stesso. Quindi $P(g\mid do(c))=P(g\mid c)$, che si ricava sommando sui valori di $R$ e $S$ (indipendenti dato $C$): $$P(g\mid c)=\sum_{r,s}P(r\mid c)P(s\mid c)P(g\mid r,s)=0.8\cdot0.1\cdot0.95+0.8\cdot0.9\cdot0.90+0.2\cdot0.1\cdot0.90+0.2\cdot0.9\cdot0.10=0.076+0.648+0.018+0.018=0.76.$$

P_g_c = [sum(P_R_C[r, c] * P_S_C[s, c] * P_G_SR[t, s, r] for r in (t, f) for s in (t, f)) for c in (t, f)]
print(P_g_c)       # [0.76, 0.5425]: P(g | c) e P(g | no c)

3. Pioggia o irrigatore?

$R$ e $S$ hanno lo stesso genitore $C$ e lo stesso figlio $G$: il calcolo di $P(g\mid do(s))$ è identico con $r$ e $s$ scambiati. $$\mathbf P(G\mid s,c)=P(r\mid c)\mathbf P(G\mid s,r)+P(\neg r\mid c)\mathbf P(G\mid s,\neg r)=0.8\langle0.95,0.05\rangle+0.2\langle0.90,0.10\rangle=\langle0.94,0.06\rangle,$$ $$\mathbf P(G\mid s,\neg c)=0.1\langle0.95,0.05\rangle+0.9\langle0.90,0.10\rangle=\langle0.905,0.095\rangle,\qquad P(g\mid do(s))=0.5\cdot0.94+0.5\cdot0.905=0.9225.$$

P_G_sc = P_R_C[t, t] * P_G_SR[:, t, t] + P_R_C[f, t] * P_G_SR[:, t, f]     # P(G | s, c)
P_G_snc = P_R_C[t, f] * P_G_SR[:, t, t] + P_R_C[f, f] * P_G_SR[:, t, f]    # P(G | s, no c)
p_g_do_s = P_C[t] * P_G_sc[t] + P_C[f] * P_G_snc[t]
print(P_G_sc, P_G_snc, p_g_do_s)

Output: [0.94 0.06] [0.905 0.095] 0.9225.

Confronto: $P(g\mid do(s))=0.9225>P(g\mid do(r))=0.915$, ma questo confronta solo la probabilità di erba bagnata quando la causa è forzata a vero, non quanto la causa cambia il risultato. Per «effetto più forte» conta di quanto cambia $P(g)$ passando da causa assente a causa presente, cioè l'ACE. Quello dell'irrigatore è $0.4625$ (punto 4); quello della pioggia, con lo stesso calcolo, $P(g\mid do(\neg r))=0.34$ e $ACE_R=0.915-0.34=0.575$. Quindi la pioggia ha l'effetto causale più forte, nonostante $P(g\mid do(r))<P(g\mid do(s))$. Il motivo: ciascuna causa cambia molto l'esito solo se l'altra è assente. L'effetto della pioggia è $0.05$ se l'irrigatore è acceso ($0.95-0.90$) e $0.80$ se è spento ($0.90-0.10$), con $P(s)=0.3$: $0.3\cdot0.05+0.7\cdot0.80=0.575$. L'effetto dell'irrigatore è $0.05$ se piove e $0.80$ se non piove, con $P(r)=0.45$: $0.45\cdot0.05+0.55\cdot0.80=0.4625$. L'irrigatore è acceso meno spesso della pioggia che cade, quindi nel caso medio «l'altra causa è assente» è più probabile per la pioggia. Nota: la soluzione del corso calcola solo $P(g\mid do(s))$ e non dà una conclusione esplicita sul confronto.

4. ACE dell'irrigatore

$$ACE=P(g\mid do(s))-P(g\mid do(\neg s)).$$ Per $do(\neg s)$: $$\mathbf P(G\mid\neg s,c)=0.8\langle0.90,0.10\rangle+0.2\langle0.10,0.90\rangle=\langle0.74,0.26\rangle,\qquad\mathbf P(G\mid\neg s,\neg c)=0.1\langle0.90,0.10\rangle+0.9\langle0.10,0.90\rangle=\langle0.18,0.82\rangle,$$ $$P(g\mid do(\neg s))=0.5\cdot0.74+0.5\cdot0.18=0.46,\qquad ACE=0.9225-0.46=0.4625.$$

P_G_nsc = P_R_C[t, t] * P_G_SR[:, f, t] + P_R_C[f, t] * P_G_SR[:, f, f]    # P(G | no s, c)
P_G_nsnc = P_R_C[t, f] * P_G_SR[:, f, t] + P_R_C[f, f] * P_G_SR[:, f, f]   # P(G | no s, no c)
p_g_do_ns = P_C[t] * P_G_nsc[t] + P_C[f] * P_G_nsnc[t]
print(P_G_nsc, P_G_nsnc, p_g_do_ns, p_g_do_s - p_g_do_ns)

Output: [0.74 0.26] [0.18 0.82] 0.46 0.4625 (a meno di errori di arrotondamento in virgola mobile). $ACE=0.4625>0$: accendere l'irrigatore aumenta la probabilità di erba bagnata di circa 46 punti percentuali. La simulazione dell'SCM con $S$ forzata dà $0.9228$ e $0.4597$.

Nota sul codice del corso. La soluzione ufficiale indicizza P_G_SR[:, t, f] dove qui compare P_G_SR[:, f, t]: i due valori coincidono ($P(g\mid s,\neg r)=P(g\mid\neg s,r)=0.90$), quindi i risultati sono identici.