Stratifier, et voir la pente s'inverser

• Ajustement global. On a

$$\bar{x}=\frac{276}{12}=23,\bar{y}=\frac{222}{12}=18,5,$$

puis, par le calcul détaillé en Python ci-dessous,

$$\sum (x_i-\bar{x}{)}^2=860,\sum (x_i-\bar{x})(y_i-\bar{y})=930,\sum (y_i-\bar{y}{)}^2\approx 1167,4.$$

D'où

$$a=\frac{930}{860}=\frac{93}{86}\approx 1,0814,b=18,5-\frac{93}{86}\times 23\approx -6,3721.$$

Contrôles : $1,0814\times 23-6,3721=18,50=\bar{y}$, la droite passe par $G(23;18,5)$ ; et la somme des résidus, calculée à la question 3, est nulle. Enfin

$$r=\frac{930}{\sqrt{860}\times \sqrt{1167,4}}\approx 0,928.$$

Voilà un ajustement irréprochable, avec une pente franchement positive.

• Ajustement par strate. Prenons les cas légers : $\bar{x}=13$, $\bar{y}=\frac{26}{4}=6,5$, et

$$\begin{array}{lccccc}{x}_i-\bar{x}& -3& -1& 1& 3& 0\\ y_i-\bar{y}& 1,5& 0,4& -0,3& -1,6& 0\\ \text{produit}& -4,5& -0,4& -0,3& -4,8& -10\\ (x_i-\bar{x}{)}^2& 9& 1& 1& 9& 20\end{array}$$

d'où $a=\frac{-10}{20}=-0,5$ et $b=6,5+0,5\times 13=13$. Les valeurs prédites sont $8,0$ ; $7,0$ ; $6,0$ ; $5,0$, donc les résidus valent $0$ ; $-0,1$ ; $+0,2$ ; $-0,1$, de somme nulle.

Les deux autres strates donnent exactement la même pente :

$$\text{modeˊreˊs}:y=-0,5x+30,\text{seˊveˋres}:y=-0,5x+47,$$

avec des résidus $-0,1$ ; $+0,3$ ; $-0,3$ ; $+0,1$ puis $+0,2$ ; $-0,2$ ; $-0,2$ ; $+0,2$, tous de somme nulle, et des droites qui passent par $G_2(23;18,5)$ et $G_3(33;30,5)$.

• Les résidus globaux. Ils valent, dans l'ordre des douze groupes,

$$\begin{array}{lcccc}\text{leˊgers}& +3,56& +0,30& -2,57& -6,03\\ \text{modeˊreˊs}& +4,64& +1,88& -1,88& -4,64\\ \text{seˊveˋres}& +6,13& +2,57& -0,60& -3,36\end{array}$$

Leur somme est nulle, comme il se doit — mais ils ne sont pas répartis au hasard : à l'intérieur de chaque strate ils décroissent régulièrement, puis remontent brutalement au passage à la strate suivante. Ce motif en dents de scie est la signature d'une variable oubliée. Le coefficient $r\approx 0,93$ ne l'a pas vu : il mesure la dispersion, pas la structure.

• La pente qui répond. Celle des strates, $-0,5$. La question « faut-il augmenter la dose ? » est une question d'intervention : elle porte sur un patient donné, dont la gravité est ce qu'elle est. C'est donc à gravité fixée qu'il faut comparer — et à gravité fixée, augmenter la dose fait baisser la mortalité de $0,5$ point par milligramme. La pente globale $+1,08$ répond, elle, à la question d'observation : « si l'on m'annonce une dose de $30$ mg, que dois-je craindre ? » — et la réponse, le pire, est juste : cette dose signale un cas grave.

def moyenne(v):
    return sum(v) / len(v)

def ajustement(x, y):
    xb, yb = moyenne(x), moyenne(y)
    sxy = sum((x[i] - xb) * (y[i] - yb) for i in range(len(x)))
    sxx = sum((t - xb) ** 2 for t in x)
    a = sxy / sxx
    b = yb - a * xb
    residus = [y[i] - (a * x[i] + b) for i in range(len(x))]
    assert abs(a * xb + b - yb) < 1e-9          # passe par le point moyen
    assert abs(sum(residus)) < 1e-9             # residus de somme nulle
    return a, b, residus

strates = [([10, 12, 14, 16], [8.0, 6.9, 6.2, 4.9]),
           ([20, 22, 24, 26], [19.9, 19.3, 17.7, 17.1]),
           ([30, 32, 34, 36], [32.2, 30.8, 29.8, 29.2])]

for x, y in strates:
    a, b, e = ajustement(x, y)
    print(round(a, 6), round(b, 4))     # -0.5 13.0 / -0.5 30.0 / -0.5 47.0

xg = [t for x, y in strates for t in x]
yg = [t for x, y in strates for t in y]
a, b, e = ajustement(xg, yg)
print(round(a, 6), round(b, 4))         # 1.081395 -6.3721
print([round(t, 2) for t in e])         # motif decroissant dans chaque strate