Trente-huit cas, pas un rapport de cotes

Reconvertissez le modèle en probabilités

Un coefficient logistique vit sur l’échelle du log des cotes, où personne ne pense. Le remède tient en une ligne et c’est la même à chaque fois : prédire, et moyenner.

import pandas as pd
import numpy as np
import statsmodels.formula.api as smf

model = smf.logit("completed ~ disability + case_category + age_band + sex"
                  " + service_requested + admin1", data=d).fit(disp=0)

everyone_with = d.assign(disability=1)
everyone_without = d.assign(disability=0)

p1 = model.predict(everyone_with).mean()
p0 = model.predict(everyone_without).mean()
print(f"{p1:.4f} vs {p0:.4f}   average marginal effect {p1 - p0:+.4f}")

library(marginaleffects)

avg_comparisons(model, variables = "disability")

26,11 % contre 45,52 %, une différence de −19,41 points.

C’est l’effet marginal moyen : prendre chaque cas des données, demander au modèle ce qu’il prédit si ce cas signalait un handicap, redemander s’il n’en signalait pas, et moyenner l’écart. C’est la réponse du modèle dans l’unité où la question a été posée.

Comparez-la à l’écart brut de la leçon précédente : −19,05 points. Le rapport de cotes ajusté est passé de 0,430 à 0,388 et la différence de risques ajustée n’a presque pas bougé. L’effet marginal moyen est le nombre qui se comporte comme un lecteur attend qu’une estimation ajustée se comporte.

Pourquoi un rapport de cotes fait plusieurs différences de risques

C’est la propriété qui rend l’effet marginal nécessaire plutôt que simplement commode.

profile = dict(case_category="child-protection", age_band="0-11", sex="f",
               admin1="Artibonite")
rows = []
for service in ["health", "psychosocial", "legal", "livelihood-support"]:
    with_d = pd.DataFrame([{**profile, "service_requested": service, "disability": 1}])
    without = pd.DataFrame([{**profile, "service_requested": service, "disability": 0}])
    rows.append((service, model.predict(without)[0], model.predict(with_d)[0]))

print(pd.DataFrame(rows, columns=["service", "no disability", "disability"]).round(3))

predictions(model, newdata = datagrid(service_requested = unique(d$service_requested),
                                      disability = 0:1))

Service demandé	Sans handicap	Handicap signalé	Écart
Santé	69,2 %	46,5 %	−22,7 pts
Psychosocial	64,5 %	41,3 %	−23,2 pts
Juridique	41,0 %	21,2 %	−19,8 pts
Soutien aux moyens d’existence	28,7 %	13,5 %	−15,2 pts

Le rapport de cotes vaut 0,39 sur chacune de ces lignes. Le modèle a été ajusté avec un seul coefficient de handicap, donc par construction le rapport de cotes ne varie pas — et la différence de risques varie tout de même de 15,2 à 23,2 points.

Un rapport de cotes constant n’est pas un effet constant. Là où un résultat est déjà fréquent, le même rapport de cotes déplace plus de points de pourcentage ; là où il est rare, moins. Donc « l’effet du handicap » n’a pas de réponse unique en probabilité sans dire à quel niveau de base — ce que fait précisément l’effet marginal moyen, en moyennant sur les niveaux que la charge de cas possède réellement.

Un intervalle sur l’effet marginal

L’effet marginal est une fonction de tous les coefficients, si bien que son intervalle n’est pas dans le tableau des coefficients. Deux façons de l’obtenir, et la seconde est celle vers laquelle se tourner.

# statsmodels: delta method
margins = model.get_margeff(at="overall")
print(margins.summary())

avg_comparisons(model, variables = "disability")   # delta method, with CI

# bootstrap, when the delta method is awkward or the estimator is custom
rng = np.random.default_rng(20260729)
draws = []
for _ in range(400):
    sample = d.sample(len(d), replace=True, random_state=int(rng.integers(1e9)))
    m = smf.logit(model.model.formula, data=sample).fit(disp=0)
    draws.append(m.predict(sample.assign(disability=1)).mean()
                 - m.predict(sample.assign(disability=0)).mean())
