Esercizio - Paradosso di Simpson ed effetto causale della medicina (slide del corso)

Testo (slide del corso di Intelligenza Artificiale, Ing. dell'Informazione UniPD, Causal Inference parte 1, slide 8, 36 e 37). Una nuova medicina è stata offerta a 700 pazienti (357 uomini e 343 donne): 350 hanno scelto di prenderla, 350 no.

La medicina funziona nei due sottogruppi, uomini e donne, ma non nella popolazione intera: come è possibile?

1. Spiegare l'apparente contraddizione con i numeri.
2. Si assuma il modello causale $Z\to X$, $Z\to Y$, $X\to Y$ con $X$ = medicina, $Y$ = guarigione, $Z$ = sesso. Calcolare con la formula di aggiustamento $P(Y=1\mid do(X=1))$, $P(Y=1\mid do(X=0))$ e l'effetto causale medio (ACE); concludere se la medicina migliora la guarigione.
3. Cosa cambierebbe se $Z$ non fosse il sesso ma una conseguenza della medicina (modello $X\to Z\to Y$)?
4. Controllare il risultato con una simulazione del modello.

Teoria: 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à → e 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 →. La versione con la libreria pyAgrum è in Esercizio - Paradosso di Simpson con pyAgrum (laboratorio 7).

Soluzione

1. Perché il paradosso non è una contraddizione

Le percentuali nei sottogruppi sono rapporti tra numeri diversi, e i totali sono medie pesate con pesi diversi nei due gruppi. Il totale di chi prende la medicina è $$\frac{273}{350}=\frac{87}{350}\cdot\frac{81}{87}+\frac{263}{350}\cdot\frac{192}{263}=0.2486\cdot0.931+0.7514\cdot0.730=0.78,$$ quello di chi non la prende $$\frac{289}{350}=\frac{270}{350}\cdot\frac{234}{270}+\frac{80}{350}\cdot\frac{55}{80}=0.7714\cdot0.867+0.2286\cdot0.6875=0.8257.$$ Tra chi prende la medicina gli uomini sono solo il $25\%$ (87 su 350); tra chi non la prende sono il $77\%$ (270 su 350). Gli uomini guariscono di più comunque (87-93% contro 69-73%), quindi il gruppo senza medicina, pieno di uomini, parte avvantaggiato e compensa (e supera) il vantaggio che la medicina dà nei sottogruppi. Il sesso è un confondente: influenza la scelta della medicina e la guarigione.

2. Aggiustamento sul sesso

$Z$ è l'unico genitore di $X$, quindi $P(y\mid do(x))=\sum_zP(y\mid x,z)P(z)$ con i pesi della popolazione intera: $P(Z=\text{uomo})=\frac{87+270}{700}=\frac{357}{700}=0.51$, $P(Z=\text{donna})=0.49$ (uguali per entrambi i valori di $x$).

$$P(Y{=}1\mid do(X{=}1))=\frac{81}{87}\cdot0.51+\frac{192}{263}\cdot0.49=0.9310\cdot0.51+0.7300\cdot0.49=0.4748+0.3577=0.8325,$$ $$P(Y{=}1\mid do(X{=}0))=\frac{234}{270}\cdot0.51+\frac{55}{80}\cdot0.49=0.8667\cdot0.51+0.6875\cdot0.49=0.4420+0.3369=0.7789.$$ $$ACE=0.8325-0.7789=0.0537>0.$$

La medicina migliora la guarigione di circa 5.4 punti percentuali: con $do(X=1)$ si dà la medicina a tutti, uomini e donne nelle proporzioni naturali, e si confrontano gli stessi pesi.

Le slide arrotondano a due cifre ($0.93$, $0.73$, $0.87$, $0.69$) e ottengono $0.832$, $0.7818$ e $ACE=0.0502$: stessa conclusione, differenza dovuta solo all'arrotondamento.

3. Se $Z$ fosse un mediatore

Si supponga che la variabile che divide in gruppi sia dopo la medicina e ne dipenda, per esempio la pressione del sangue dopo il trattamento: $X\to Z\to Y$ (senza confondente). Non ci sono cammini backdoor da $X$ a $Y$, quindi $P(y\mid do(x))=P(y\mid x)$ e vanno usati i dati aggregati: $$ACE=\frac{273}{350}-\frac{289}{350}=-\frac{16}{350}=-0.0457<0.$$ In questo modello la medicina peggiora la guarigione (la parte di effetto che passa per $Z$ è un effetto della medicina e non va «controllata»). I dati sono gli stessi: solo il grafo causale decide quale tabella usare. È questo il senso del paradosso: la matematica delle tabelle non risolve la domanda.

4. Verifica

from fractions import Fraction as F
import numpy as np

# dati: (sesso, medicina) -> (guariti, pazienti); medicina 1 = presa, 0 = non presa
dati = {("M", 1): (81, 87), ("M", 0): (234, 270), ("F", 1): (192, 263), ("F", 0): (55, 80)}
N = sum(n for _, n in dati.values())                                    # 700

# a) percentuali per sottogruppo e totali
for (z, x), (r, n) in dati.items():
    print(z, x, f"{r}/{n} = {r / n:.4f}")
for x in (1, 0):
    r = sum(dati[z, x][0] for z in "MF"); n = sum(dati[z, x][1] for z in "MF")
    print("totale, medicina =", x, f"{r}/{n} = {r / n:.4f}", " P(uomo | x) =", round(dati["M", x][1] / n, 4))

# b) aggiustamento su Z (sesso): P(y | do(x)) = sum_z P(y | x, z) P(z)
Pz = {z: F(sum(dati[z, x][1] for x in (0, 1)), N) for z in "MF"}        # P(uomo) = 357/700
do = {x: sum(F(*dati[z, x]) * Pz[z] for z in "MF") for x in (0, 1)}
print("P(uomo) =", Pz["M"], "  do(1) =", float(do[1]), " do(0) =", float(do[0]), " ACE =", float(do[1] - do[0]))

# c) se Z fosse un mediatore (X -> Z -> Y) non si aggiusta: P(y | do(x)) = P(y | x)
tot = {x: F(sum(dati[z, x][0] for z in "MF"), 350) for x in (0, 1)}
print("mediatore: ACE =", float(tot[1] - tot[0]))

# d) controllo con una simulazione dell'SCM Z -> X, Z -> Y, X -> Y
rng = np.random.default_rng(1)
M = 2_000_000
Z = rng.random(M) < float(Pz["M"])                                       # True = uomo
X = rng.random(M) < np.where(Z, 87 / 357, 263 / 343)                    # P(medicina | sesso)
Uy = rng.random(M)
p_guarigione = lambda X: np.where(Z, np.where(X, 81 / 87, 234 / 270), np.where(X, 192 / 263, 55 / 80))
print("osservato:", (Uy < p_guarigione(X))[X].mean().round(4), (Uy < p_guarigione(X))[~X].mean().round(4))
for x in (True, False):
    print("do(x =", x, ") =", (Uy < p_guarigione(np.full(M, x))).mean().round(4))

Output:

M 1 81/87 = 0.9310
M 0 234/270 = 0.8667
F 1 192/263 = 0.7300
F 0 55/80 = 0.6875
totale, medicina = 1 273/350 = 0.7800  P(uomo | x) = 0.2486
totale, medicina = 0 289/350 = 0.8257  P(uomo | x) = 0.7714
P(uomo) = 51/100   do(1) = 0.8325462173856037  do(0) = 0.778875  ACE = 0.05367121738560378
mediatore: ACE = -0.045714285714285714
osservato: 0.7802 0.8255
do(x = True ) = 0.8326
do(x = False ) = 0.7788

La simulazione del modello strutturale ($Z$ con $P=0.51$, $X$ che dipende dal sesso con $P(x\mid\text{uomo})=87/357$ e $P(x\mid\text{donna})=263/343$, $Y$ dalle tabelle di guarigione) riproduce le percentuali osservate ($0.78$ e $0.8257$) e, forzando $X$, i valori della formula di aggiustamento ($0.8326$ e $0.7788$).