Esercizio - Effetto causale della pioggia con pyAgrum (esempio sprinkler)

Testo (taccuino del corso Sprinkler example pyAgrum, laboratorio 7). Si consideri la rete dello sprinkler (nuvoloso, pioggia, irrigatore, erba bagnata) come descrizione affidabile delle relazioni causali. Non si può modificare fisicamente la pioggia $R$, ma si può usare l'osservazione del tempo per calcolare il suo effetto causale sul prato $G$ «come se» si potesse intervenire. Si calcoli $P(G=\text{true}\mid do(R=\text{true}))$ a mano con la formula di aggiustamento e con pyAgrum, e si confrontino i due risultati.

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 →. Lo stesso calcolo, in numpy, per tutte e quattro le richieste del laboratorio, è in Esercizio - Effetto causale di pioggia e irrigatore sul prato bagnato (laboratorio 7).

Soluzione

A mano

L'unico genitore di $R$ è $\text{Cloudy}$, quindi $$P(g\mid do(r))=P(g\mid r,c)P(c)+P(g\mid r,\neg c)P(\neg c),$$ con $P(c)=0.5$ dalla rete e, sommando sulla variabile nascosta $S$, $$\mathbf P(G\mid r,c)=\alpha'\big[P(s\mid c)\mathbf P(G\mid r,s)+P(\neg s\mid c)\mathbf P(G\mid r,\neg s)\big].$$

• Per $c$: $0.1\langle0.95,0.05\rangle+0.9\langle0.90,0.10\rangle=\langle0.905,0.095\rangle$.
• Per $\neg c$: $0.5\langle0.95,0.05\rangle+0.5\langle0.90,0.10\rangle=\langle0.925,0.075\rangle$.
• $P(g\mid do(r))=0.905\cdot0.5+0.925\cdot0.5=0.915$.

Con pyAgrum

import pyagrum as gum
import pyagrum.lib.notebook as gnb
import pyagrum.causal as csl
import pyagrum.causal.notebook as cslnb

m1 = gum.fastBN("Cloudy{T|F}->Rain{T|F}->GrassWet{T|F}<-Sprinkler{T|F}<-Cloudy{T|F}")
m1.cpt("Cloudy")[:] = [0.5, 0.5]
m1.cpt("Rain")[:] = [[0.8, 0.2],              # Cloudy = T
                     [0.1, 0.9]]              # Cloudy = F
m1.cpt("Sprinkler")[:] = [[0.1, 0.9],         # Cloudy = T
                          [0.5, 0.5]]         # Cloudy = F
m1.cpt("GrassWet")[{"Sprinkler": "T", "Rain": "T"}] = [0.95, 0.05]
m1.cpt("GrassWet")[{"Sprinkler": "T", "Rain": "F"}] = [0.90, 0.10]
m1.cpt("GrassWet")[{"Sprinkler": "F", "Rain": "T"}] = [0.90, 0.10]
m1.cpt("GrassWet")[{"Sprinkler": "F", "Rain": "F"}] = [0.10, 0.90]

d1 = csl.CausalModel(m1)                       # modello causale = rete + ipotesi che le frecce siano cause
cslnb.showCausalImpact(d1, "GrassWet", doing="Rain", values={"Rain": "T"})   # formula, valore, spiegazione
prob = cslnb.getCausalImpact(d1, "GrassWet", doing="Rain", values={"Rain": "T"})

Le tabelle sono scritte con la riga = valore del genitore e le due colonne = $[\text{T},\text{F}]$ del figlio. showCausalImpact mostra tre cose: la formula trovata (qui l'aggiustamento sul genitore $\text{Cloudy}$, spiegazione «backdoor found»), il valore e l'impatto; getCausalImpact restituisce la tabella come oggetto HTML. Il corso ne ricava il numero leggendo l'HTML con pandas.read_html e prendendo la prima cella:

import pandas as pd
df = pd.read_html(prob.data)[0]
print(df)               # colonne GrassWet = T, F
print(df.loc[0][0])     # valore per T

Output del taccuino del corso:

  GrassWet
         T      F
0    0.915  0.085
0.915

La riga $\langle0.915,0.085\rangle$ è la distribuzione $\mathbf P(G\mid do(r))$: la seconda cella vale $1-0.915$. Il risultato coincide con il calcolo a mano.

Verifica con numpy (pyAgrum non è installato qui)

Si calcola sia l'osservazione sia l'intervento con la chirurgia (togliendo il fattore $P(r\mid c)$ dalla congiunta e fissando $R$), che è un secondo modo di arrivare allo stesso numero senza usare la formula di aggiustamento:

import numpy as np

P_C = np.array([0.5, 0.5])                                # indici: 0 = T, 1 = F
P_R_C = np.array([[0.8, 0.1], [0.2, 0.9]])                # [r][c]
P_S_C = np.array([[0.1, 0.5], [0.9, 0.5]])                # [s][c]
P_G_SR = np.array([[[0.95, 0.90], [0.90, 0.10]],
                   [[0.05, 0.10], [0.10, 0.90]]])         # [g][s][r]

# osservazione: P(G=T | R=T) dalla congiunta completa
num = sum(P_C[c] * P_R_C[0, c] * P_S_C[s, c] * P_G_SR[0, s, 0] for c in range(2) for s in range(2))
den = sum(P_C[c] * P_R_C[0, c] * P_S_C[s, c] for c in range(2) for s in range(2))
print("P(g | r)     =", round(num / den, 4))

# intervento: congiunta dopo la chirurgia = P(c) P(s|c) P(g|s,r) con R fissata (manca il fattore P(r|c))
def do_rain(r):
    return sum(P_C[c] * P_S_C[s, c] * P_G_SR[0, s, r] for c in range(2) for s in range(2))
print("P(g | do(r)) =", round(do_rain(0), 4), "  P(g | do(no r)) =", round(do_rain(1), 4))

Output:

P(g | r)     = 0.9072
P(g | do(r)) = 0.915   P(g | do(no r)) = 0.34

A mano, il valore osservazionale è $P(g\mid r)=\dfrac{P(g,r)}{P(r)}$ con $P(g,r)=\sum_{c,s}P(c)P(r\mid c)P(s\mid c)P(g\mid r,s)=0.5\cdot0.8\cdot0.905+0.5\cdot0.1\cdot0.925=0.362+0.04625=0.40825$ e $P(r)=0.5\cdot0.8+0.5\cdot0.1=0.45$, quindi $0.40825/0.45=0.9072$. Qui compare il fattore $P(r\mid c)$ che l'intervento elimina: nell'osservazione i casi con pioggia sono per l'$89\%$ nuvolosi ($0.4/0.45$), e quando è nuvoloso l'irrigatore è quasi sempre spento; nell'intervento invece le nuvole hanno la frequenza $0.5$ della popolazione.

Il valore osservazionale $P(g\mid r)=0.9072$ è diverso da quello causale $0.915$: la pioggia è correlata con il tempo nuvoloso, che spegne di solito l'irrigatore; forzare la pioggia (a prescindere dalle nuvole) lascia l'irrigatore acceso in media nel $30\%$ dei casi ($P(s)=0.5\cdot0.1+0.5\cdot0.5$). Il secondo numero, $P(g\mid do(\neg r))=0.34$, dà l'ACE della pioggia $0.915-0.34=0.575$ (l'ACE dell'irrigatore è nell'esercizio del laboratorio 7). Il calcolo è confermato anche da una simulazione dell'SCM con un milione di campioni ($0.9153$).