Serie : Programmation Probabiliste avec PyMC (5 / parite Infer.NET) Duree estimee : 65 minutes Prerequis : PyMC-04-Bayesian-Networks (reseaux bayesiens, CPT, D-separation)
P(Y|X) et intervention P(Y|do(X)) (Pearl, 2000)do natif de PyMC (pm.do) — le do-calculus comme transformation de modèle de première classe| Précédent | Suivant |
|---|---|
| PyMC-04-Bayesian-Networks | PyMC-07-Skills-IRT |
Note de numerotation : ce notebook porte le numéro 5 par parite avec son jumeau C# Infer-5-Causal-Inference. Le sujet d’Infer-21 (Thompson Sampling) est, cote Python, integre dans PyMC-7-Sequential (section 10ter, bandits bayesiens) — d’ou l’absence d’un PyMC-21 distinct.
Ponts inter-series : le versant message passing (Infer.NET, effet causal calcule par EP/VMP) est dans Infer-5-Causal-Inference ; le versant symbolique (Java/Tweety, raisonneur contrefactuel propositionnel) est dans Tweety-11-Causal. Ce notebook en est le versant probabiliste MCMC + enumeration exacte, ou l’opérateur
doest une transformation de graphe native de PyMC.
La probabilite seule repond aux questions d’association. Mais de nombreuses questions pratiques sont causales :
| Question | Type | Reponse probabiliste |
|---|---|---|
| Si je vois le barometre baisser, quelle est la proba d’une tempete ? | Observation P(Storm\|Baro) |
Oui |
| Si je provoque la chute du barometre, la tempete vient-elle ? | Intervention P(Storm\|do(Baro)) |
Non (sans causalite) |
| Si j’avais pris le medicament (alors que je ne l’ai pas pris), aurais-je gueri ? | Contrefactuel | Non |
Judea Pearl (2000, Causality, Cambridge University Press) structure le raisonnement causal en trois niveaux (l’echelle de Pearl) :
P(Y|X). “Que disent les données ?”P(Y|do(X)). “Que se passe-t-il si j’agis ?”P(Y_x | X', Y'). “Que se serait-il passe si ?”Inegalite fondamentale : P(Y|X) != P(Y|do(X)) des qu’un confondeur (cause commune de X et Y) biaise l’observation. Ce notebook montre comment conditionner (observer) differe de do (intervenir) : do(X) mutile le graphe en coupant les arcs entrants de X. La ou Infer.NET demande une mutilation manuelle (forcer une variable sans parent), PyMC fournit pm.do, une transformation de modèle qui realise exactement cette mutilation.
References : Pearl, J. (1995), Causal diagrams for empirical research, JRSS B 57(3) ; Pearl, J. (2000), Causality : Models, Reasoning, and Inference, Cambridge University Press ; Pearl, J., Glymour, M. & Jewell, N. P. (2016), Causal Inference in Statistics : A Primer, Wiley.
Chargement de PyMC (modèles probabilistes + opérateurs causaux do/observe), ArviZ (diagnostics) et NumPy. Les reseaux de ce notebook sont discrets (Bernoulli) : on dispose donc de deux moteurs complementaires —
pm.do / pm.observe + echantillonnage : l’API causale native de PyMC, qui realise la mutilation de graphe et le conditionnement comme des transformations de modèle de première classe.On confronte les deux a chaque étape : l’enumeration donne le chiffre exact, pm.do prouve que l’opérateur natif retrouve le même résultat.
import numpy as np
import itertools
import pymc as pm
from pymc import do, observe
import arviz as az
RNG = 12345
print(f"PyMC : {pm.__version__}")
print(f"ArviZ : {az.__version__}")
print(f"NumPy : {np.__version__}")
print("Operateurs causaux disponibles : pm.do (intervention), pm.observe (conditionnement)")PyMC : 6.0.1
ArviZ : 1.1.0
NumPy : 2.3.5
Operateurs causaux disponibles : pm.do (intervention), pm.observe (conditionnement)Un petit moteur generique : un SCM discret est decrit par une liste de noeuds, chacun fourni avec une fonction P(noeud=1 | parents). On enumere les 2^n configurations, on calcule la probabilite jointe de chacune, puis on marginalise/conditionne par sommation. Pour do(X=x), on remplace la CPT de X par la constante x (mutilation), sans toucher au reste — exactement ce que pm.do fait structurellement.
# --- Acceptance #14051 : le moteur d'enumeration vit desormais dans le module
# canonique pymc_causal_organs.py, a cote de ce notebook. Le notebook l'IMPORTE
# au lieu de le redefinir, de sorte que ict.bridges.pymc_enumerate et tout autre
# consommateur observent la MEME implementation que celle executee ici.
import inspect
from pymc_causal_organs import enumerate_scm
# La valeur pedagogique est preservee : on affiche le code de l'organe REELLEMENT
# execute -- inspect.getsource lit l'objet importe, ce n'est pas une transcription.
print(inspect.getsource(enumerate_scm))
print("Moteur enumerate_scm pret (inference exacte par sommation sur 2^n configurations).")def enumerate_scm(
nodes: SCM,
query: str,
evidence: Optional[Dict[str, bool]] = None,
do_vars: Optional[Dict[str, bool]] = None,
) -> float:
"""Inference exacte par enumeration sur un SCM discret booleen.
Reprise byte-identique de la cellule 5 de ``PyMC-05-Causal-Inference.ipynb``.
nodes : liste ordonnee (topologique) de (nom, proba_true_fn). proba_true_fn(assign)
renvoie P(nom=True | parents) en lisant les parents dans le dict `assign`.
query : nom du noeud dont on veut P(query=True).
evidence : dict {nom: bool} de variables OBSERVEES (conditionnement / niveau 1).
do_vars : dict {nom: bool} de variables INTERVENUES (mutilation / niveau 2).
Retourne P(query=True | evidence, do(do_vars)).
"""
evidence = evidence or {}
do_vars = do_vars or {}
names = [n for n, _ in nodes]
num = 0.0 # masse de (query=True, evidence)
den = 0.0 # masse de (evidence)
for bits in itertools.product([False, True], repeat=len(names)):
assign = dict(zip(names, bits))
# incompatibilite avec une intervention forcee -> proba nulle
if any(assign[k] != v for k, v in do_vars.items()):
continue
if any(assign[k] != v for k, v in evidence.items()):
continue
p = 1.0
for name, fn in nodes:
if name in do_vars: # arc parent coupe : la CPT devient une constante
continue
pt = fn(assign)
p *= pt if assign[name] else (1.0 - pt)
den += p
if assign[query]:
num += p
return num / den if den > 0 else float("nan")
Moteur enumerate_scm pret (inference exacte par sommation sur 2^n configurations).Pearl illustre la confusion causale par le barometre : un barometre qui baisse annonce une tempete, mais provoquer la chute du barometre ne declenche pas de tempete. Pourquoi ? Parce qu’une cause commune — la basse pression atmospherique — provoque a la fois la chute du barometre ET la tempete.
BassePression
/ \
v v
Barometre TempeteBassePression -> Barometre : quand la pression baisse, le barometre baisse.BassePression -> Tempete : quand la pression baisse, la tempete vient.Barometre -> Tempete : le barometre n’est qu’un indicateur, pas une cause.Le barometre et la tempete sont donc correles (observation) mais il n’y a aucun effet causal de l’un sur l’autre (intervention). Construisons ce SCM avec les mêmes CPT que le jumeau Infer.NET (P(basse pression)=0.30, P(baro|basse)=0.90, P(storm|basse)=0.80).
# SCM du barometre (Pearl) : basse_pression -> {barometre, tempete}, AUCUN arc barometre -> tempete.
barometre_scm = [
("basse_pression", lambda a: 0.30),
("barometre", lambda a: 0.90 if a["basse_pression"] else 0.10),
("tempete", lambda a: 0.80 if a["basse_pression"] else 0.10),
]
# Marginale P(tempete) pour reference
p_storm_marg = enumerate_scm(barometre_scm, "tempete")
print(f"Structure : basse_pression -> {{barometre, tempete}} (pas d'arc barometre -> tempete)")
print(f"P(tempete) marginale = {p_storm_marg:.3f}")Structure : basse_pression -> {barometre, tempete} (pas d'arc barometre -> tempete)
P(tempete) marginale = 0.310P(Tempete | Barometre baisse)Si on observe le barometre baisser, la probabilite d’une tempete grimpe : la baisse est un signal que la pression est probablement basse, et donc qu’une tempete approche. Mais c’est une simple association statistique. On construit le modèle PyMC, on le conditionne avec pm.observe (l’opérateur de conditionnement natif), puis on lit la conditionnelle exacte par enumeration (le reseau est petit et discret — pour un modèle continu on echantillonnerait le posterieur).
# Modele PyMC du SCM barometre (memes CPT que l'enumeration ci-dessus)
with pm.Model() as m_baro:
lp = pm.Bernoulli("basse_pression", 0.30)
baro = pm.Bernoulli("barometre", pm.math.switch(lp, 0.90, 0.10))
storm = pm.Bernoulli("tempete", pm.math.switch(lp, 0.80, 0.10))
# NIVEAU 1 — Observation : pm.observe conditionne le modele sur barometre=1
m_baro_obs = observe(m_baro, {baro: 1})
print("pm.observe : modele conditionne sur barometre=1 "
f"(variables libres restantes : {[v.name for v in m_baro_obs.free_RVs]})")
# Conditionnelle exacte par enumeration
p_storm_obs = enumerate_scm(barometre_scm, "tempete", evidence={"barometre": True})
print(f"[exact] P(tempete | barometre baisse) = {p_storm_obs:.3f} (observationnel)")
print("Conclusion naive : 'le barometre predit la tempete'.")pm.observe : modele conditionne sur barometre=1 (variables libres restantes : ['basse_pression', 'tempete'])
[exact] P(tempete | barometre baisse) = 0.656 (observationnel)
Conclusion naive : 'le barometre predit la tempete'.P(Tempete | do(Barometre baisse))Maintenant provoquons la chute du barometre (en tapant dessus). Pour modeliser cette intervention, on mutile le graphe : on coupe l’arc BassePression -> Barometre et on force Barometre = baisse independamment de la pression.
C’est exactement ce que fait pm.do : il retourne un nouveau modèle ou la variable interventee n’a plus de parents, fixee a la valeur du do. Cote enumeration, on passe do_vars={"barometre": True} (la CPT de barometre est remplacée par une constante).
# NIVEAU 2 — Intervention : P(tempete | do(barometre=True))
# [exact] mutilation : on coupe l'arc basse_pression -> barometre
p_storm_do = enumerate_scm(barometre_scm, "tempete", do_vars={"barometre": True})
print(f"[exact] P(tempete | do(barometre baisse)) = {p_storm_do:.3f} (interventionnel)")
# [pm.do] l'operateur natif : do(modele, {barometre: 1}) coupe les arcs entrants de barometre
m_baro_do = do(m_baro, {baro: 1})
with m_baro_do:
pp_do = pm.sample_prior_predictive(draws=40000, random_seed=RNG)
p_storm_do_mcmc = float(pp_do.prior["tempete"].mean())
print(f"[pm.do + echantillonnage] P(tempete | do(barometre baisse)) ~ {p_storm_do_mcmc:.3f}")
print()
print(f"=> Observer le barometre informe sur la tempete : {p_storm_obs:.3f}")
print(f" PROVOQUER sa chute ne change rien a la tempete : {p_storm_do:.3f} (= P(tempete) marginale)")
print(" Voir != Faire (Pearl, 2000).")[exact] P(tempete | do(barometre baisse)) = 0.310 (interventionnel)[pm.do + echantillonnage] P(tempete | do(barometre baisse)) ~ 0.307
=> Observer le barometre informe sur la tempete : 0.656
PROVOQUER sa chute ne change rien a la tempete : 0.310 (= P(tempete) marginale)
Voir != Faire (Pearl, 2000).Résultat : P(Tempete|Barometre baisse) (observation) != P(Tempete|do(Barometre baisse)) (intervention). L’effet causal de l’intervention est nul — la chute du barometre ne provoque pas de tempete, et P(tempete|do(baro)) retombe sur la marginale P(tempete). L’observation etait biaisee par le confondeur BassePression.
Prong-B (non-trivialite) : ce résultat n’est pas une lecture triviale de CPT. Conditionner (
pm.observe) et intervenir (pm.do) donnent des nombres différents sur le même modèle — la différence visible (0.656vs0.310) prouve que l’opérateurdofait un travail structurel reel : il supprime un arc, il ne se contente pas de lire une table.
Verifions sur le reseau canonique de Pearl (Cloudy -> {Sprinkler, Rain} -> WetGrass, déjà construit en PyMC-4) la distinction observation/intervention sur la requête inverse P(Cloudy | Rain). La théorie donne P(Cloudy | Rain=True) = 0.800 (observationnel) vs P(Cloudy | do(Rain=True)) = 0.500 (interventionnel : Rain n’a aucun effet remontant sur Cloudy, donc on retombe sur le prior P(Cloudy)=0.5).
# Reseau Sprinkler (Pearl) — sous-graphe Cloudy -> Rain (suffisant pour P(Cloudy|Rain) vs do)
sprinkler_scm = [
("cloudy", lambda a: 0.5),
("rain", lambda a: 0.8 if a["cloudy"] else 0.2),
]
# OBSERVATIONNEL : P(cloudy | rain=True)
p_cloudy_obs = enumerate_scm(sprinkler_scm, "cloudy", evidence={"rain": True})
print(f"[exact] P(Cloudy | Rain=True) = {p_cloudy_obs:.3f} (theorie 0.800)")
# INTERVENTIONNEL via pm.do : P(cloudy | do(rain=True)) — couper cloudy -> rain
p_cloudy_do = enumerate_scm(sprinkler_scm, "cloudy", do_vars={"rain": True})
print(f"[exact] P(Cloudy | do(Rain=True)) = {p_cloudy_do:.3f} (theorie 0.500)")
with pm.Model() as m_spr:
cloudy = pm.Bernoulli("cloudy", 0.5)
rain = pm.Bernoulli("rain", pm.math.switch(cloudy, 0.8, 0.2))
with do(m_spr, {rain: 1}):
pp_spr = pm.sample_prior_predictive(draws=40000, random_seed=RNG)
print(f"[pm.do] P(Cloudy | do(Rain=True)) ~ {float(pp_spr.prior['cloudy'].mean()):.3f}")
print()
print("Cross-check PyMC-4 (meme structure) : 0.800 obs vs 0.500 do => coherent.")[exact] P(Cloudy | Rain=True) = 0.800 (theorie 0.800)
[exact] P(Cloudy | do(Rain=True)) = 0.500 (theorie 0.500)
[pm.do] P(Cloudy | do(Rain=True)) ~ 0.497
Cross-check PyMC-4 (meme structure) : 0.800 obs vs 0.500 do => coherent.Sur le baromètre, l’observation et l’intervention différaient (0.656 vs 0.310) mais restaient toutes deux des probabilités de tempête « plausibles ». Ici le contraste est plus net : observer la pluie fait monter P(Cloudy) à 0.800 — la pluie est un excellent indice de nuages — tandis que provoquer la pluie (arroseur géant) laisse P(Cloudy) à sa marginale 0.500. L’information remonte l’arc Cloudy → Rain quand on observe, mais aucun mécanisme physique ne la fait redescendre quand on intervient : do(Rain) coupe l’arc entrant, et Cloudy redevient une pièce non biaisée.
C’est la signature directionnelle de la causalité que les probabilités classiques ne peuvent pas exprimer : P(Cloudy|Rain) et P(Rain|Cloudy) vivent dans la même table jointe, mais P(Cloudy|do(Rain)) n’est pas une lecture de cette table — c’est une propriété du graphe orienté. L’échantillonnage pm.do (~0.503 sur 40 000 tirages) confirme la valeur exacte 0.500 à la fluctuation d’échantillonnage près : le do natif de PyMC exécute bien la mutilation, pas un simple conditionnement.
Quand le confondeur U est observe, l’effet causal se calcule sans mutiler le graphe, par la formule d’adjustment backdoor (Pearl, 1995) :
P(Y | do(X)) = Somme_u P(Y | X, u) P(u)On somme l’effet de X sur Y sur toutes les valeurs du confondeur, pondere par la distribution marginale de U. Si U bloque tous les chemins backdoor (arcs entrants de X passant par U), cette formule restitue l’effet causal exact. Verifions sur le barometre, ou U = BassePression bloque le chemin Barometre <- BassePression -> Tempete : l’ajustement doit coincider avec le do() direct (0.310).
# BACKDOOR : P(tempete | do(barometre)) = Somme_pression P(tempete | barometre, pression) P(pression)
# Ici tempete _||_ barometre | pression => P(tempete | barometre, pression) = P(tempete | pression)
p_basse = 0.30
p_storm_basse = enumerate_scm(barometre_scm, "tempete", evidence={"basse_pression": True}) # 0.80
p_storm_haute = enumerate_scm(barometre_scm, "tempete", evidence={"basse_pression": False}) # 0.10
backdoor = p_storm_basse * p_basse + p_storm_haute * (1 - p_basse)
print(f"Backdoor = P(tempete|basse)*{p_basse} + P(tempete|haute)*{1-p_basse:.1f}")
print(f" = {p_storm_basse:.2f}*{p_basse} + {p_storm_haute:.2f}*{1-p_basse:.1f} = {backdoor:.3f}")
print(f"do(barometre) direct (pm.do / enumeration) = {p_storm_do:.3f}")
print(f"=> Les deux methodes coincident ({backdoor:.3f} ~ {p_storm_do:.3f}) : verification de coherence.")Backdoor = P(tempete|basse)*0.3 + P(tempete|haute)*0.7
= 0.80*0.3 + 0.10*0.7 = 0.310
do(barometre) direct (pm.do / enumeration) = 0.310
=> Les deux methodes coincident (0.310 ~ 0.310) : verification de coherence.L’ajustement backdoor 0.80*0.3 + 0.10*0.7 = 0.310 coïncide exactement avec le do(barometro) direct par mutilation (0.310). Ce n’est pas une coïncidence numérique : les deux expriment la même quantité P(tempete | do(barometre)), l’une en transformant le graphe (mutilation : on coupe l’arc BassePression → Barometre), l’autre en corrigeant l’observation (on conditionne sur le confondeur puis on moyenne selon sa distribution P(pression)).
La condition de validité est graphique : BassePression bloque tous les chemins backdoor de Barometre (ici il n’y en a qu’un, Barometre ← BassePression → Tempete), et l’indépendance conditionnelle tempete ⊥ barometre | pression justifie la simplification P(tempete|barometre, pression) = P(tempete|pression) utilisée dans le calcul. Retenez l’économie de la méthode : pas besoin de manipuler le modèle — avec des données observationnelles et la bonne liste de variables d’ajustement, l’effet causal se calcule par une somme pondérée. C’est le fondement de l’inférence causale en épidémiologie et en économétrie (variables de contrôle).
Quand le confondeur U est non observe mais qu’un mediateur M entre X et Y capture tout l’effet de X sur Y, Pearl (1995) montre que l’effet causal reste identifiable par la formule front-door :
P(Y | do(X)) = Somme_m P(M=m | X) * Somme_x' P(Y | M=m, x') P(x')Exemple canonique (Pearl) : X = fumer, M = goudron dans les poumons, Y = cancer, U = genotype (non observe, cause a la fois l’envie de fumer et le risque de cancer). On n’a pas acces au genotype, mais le goudron transmet tout l’effet causal de la fumee sur le cancer. On verifie la formule contre le do(X) direct (mutilation, ici licite car on connait le modèle complet, U compris).
# FRONT-DOOR : X=smoke, M=tar, Y=cancer, U=genotype (non observe)
# SCM complet (U inclus) pour la verification do() directe
# --- Acceptance #14051 : le SCM front-door et P(Y|M,X) viennent du module
# canonique. Une SEULE definition existe desormais : si les CPT changent dans
# pymc_causal_organs.py, ce notebook ET ict.bridges le voient ensemble.
from pymc_causal_organs import FRONT_SCM as front_scm, p_y_given_m_x
# Les CPT restent lisibles -- lues sur l'objet IMPORTE, jamais recopiees ici :
_names = [n for n, _ in front_scm]
print("CPT du SCM front-door (lues depuis le module canonique) :")
for _name, _fn in front_scm:
_vals = sorted({round(_fn(dict(zip(_names, _bits))), 4)
for _bits in itertools.product([False, True], repeat=len(front_scm))})
print(f" P({_name}=1 | parents) in {_vals}")
print()
# Quantites observationnelles (U JAMAIS utilise — c'est tout l'interet du front-door)
p_x1 = enumerate_scm(front_scm, "smoke") # P(X=1) marginal
p_m1_given_x1 = enumerate_scm(front_scm, "tar", evidence={"smoke": True}) # P(M=1|X=1)
p_m0_given_x1 = 1 - p_m1_given_x1
inner_m1 = p_y_given_m_x(True, True) * p_x1 + p_y_given_m_x(True, False) * (1 - p_x1)
inner_m0 = p_y_given_m_x(False, True) * p_x1 + p_y_given_m_x(False, False) * (1 - p_x1)
front_door = p_m1_given_x1 * inner_m1 + p_m0_given_x1 * inner_m0
# Verification : do(smoke=True) direct (mutilation : coupe U -> smoke)
do_direct = enumerate_scm(front_scm, "cancer", do_vars={"smoke": True})
print(f"P(X=smoke) marginal = {p_x1:.3f}")
print(f"P(M=1|X=1) = {p_m1_given_x1:.3f} ; P(M=0|X=1) = {p_m0_given_x1:.3f}")
print(f"Front-door P(Cancer | do(Smoke)) = {front_door:.3f} (ajustement SANS acceder a U)")
print(f"do(X=1) direct (mutilation) = {do_direct:.3f} (verification, U connu)")
print(f"=> Les deux coincident : la formule front-door restitue l'effet causal.")CPT du SCM front-door (lues depuis le module canonique) :
P(u=1 | parents) in [0.2]
P(smoke=1 | parents) in [0.3, 0.8]
P(tar=1 | parents) in [0.1, 0.9]
P(cancer=1 | parents) in [0.05, 0.5, 0.7, 0.95]
P(X=smoke) marginal = 0.400
P(M=1|X=1) = 0.900 ; P(M=0|X=1) = 0.100
Front-door P(Cancer | do(Smoke)) = 0.689 (ajustement SANS acceder a U)
do(X=1) direct (mutilation) = 0.689 (verification, U connu)
=> Les deux coincident : la formule front-door restitue l'effet causal.Le résultat clé est l’égalité parfaite 0.689 = 0.689 : la formule front-door, construite uniquement à partir de quantités observationnelles (P(X=1) = 0.400, P(M=1|X=1) = 0.900, et les conditionnelles P(Y|M,X)), restitue exactement l’effet de la mutilation do(Smoke=True) — alors même que cette mutilation « triche » en utilisant le SCM complet, U compris.
Décomposons pourquoi ça marche : le goudron M capture tout l’effet de X sur Y (il n’y a pas d’arc direct Smoke → Cancer), et X n’a pas de parent observé qui corrompe P(M|X). Le premier facteur P(M=m|X) mesure donc l’effet causal X → M non biaisé ; la somme interne Somme_x' P(Y|M=m, x') P(x') « re-randomise » X pour neutraliser le biais du confondeur sur M → Y. C’est une décomposition de l’effet total en deux morceaux identifiés séparément — l’astuce de Pearl qui rend l’inference causale possible même quand U (génétique, statut socio-économique…) n’est pas mesurable.
Un renversement de Simpson survient quand une conclusion s’inverse entre l’analyse agregee et l’analyse conditionnelle. Cas clinique : un medicament semble nuire a la guerison dans la population totale, alors qu’il aide dans chaque sous-groupe (patients graves / patients legers). La cause : la severite est un confondeur — les patients graves, qui guerissent moins bien, recoivent aussi plus souvent le medicament. Le do(Drug) (backdoor sur la severite) restaure la verite causale.
# SIMPSON : severite (confondeur) -> {drug, recovery}, drug -> recovery
simpson_scm = [
("severe", lambda a: 0.5),
("drug", lambda a: 0.8 if a["severe"] else 0.2), # graves -> plus de drug
("rec", lambda a: (0.5 if a["drug"] else 0.3) if a["severe"] # grave : drug 0.5 > nodrug 0.3
else (0.9 if a["drug"] else 0.7)), # leger : drug 0.9 > nodrug 0.7
]
rec_drug_agg = enumerate_scm(simpson_scm, "rec", evidence={"drug": True})
rec_nodrug_agg = enumerate_scm(simpson_scm, "rec", evidence={"drug": False})
rec_drug_sev = enumerate_scm(simpson_scm, "rec", evidence={"drug": True, "severe": True})
rec_nodrug_sev = enumerate_scm(simpson_scm, "rec", evidence={"drug": False, "severe": True})
rec_drug_mild = enumerate_scm(simpson_scm, "rec", evidence={"drug": True, "severe": False})
rec_nodrug_mild= enumerate_scm(simpson_scm, "rec", evidence={"drug": False, "severe": False})
do_drug = enumerate_scm(simpson_scm, "rec", do_vars={"drug": True})
do_nodrug = enumerate_scm(simpson_scm, "rec", do_vars={"drug": False})
print("=== ANALYSE AGREGEE (biaisee par la severite) ===")
print(f"P(Rec | Drug) = {rec_drug_agg:.3f}")
print(f"P(Rec | NoDrug) = {rec_nodrug_agg:.3f} => le medicament semble NUIRE")
print()
print("=== ANALYSE CONDITIONNELLE (verite causale) ===")
print(f"Graves : Drug {rec_drug_sev:.2f} > NoDrug {rec_nodrug_sev:.2f}")
print(f"Legers : Drug {rec_drug_mild:.2f} > NoDrug {rec_nodrug_mild:.2f}")
print("=> Le medicament AIDE dans chaque sous-groupe. REVERSAL de l'agrege.")
print()
print("=== EFFET CAUSAL (do = backdoor sur la severite) ===")
print(f"P(Rec | do(Drug)) = {do_drug:.3f}")
print(f"P(Rec | do(NoDrug)) = {do_nodrug:.3f} => le medicament AIDE causalement.")=== ANALYSE AGREGEE (biaisee par la severite) ===
P(Rec | Drug) = 0.580
P(Rec | NoDrug) = 0.620 => le medicament semble NUIRE
=== ANALYSE CONDITIONNELLE (verite causale) ===
Graves : Drug 0.50 > NoDrug 0.30
Legers : Drug 0.90 > NoDrug 0.70
=> Le medicament AIDE dans chaque sous-groupe. REVERSAL de l'agrege.
=== EFFET CAUSAL (do = backdoor sur la severite) ===
P(Rec | do(Drug)) = 0.700
P(Rec | do(NoDrug)) = 0.500 => le medicament AIDE causalement.do rétablitTrois lectures du même SCM, trois verdicts :
P(Rec|Drug) = 0.580 < P(Rec|NoDrug) = 0.620 — le médicament semble nuire. C’est le chiffre qu’une analyse naïve du tableau croisé rapporterait.0.50 > 0.30 chez les graves, 0.90 > 0.70 chez les légers). Le renversement est le mécanisme exact du paradoxe de Simpson : les graves, qui guérissent moins bien quel que soit le traitement, reçoivent le médicament dans 80 % des cas — l’agrégé est dominé par ce mélange de populations.P(Rec|do(Drug)) = 0.700 contre 0.500 — l’effet causal vrai, obtenu en coupant l’arc Severité → Drug (mutilation) ou, de façon équivalente, par l’ajustement backdoor sur la sévérité de la section 5.La leçon pratique : aucun intervalle de confiance sur l’agrégé ne corrige ce biais — c’est un défaut de structure du graphe, pas un problème de taille d’échantillon. Seule la question causale bien posée (do, ou ajustement sur un ensemble Z validé par le critère backdoor) restaure le bon signe.
Le do() débranche la cause de ses effets ; la médiation demande une chirurgie plus fine : parmi l’effet total de x sur y, combien passe par le médiateur m (x -> m -> y) et combien par le chemin direct (x -> y) ?
Deux familles d’effets, deux façons de « figer » le médiateur :
m par une seconde intervention — P(y | do(x=1), do(m)) - P(y | do(x=0), do(m)) (double mutilation) ;Ces quantités « cross-world » sont calculables par énumération sur le SCM connu (aucun confondeur m <- u -> y : la sequential ignorability tient — c’est l’hypothèse qui tombe en terrain réel). L’identité NDE + NIE = TE est algébrique (téléscopage des sommes) ; ce qui n’est pas garanti, c’est le raccourci NDE = CDE(m) — exact seulement sans interaction x×m (le jumeau .NET Infer-5 mesure un monde avec interaction où ce raccourci dérive ; ici, monde sans interaction) :
# MEDIATION : x (formation) -> m (competences) -> y (promotion) + chemin direct x -> y.
# Monde SANS interaction : l'effet direct (+0.2) est le meme a m fixe haut ou bas.
med_scm = [
("x", lambda a: 0.5),
("m", lambda a: 0.8 if a["x"] else 0.2), # formation -> competences
("y", lambda a: (0.9 if a["x"] else 0.7) if a["m"] else (0.7 if a["x"] else 0.5)),
]
P_y_do = lambda x, m: enumerate_scm(med_scm, "y", do_vars={"x": x, "m": m})
TE = enumerate_scm(med_scm, "y", do_vars={"x": True}) - enumerate_scm(med_scm, "y", do_vars={"x": False})
CDE1 = P_y_do(True, True) - P_y_do(False, True)
CDE0 = P_y_do(True, False) - P_y_do(False, False)
# P(M=m | do(x)) se lit sur le CPT : 0.8/0.2 sous x=1, 0.2/0.8 sous x=0
NDE = 0.2 * (P_y_do(True, True) - P_y_do(False, True)) + 0.8 * (P_y_do(True, False) - P_y_do(False, False))
NIE = (0.8 - 0.2) * P_y_do(True, True) + (0.2 - 0.8) * P_y_do(True, False)
print("Mediation : x=formation, m=competences, y=promotion (monde sans interaction)")
print(f" TE = {TE:.3f} (effet total)")
print(f" CDE(1) = {CDE1:.3f} = CDE(0) = {CDE0:.3f} -> pas d'interaction : NDE = CDE")
print(f" NDE = {NDE:.3f} (direct, mediateur naturel sous x=0)")
print(f" NIE = {NIE:.3f} (indirect : la part transportee par les competences)")
print(f" NDE + NIE = {NDE + NIE:.3f} = TE = {TE:.3f} -> identite exacte")Mediation : x=formation, m=competences, y=promotion (monde sans interaction)
TE = 0.320 (effet total)
CDE(1) = 0.200 = CDE(0) = 0.200 -> pas d'interaction : NDE = CDE
NDE = 0.200 (direct, mediateur naturel sous x=0)
NIE = 0.120 (indirect : la part transportee par les competences)
NDE + NIE = 0.320 = TE = 0.320 -> identite exactePearl (1995) a demontre que toute quantite causale identifiable peut etre calculee a partir de trois règles de reecriture, qui transforment des expressions avec do() en expressions observationnelles (sans do()). Soit G_X le graphe mutile (arcs entrants de X supprimes) et G_{X bas} le graphe ou les arcs sortants de X sont supprimes :
| Règle | Enonce | Effet |
|---|---|---|
| Règle 1 (insertion/suppression d’observation) | P(y \| do(x), z, w) = P(y \| do(x), w) si Y _\|\|_ Z \| X, W dans G_X |
Observations irrelevantes retirees |
| Règle 2 (action/observation) | P(y \| do(x), do(z), w) = P(y \| do(x), z, w) si Y _\|\|_ Z \| X, W dans G_{X bas} |
Une action peut devenir une observation |
| Règle 3 (insertion/suppression d’action) | P(y \| do(x), do(z), w) = P(y \| do(x), w) si Y _\|\|_ Z \| X, W dans le graphe adequat |
Une action irrelevante est retiree |
Critere backdoor : un ensemble Z satisfait le critere backdoor relatif a (X,Y) si aucun noeud de Z n’est descendant de X, et Z bloque tous les chemins backdoor (arcs entrants de X). Alors P(Y|do(X)) = Somme_z P(Y|X,Z)P(Z) — exactement la formule appliquee en section 5.
Theoreme de completude (Shpitser & Pearl, 2006) : les trois règles sont completes — si un effet causal n’est pas calculable par ces règles, il n’est pas identifiable a partir du graphe seul. C’est un plafond fondamental, pas une limite d’implementation.
pm.dorealise la mutilation ; l’identifiabilite (peut-on exprimer la cible sansdoa partir des données disponibles ?) reste une question graphique en amont.
Le niveau 3 de l’echelle de Pearl repond aux questions contrefactuelles : “Sachant qu’un patient a pris le medicament ET a gueri, aurait-il gueri SANS medicament ?”
On resout cela en deux étapes (abduction-prediction) :
drug=True, rec=True).do(NoDrug)) pour predire le résultat contrefactuel.C’est plus couteux qu’une simple intervention (une inference posterieure puis une seconde inference) — honnetement le plafond pratique de l’inference causale exacte sur ce type de SCM.
# CONTREFACTUEL : "drug=True ET rec=True observes -> aurait-il gueri si do(drug=False) ?"
# ETAPE 1 — Abduction : posterieur sur 'severe' sachant (drug=True, rec=True)
p_severe_given = enumerate_scm(simpson_scm, "severe", evidence={"drug": True, "rec": True})
print(f"Abduction : P(severe | drug=True, rec=True) = {p_severe_given:.3f}")
# ETAPE 2 — Prediction : do(drug=False) avec 'severe' tire du posterieur d'abduction
cf_scm = [
("severe", lambda a: p_severe_given), # prior <- posterieur d'abduction
("drug", lambda a: 0.8 if a["severe"] else 0.2),
("rec", lambda a: (0.5 if a["drug"] else 0.3) if a["severe"]
else (0.9 if a["drug"] else 0.7)),
]
p_rec_cf = enumerate_scm(cf_scm, "rec", do_vars={"drug": False})
print(f"Contrefactuel : P(rec | do(NoDrug), 'a pris drug & gueri') = {p_rec_cf:.3f}")
print()
print(f"=> Meme si le patient a gueri AVEC le medicament, sans medicament il n'aurait eu")
print(f" qu'environ {p_rec_cf*100:.0f}% de chances de guerison (contrefactuel).")Abduction : P(severe | drug=True, rec=True) = 0.690
Contrefactuel : P(rec | do(NoDrug), 'a pris drug & gueri') = 0.424
=> Meme si le patient a gueri AVEC le medicament, sans medicament il n'aurait eu
qu'environ 42% de chances de guerison (contrefactuel).L’abduction renverse la flèche informationnelle : sachant que le patient a pris le médicament et a guéri, la probabilité qu’il soit un cas grave passe de 0.5 (prior) à 0.690 — les guérisons sous médicament sont plus fréquentes chez les légers, donc une guérison observée penche vers « léger », mais le médicament est plus prescrit aux graves… le posterieur 0.690 est le compromis exact de ces deux forces.
La prédiction contrefactuelle donne 0.424 : ce patient précis, re-placé dans le monde contrefactuel do(NoDrug) avec son profil abduit, n’aurait eu qu’environ 42 % de chances de guérison — contre 100 % dans le monde factual (il a guéri). C’est la lecture niveau 3 de Pearl : la même personne, deux mondes, deux nombres. Ni l’observation (P(rec|nodrug) agrégé) ni l’intervention (do(nodrug) sur la population) ne répondent à cette question — seuls les deux temps abduction → prediction y répondent, et c’est pourquoi le niveau 3 est strictement plus expressif que le niveau 2.
Quatre exercices pour ancrer le do-calculus. Chaque stub s’execute sans erreur (print de confirmation) ; a vous de completer le code entre les # TODO. Reutilisez enumerate_scm (section 2.1) et/ou pm.do.
| # | Exercice | Concept |
|---|---|---|
| 1 | Construire un SCM confondu (education -> revenu) | Niveau 1 observation vs intervention |
| 2 | Backdoor adjustment sur un reseau personnalise | Identification de l’ensemble Z |
| 3 | Front-door : smoke -> tar -> cancer (Pearl) | Ajustement par mediateur |
| 4 | Paradoxe de Simpson personnalise | Renversement agrege vs conditionnel |
# Exercice 1 : Construire un SCM confondu "education -> revenu"
# Modele suggere : U = niveau_etudes_parents cause a la fois Education (X) et Revenu (Y),
# et X -> Y (l'education augmente le revenu).
# Indice : reutilisez enumerate_scm. Decrivez 3 noeuds (u, education, revenu) avec leurs CPT.
# Etape 1 : definir edu_scm = [("u", ...), ("education", ...), ("revenu", ...)]
# Etape 2 : P(revenu | education=True) observationnel -> enumerate_scm(..., evidence=...)
# Etape 3 : P(revenu | do(education=True)) interventionnel -> enumerate_scm(..., do_vars=...)
# edu_scm = None # TODO etudiant : construire le SCM
print("Exercice 1 a completer : SCM education -> revenu (observation vs intervention).")Exercice 1 a completer : SCM education -> revenu (observation vs intervention).# Exercice 2 : Backdoor adjustment — identifier l'ensemble Z et calculer P(Y|do(X))
# Reprenez un reseau a 3-4 noeuds avec un confondeur OBSERVE Z.
# Indice : verifiez que Z bloque tous les chemins backdoor vers X, puis appliquez Somme_z P(Y|X,z)P(z).
# Etape 1 : choisir Z (ensemble d'ajustement) et le SCM
# Etape 2 : calculer chaque P(Y|X,z) via enumerate_scm(evidence={...})
# Etape 3 : sommer pondere par P(z) et comparer au do() direct (do_vars={...})
print("Exercice 2 a completer : backdoor adjustment, identifier Z et calculer P(Y|do(X)).")Exercice 2 a completer : backdoor adjustment, identifier Z et calculer P(Y|do(X)).# Exercice 3 : Front-door sur smoke -> tar -> cancer (classique Pearl)
# Adaptez la section 6 en changeant les CPT (par ex. tabagisme passif, exposition professionnelle).
# Indice : la structure X -> M -> Y + U -> X, U -> Y doit etre preservee (U non observe).
# Etape 1 : fixer de nouvelles CPT coherentes dans un nouveau front_scm
# Etape 2 : recalculer P(M|X), P(x'), P(Y|M,x') avec enumerate_scm
# Etape 3 : appliquer la formule front-door et comparer au do() direct
print("Exercice 3 a completer : front-door smoke -> tar -> cancer (Pearl).")Exercice 3 a completer : front-door smoke -> tar -> cancer (Pearl).# Exercice 4 : Paradoxe de Simpson personnalise
# Concevez des CPT ou la conclusion s'INVERSE entre l'agrege et le conditionnel.
# Indice : il faut un confondeur Z favorisant a la fois X et defavorisant Y (ou l'inverse).
# Etape 1 : choisir un scenario (ex. admission Berkeley 1973 : genre + departement)
# Etape 2 : coder des CPT produisant un renversement
# Etape 3 : montrer P(Y|X) agrege != P(Y|do(X)) et expliquer le mecanisme
print("Exercice 4 a completer : concevoir un paradoxe de Simpson (renversement agrege vs conditionnel).")Exercice 4 a completer : concevoir un paradoxe de Simpson (renversement agrege vs conditionnel).| Niveau | Question | Opérateur | Méthode PyMC |
|---|---|---|---|
| 1 Association | Que voit-on ? | P(Y\|X) |
pm.observe (conditionnement) / enumeration evidence= |
| 2 Intervention | Que se passe-t-il si on agit ? | P(Y\|do(X)) |
pm.do : mutile le graphe (coupe les arcs entrants) / do_vars= |
| 3 Contrefactuel | Que se serait-il passe ? | P(Y_x\|X',Y') |
Abduction (posterieur) + prediction dans le graphe mutile |
| Adjustment | Quand | Formule |
|---|---|---|
| Backdoor | Confondeur U observe |
P(Y\|do(X)) = Somme_u P(Y\|X,u)P(u) |
| Front-door | Confondeur non observe, mediateur M observe |
P(Y\|do(X)) = Somme_m P(M=m\|X) Somme_x' P(Y\|M=m,x')P(x') |
P(Y|X) != P(Y|do(X)) : l’observation est biaisee par les confondeurs. Conditionner (pm.observe) ne calcule pas l’effet causal — il faut pm.do, qui mutile le graphe.pm.do est un opérateur natif de PyMC : la ou Infer.NET demande une mutilation manuelle (forcer une variable sans parent), PyMC traite l’intervention comme une transformation de modèle de première classe. La différence visible entre observation et intervention (0.656 vs 0.310) prouve que le do() fait un travail structurel (Prong-B).dowhy (identification backdoor/front-door + estimation + refutation), qui relie cette serie a Tweety-11, Infer-5 et ICT-5.Ce notebook traite le do() en probabiliste MCMC + enumeration exacte (pm.do). Le même opérateur structure trois autres series du depot, chacune dans un paradigme différent :
do() et contrefactuels par raisonnement propositionnel (scm.intervene), sans probabilites.do(X_t = x) devient une distribution uniforme p(C) sur les etats ; le do-calculus y mesure a quelle echelle la causalite est la plus effective (EI / CP).P(Cloudy|Rain) vs P(Cloudy|do(Rain)).Synthese transversale : voir le README IIT, section « Ponts causaux : le do-calculus de Pearl a travers les paradigmes », qui aligne l’echelle de Pearl (association -> intervention -> contrefactuel) sur les quatre instanciations de l’opérateur
do.
← PyMC-04-Bayesian-Networks | PyMC-07-Skills-IRT | Retour a la serie