Esercizio - Paradosso di Simpson con pyAgrum (laboratorio 7)

Testo (laboratorio 7, Causal inference - Part A, paradosso di Simpson con pyAgrum). In uno studio statistico su un farmaco si valuta l'efficacia del farmaco su una popolazione di uomini e donne. Le variabili sono Drug (assunzione del farmaco: Without o With), Patient (paziente guarito: Sick o Healed) e Gender (sesso: F o M). Il modello costruito dai dati osservati ha grafo $\text{Gender}\to\text{Drug}$, $\text{Gender}\to\text{Patient}$, $\text{Drug}\to\text{Patient}$ e tabelle (nota del corso: dati di ingresso diversi da quelli della lezione):

1. Calcolare con pyAgrum $P(\text{Healed}\mid\text{Drug})$ osservazionale, in generale e per ogni sesso.
2. Verificare che si ottiene il paradosso di Simpson.
3. «Dare il farmaco» non è un'osservazione ma un intervento: calcolare $P(\text{Healed}\mid do(\text{Drug}))$ con pyAgrum e concludere se prescrivere il farmaco.
4. (Esercizio finale del laboratorio) Risolvere lo stesso problema con i dati dei 700 pazienti della lezione (357 uomini e 343 donne; 350 hanno preso la medicina) e confrontare con il calcolo a mano delle slide.

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à → (il paradosso) 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 → (l'intervento).

Soluzione

Il modello in 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("Gender{F|M}->Drug{Without|With}->Patient{Sick|Healed}<-Gender")
m1.cpt("Gender")[:] = [0.5, 0.5]
m1.cpt("Drug")[:] = [[0.25, 0.75],        # Gender=F: Without 0.25, With 0.75
                     [0.75, 0.25]]        # Gender=M
m1.cpt("Patient")[{"Drug": "Without", "Gender": "F"}] = [0.2, 0.8]    # [Sick, Healed]
m1.cpt("Patient")[{"Drug": "Without", "Gender": "M"}] = [0.6, 0.4]
m1.cpt("Patient")[{"Drug": "With", "Gender": "F"}] = [0.3, 0.7]
m1.cpt("Patient")[{"Drug": "With", "Gender": "M"}] = [0.8, 0.2]

(I commenti del codice del corso scambiano «Male» e «Female» rispetto alle etichette $F$ e $M$ delle chiavi; i valori sono quelli della tabella sopra, coerenti con le didascalie dei risultati.)

1-2. Probabilità osservazionali e paradosso

pyAgrum calcola la distribuzione a posteriori con gum.getPosterior(m1, target="Patient", evs={...}); il corso la usa per tre casi: nessuna evidenza di sesso, $\text{Gender}=F$, $\text{Gender}=M$, ciascuno con $\text{Drug}=\text{Without}$ e $\text{With}$ (si prende l'indice [1], cioè Healed).

def getCuredObservedProba(m1, evs):
    evs0, evs1 = dict(evs), dict(evs)
    evs0["Drug"], evs1["Drug"] = "Without", "With"
    return gum.Potential().add(m1.variableFromName("Drug")).fillWith([
        gum.getPosterior(m1, target="Patient", evs=evs0)[1],
        gum.getPosterior(m1, target="Patient", evs=evs1)[1]])

Risultati osservazionali, che si controllano a mano con la congiunta $P(g)P(d\mid g)P(h\mid g,d)$:

• $P(\text{Drug}=\text{Without})=0.5\cdot0.25+0.5\cdot0.75=0.5$ e $P(\text{Healed},\text{Without})=0.5\cdot0.25\cdot0.8+0.5\cdot0.75\cdot0.4=0.1+0.15=0.25$, quindi $P(\text{Healed}\mid\text{Without})=0.25/0.5=0.5$;
• $P(\text{Healed},\text{With})=0.5\cdot0.75\cdot0.7+0.5\cdot0.25\cdot0.2=0.2625+0.025=0.2875$ e $P(\text{With})=0.5$, quindi $P(\text{Healed}\mid\text{With})=0.575$;
• per sesso: $F$: $0.8$ senza farmaco, $0.7$ con; $M$: $0.4$ senza, $0.2$ con.

Quindi $$P(C\mid\neg D)=0.5<P(C\mid D)=0.575,\qquad P(C\mid\neg D,F)=0.8>P(C\mid D,F)=0.7,\qquad P(C\mid\neg D,M)=0.4>P(C\mid D,M)=0.2.$$ È il paradosso di Simpson, in direzione opposta a quello della lezione: nei sottogruppi il farmaco peggiora la guarigione, sull'insieme sembra migliorarla. Perché: le donne, che guariscono molto di più, prendono il farmaco nel $75\%$ dei casi e gli uomini nel $25\%$, quindi tra chi prende il farmaco ci sono tre donne ogni uomo ($P(F\mid\text{With})=0.75$): il farmaco appare efficace perché è preso soprattutto da chi guarisce comunque.

3. Intervento

Il farmaco si decide: serve $P(\text{Healed}\mid do(\text{Drug}))$. Il modello causale proposto è lo stesso grafo, con $\text{Gender}$ confondente (causa di $\text{Drug}$ e di $\text{Patient}$). pyAgrum trova da solo l'insieme di aggiustamento («backdoor ['Gender'] found»):

d1 = csl.CausalModel(m1)
cslnb.showCausalImpact(d1, "Patient", doing="Drug", values={"Drug": "Without"})
cslnb.showCausalImpact(d1, "Patient", doing="Drug", values={"Drug": "With"})

La formula è la formula di aggiustamento su $\text{Gender}$: $$P(\text{Healed}\mid do(\text{Drug}))=\sum_{g}P(\text{Healed}\mid\text{Drug},g)\,P(g).$$

• senza farmaco: $0.8\cdot0.5+0.4\cdot0.5=0.6$;
• con il farmaco: $0.7\cdot0.5+0.2\cdot0.5=0.45$.

Dunque $P(\text{Healed}\mid do(\text{Without}))=0.6>P(\text{Healed}\mid do(\text{With}))=0.45$ (ACE del farmaco $=-0.15$): il farmaco non migliora la guarigione e conviene non prescriverlo. Con l'aggiustamento si pesano i due sottogruppi con le frequenze della popolazione ($0.5$ e $0.5$) anziché con quelle, molto diverse, di chi prende o no il farmaco.

Verifica con numpy

pyAgrum non è installato qui: i numeri sono ricalcolati con la funzione seguente (enumerazione della congiunta per le probabilità osservazionali e aggiustamento su $\text{Gender}$ per l'intervento).

import numpy as np

def analizza(p_gender, p_with_given_g, p_heal, nome):
    """p_gender[g]: P(Gender=g); p_with_given_g[g]: P(Drug=With | g); p_heal[g][d]: P(Healed | g, d), d=0 Without, d=1 With"""
    J = np.zeros((2, 2, 2))                                   # congiunta [g][d][h], h=1 guarito
    for g in range(2):
        for d in range(2):
            pd = p_with_given_g[g] if d == 1 else 1 - p_with_given_g[g]
            for h in range(2):
                ph = p_heal[g][d] if h == 1 else 1 - p_heal[g][d]
                J[g, d, h] = p_gender[g] * pd * ph
    print(nome)
    for d in range(2):                                         # osservazione: P(Healed | Drug)
        print("  P(Healed | Drug=%s) = %.4f" % (["Without", "With"][d], J[:, d, 1].sum() / J[:, d, :].sum()))
    for g in range(2):                                         # osservazione per sottogruppo
        print("  Gender=%s: " % "FM"[g] + ", ".join("P(Healed | %s) = %.4f" % (["Without", "With"][d], J[g, d, 1] / J[g, d, :].sum()) for d in range(2)))
    do = [sum(p_gender[g] * p_heal[g][d] for g in range(2)) for d in range(2)]   # aggiustamento su Gender
    print("  do(Without) = %.4f, do(With) = %.4f, ACE = %.4f" % (do[0], do[1], do[1] - do[0]))

analizza([0.5, 0.5], [0.75, 0.25], [[0.8, 0.7], [0.4, 0.2]], "Modello del laboratorio")     # indici: 0 = F, 1 = M

Output:

Modello del laboratorio
  P(Healed | Drug=Without) = 0.5000
  P(Healed | Drug=With) = 0.5750
  Gender=F: P(Healed | Without) = 0.8000, P(Healed | With) = 0.7000
  Gender=M: P(Healed | Without) = 0.4000, P(Healed | With) = 0.2000
  do(Without) = 0.6000, do(With) = 0.4500, ACE = -0.1500

4. I 700 pazienti della lezione

Dai conteggi della tabella (uomini: 87 con medicina, 81 guariti; 270 senza, 234 guariti; donne: 263 con, 192 guariti; 80 senza, 55 guariti) si costruiscono le tabelle del modello: $P(M)=357/700=0.51$, $P(F)=343/700=0.49$; $P(\text{With}\mid M)=87/357=0.2437$, $P(\text{With}\mid F)=263/343=0.7668$; $P(\text{Healed}\mid\text{With},M)=81/87$, $P(\text{Healed}\mid\text{Without},M)=234/270$, $P(\text{Healed}\mid\text{With},F)=192/263$, $P(\text{Healed}\mid\text{Without},F)=55/80$.

m2 = gum.fastBN("Gender{F|M}->Drug{Without|With}->Patient{Sick|Healed}<-Gender")
m2.cpt("Gender")[:] = [343 / 700, 357 / 700]
m2.cpt("Drug")[:] = [[80 / 343, 263 / 343], [270 / 357, 87 / 357]]
m2.cpt("Patient")[{"Drug": "Without", "Gender": "F"}] = [1 - 55 / 80, 55 / 80]
m2.cpt("Patient")[{"Drug": "Without", "Gender": "M"}] = [1 - 234 / 270, 234 / 270]
m2.cpt("Patient")[{"Drug": "With", "Gender": "F"}] = [1 - 192 / 263, 192 / 263]
m2.cpt("Patient")[{"Drug": "With", "Gender": "M"}] = [1 - 81 / 87, 81 / 87]

Si esegue sul secondo insieme di dati la stessa funzione analizza (qui con indici $0=F$, $1=M$):

n = {("F", 1): 263, ("F", 0): 80, ("M", 1): 87, ("M", 0): 270}     # pazienti per (sesso, medicina)
r = {("F", 1): 192, ("F", 0): 55, ("M", 1): 81, ("M", 0): 234}     # guariti
pg = [343 / 700, 357 / 700]
pw = [n["F", 1] / (n["F", 1] + n["F", 0]), n["M", 1] / (n["M", 1] + n["M", 0])]
ph = [[r["F", d] / n["F", d] for d in (0, 1)], [r["M", d] / n["M", d] for d in (0, 1)]]
analizza(pg, pw, ph, "Dati dei 700 pazienti")

Output:

Dati dei 700 pazienti
  P(Healed | Drug=Without) = 0.8257
  P(Healed | Drug=With) = 0.7800
  Gender=F: P(Healed | Without) = 0.6875, P(Healed | With) = 0.7300
  Gender=M: P(Healed | Without) = 0.8667, P(Healed | With) = 0.9310
  do(Without) = 0.7789, do(With) = 0.8325, ACE = 0.0537

Con questi dati il paradosso va in senso opposto al modello di prima: osservazionalmente $0.78<0.83$ ma in ogni sottogruppo la medicina migliora; l'intervento dà $P(\text{Healed}\mid do(\text{With}))=0.8325>0.7789$, ACE $=+0.0537$: la medicina funziona. La slide 37 del corso riporta $0.832$, $0.7818$ e $0.0502$ perché arrotonda le frazioni a due cifre prima di moltiplicare ($0.93$, $0.87$, ...): il segno e la conclusione non cambiano. Lo svolgimento a mano e l'analisi del caso in cui $\text{Gender}$ non sia un confondente sono in Esercizio - Paradosso di Simpson ed effetto causale della medicina (slide del corso).