Trois réponses honnêtes et une qui est seule de son côté

Le même coefficient, quatre erreurs-types

Le cours de statistiques a établi que la cantine avait été attribuée à 24 écoles et mesurée sur 1 200 enfants. Voici ce que produit chaque façon de traiter cela.

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

naive = smf.ols("rate ~ feeding_programme", data=per_student).fit()
robust = naive.get_robustcov_results(cov_type="cluster",
                                     groups=per_student["school_id"])

mixed = smf.mixedlm("rate ~ feeding_programme", data=per_student,
                    groups=per_student["school_id"]).fit()

by_school = (per_student.groupby(["school_id", "feeding_programme"])["rate"]
             .mean().reset_index())
aggregated = smf.ols("rate ~ feeding_programme", data=by_school).fit()

library(lme4); library(sandwich); library(lmtest)

naive      <- lm(rate ~ feeding_programme, data = per_student)
robust     <- coeftest(naive, vcov = vcovCL, cluster = ~school_id)
mixed      <- lmer(rate ~ feeding_programme + (1 | school_id), data = per_student)
aggregated <- lm(rate ~ feeding_programme, data = by_school)

Approche	Coefficient	ET	t	n
MCO naïfs	+4,93 pts	0,88	5,63	1 200
ET robustes aux grappes	+4,93 pts	1,46	3,38	1 200
Constante aléatoire	+5,10 pts	1,48	3,44	1 200
MCO au niveau école	+5,21 pts	1,53	3,41	24

Trois des quatre concordent et une non. Les coefficients s’étalent sur 0,3 point ; les trois erreurs-types honnêtes sur 0,07 ; et la naïve est 40 % plus petite que n’importe laquelle d’entre elles.

C’est la forme à attendre. Le regroupement en grappes change rarement beaucoup une estimation et change couramment beaucoup sa précision, et les trois corrections sont trois chemins vers la même destination plutôt que trois réponses concurrentes.

Ce qu’ajoute une constante aléatoire

Les erreurs-types robustes aux grappes corrigent l’inférence et ne disent rien de la structure. Une constante aléatoire l’estime.

print(mixed.summary())
print(f"between-school variance: {float(mixed.cov_re.iloc[0, 0]):.5f}")
print(f"residual variance:       {mixed.scale:.5f}")

VarCorr(mixed)

Modèle	Variance inter-écoles	Variance résiduelle	CCI
Constante seule	0,00142	0,02042	0,065
+ cantine	0,00081	0,02042	0,038

La cantine explique 42,8 % de la variance inter-écoles et rien de la variance intra-école, ce qui est exactement ce que devrait faire un programme de niveau école et vaut d’être vérifié parce que c’est une façon d’attraper un modèle mal spécifié.

Le CCI à constante seule de 0,065 est la même quantité que le cours de statistiques calculait à la main depuis les carrés moyens, avec 0,060. Deux estimateurs, deux réponses légèrement différentes, une conclusion — environ six pour cent de la variation de présence est entre écoles, ce qui, à cinquante enfants par école, suffit à quadrupler la variance d’une comparaison naïve.

Choisir parmi les trois

Les trois corrections ne sont pas interchangeables, et le choix se décide généralement par le nombre de grappes et par ce qui varie à quel niveau.

Employez	Quand	Coût
Agréger à la grappe	Peu de grappes ; l’exposition est de niveau grappe	Une grande école compte autant qu’une petite
ET robustes aux grappes	Beaucoup de grappes (30+) ; vous voulez le modèle individuel	Peu fiables en dessous d’une trentaine de grappes
Constante aléatoire	Vous voulez les composantes de variance, ou des prédicteurs aux deux niveaux	Suppose l’effet aléatoire non corrélé aux prédicteurs

Avec 24 grappes, agrégez ou ajustez la constante aléatoire, et dites laquelle. Le sandwich robuste est l’outil standard et c’est celui qui se comporte le plus mal ici : son asymptotique exige plus de groupes que ce protocole n’en a, et il s’exécutera sans vous prévenir.

Quand l’exposition ne varie qu’entre grappes, les trois convergent, ce que montre le tableau ci-dessus. Elles se séparent quand un prédicteur varie à l’intérieur des grappes — et la constante aléatoire fait alors un travail que les deux autres ne peuvent pas faire.

Des prédicteurs à deux niveaux

C’est là qu’un modèle multiniveau cesse d’être une correction et devient un modèle.

two_level = smf.mixedlm(
    "rate ~ feeding_programme + age_years + disability_reported",
    data=d, groups=d["school_id"]).fit()
print(two_level.summary().tables[1])

lmer(rate ~ feeding_programme + age_years + disability_reported +
       (1 | school_id), data = d)

feeding_programme ne varie qu’entre écoles. age_years et disability_reported varient à l’intérieur. Un seul modèle peut porter les deux, et la constante aléatoire est ce qui le permet : les prédicteurs intra-école sont estimés depuis la variation intra-école et le prédicteur inter-écoles depuis la variation inter-écoles, sans que l’un contamine l’erreur-type de l’autre.

Le mode d’échec est de mettre une exposition de niveau grappe dans un modèle sans terme de grappe. C’est la ligne naïve du premier tableau, et c’est le comportement par défaut de tout logiciel.

Trois niveaux, et quand s’arrêter

Les enfants sont dans des ménages, les ménages dans des villages, les villages dans des districts.

# Two grouping levels: households nested within enumeration areas.
smf.mixedlm("outcome ~ x", data=survey,
            groups=survey["ea_id"], re_formula="1").fit()

lmer(outcome ~ x + (1 | ea_id / household_id), data = survey)

Ajoutez un niveau quand quelque chose est attribué ou mesuré à ce niveau et que vous en avez assez d’unités. Trois niveaux avec six districts au sommet n’estimeront pas une variance de district digne d’être rapportée — le niveau supérieur a besoin d’assez de groupes pour la même raison que la leçon sur les comparaisons multiples avait besoin d’assez de tests.

N’ajoutez pas un niveau par souci d’ordre. Un niveau à quatre groupes ajoute un paramètre inestimable et le faux sentiment d’avoir traité quelque chose.

Rapportez-le en entier

Attendance and school feeding, multilevel analysis

  Linear mixed model, 1,200 students in 24 schools.
  rate ~ feeding_programme + (1 | school_id)

  Feeding programme    +5.10 points   95% CI +2.20 to +8.00

  Between-school SD    0.028      Residual SD  0.143
  ICC 0.038 with the programme term; 0.065 without it, so the programme
  accounts for 43% of the between-school variance.

  The naive single-level model gives +4.93 points with a standard error of
  0.88 rather than 1.48. It is not reported: it treats 50 children in one
  school as 50 independent observations.

  Cluster-robust standard errors on the single-level model (+4.93, SE 1.46)
  and school-level aggregation (+5.21, SE 1.53) give the same answer. With
  24 clusters the robust sandwich is at the edge of its assumptions and is
  reported as a check rather than as the headline.

  Observational. Schools were not randomised into the programme.

Rapporter les trois fait quatre lignes et coupe court à l’objection évidente — que le résultat dépendrait de la méthode. Il n’en dépend pas, et le montrer coûte moins cher que d’en débattre.

La suite

Cette leçon a accordé le modèle à la façon dont le programme a été attribué. La suivante l’accorde à la façon dont l’échantillon a été tiré, sur une enquête aux strates délibérément inégales — et trouve des coefficients qui traversent le seuil de significativité dans les deux sens quand les poids entrent.