print(np.percentile(draws, [2.5, 97.5]).round(4))

# 400 resamples is enough for a reportable interval; 2,000 for a published one.

−19,41 points, IC à 95 % −26,1 à −13,2 sur 400 rééchantillonnages bootstrap.

Fixez la graine et dites combien de rééchantillonnages. Un intervalle bootstrap qu’on ne peut pas reproduire n’est pas un intervalle, et la règle de cette plateforme sur les nombres produits vaut autant pour ceux que produit un modèle que pour ceux que produit un générateur.

Multipliez-le par la charge de cas

C’est la phrase que lit un responsable de programme, et le cours précédent a établi pourquoi : un effet doit arriver dans l’unité où se prend la décision.

n_disability = int((d["disability"] == 1).sum())
comparator = d[d["disability"] == 0]["completed"].mean()
observed = d[d["disability"] == 1]["completed"].sum()
print(f"{n_disability} cases reporting a disability")
print(f"expected at the comparator rate: {comparator * n_disability:.0f}")
print(f"observed: {int(observed)}  shortfall: {comparator * n_disability - observed:.0f}")

# Three lines, and it is the only line of the analysis a manager will quote.

197 cas signalant un handicap. 90 auraient abouti au taux de comparaison ; 52 l’ont fait. Trente-huit cas de moins.

Trente-huit est un nombre sur lequel une équipe de protection peut agir : cela fait environ trois cas par mois, cela nomme une charge de cas plutôt qu’un pourcentage, et cela se compare à ce que coûterait une correction. « RC 0,39 » ne fait rien de tout cela.

Énoncez l’arithmétique, non le seul résultat. Un lecteur qui voit 197 × 45,5 % − 52 peut vérifier ; un lecteur à qui l’on donne seulement « 38 cas » doit faire confiance.

Quatre façons de se tromper

Prédire à la moyenne au lieu de moyenner les prédictions. Fixer chaque covariable à sa moyenne et prédire une fois donne l’effet pour un cas qui n’existe pas — un cas à 0,6 femme et 0,3 juridique. Moyennez plutôt les prédictions sur les cas réels ; c’est ce que font at="overall" et avg_comparisons.

Rapporter un effet marginal issu d’un modèle avec interaction sans dire où. Si le modèle laisse l’effet du handicap varier selon le service, la moyenne est une moyenne sur une différence réelle et le tableau des profils est la sortie honnête.

Prendre le déficit pour un effet. Trente-huit cas est ce à quoi correspond l’écart observé, non ce que rapporterait sa fermeture. La comparaison reste observationnelle.

Extrapoler au-delà des données. Le modèle prédira volontiers une probabilité d’aboutissement pour un cas juridique de 80 ans dans un département où aucun n’a été enregistré. Prédisez sur des profils que les données contiennent, et dites lesquels.

Rapportez-le en entier

Referral completion by disability status, adjusted

  Logistic model, 1,581 consenting cases with complete covariates.
  Adjusted for case category, age band, sex, service requested, department.

  Average marginal effect   -19.4 points   95% CI -26.1 to -13.2
                                           (400 bootstrap resamples, seed 20260729)
  Predicted completion      26.1% with a disability reported
                            45.5% without

  Applied to the 197 cases reporting a disability, the gap is 38 cases that
  did not reach a service and would have at the comparator rate.

  The gap is not constant: it is 22.7 points for health referrals, where
  completion is common, and 15.2 points for livelihood support, where it is
  not. The odds ratio (0.39) is the same on both.

  Observational. Cases were not randomised, and the model adjusts only for
  what the register records.

L’avant-dernier paragraphe est celui qui empêche de surappliquer la moyenne. Un nombre unique est ce qui sera cité ; la phrase qui dit où il vaut et où il ne vaut pas est ce qui garde la citation honnête.

La suite

Chaque covariable jusqu’ici a été choisie parce qu’elle paraissait pertinente. La leçon suivante en ajoute une qui est pertinente, défendable et fausse — une variable située sur le chemin causal qui retire 42 % de l’écart tout en améliorant chaque statistique d’ajustement.