Ajustez une régression linéaire simple

Jusqu’ici, nous avons étudié une variable à la fois. Mais la question la plus fréquente en data science porte sur un lien : comment le loyer dépend-il de la surface ? La distance de freinage de la vitesse ? Le nombre de lettres à trier du poids du courrier ?

La régression linéaire simple modélise ce lien par une droite, tout en tenant compte du bruit. Elle réutilise tout ce que vous avez vu : estimation, lois de Student et de Fisher, intervalles et tests. Notre fil rouge : le loyer d’appartements rennais en fonction de leur surface.

Le fil rouge et le modèle

Voici 10 annonces de T1 et T2 dans un même quartier :

Surface $x$ (m²)	18	25	30	34	40	45	52	60	68	75
Loyer $y$ (€)	430	480	590	560	700	660	810	820	990	980

Le nuage de points monte clairement, à peu près en ligne droite. On pose le modèle de régression linéaire simple :

$$Y_i = \beta_0 + \beta_1 x_i + \varepsilon_i, \qquad i = 1, \dots, n$$

• $Y_i$ est la variable à expliquer (ou réponse), aléatoire : le loyer ;
• $x_i$ est la variable explicative, considérée comme fixée : la surface ;
• $\beta_0$ (ordonnée à l’origine) et $\beta_1$ (pente) sont les paramètres inconnus ;
• $\varepsilon_i$ est l’erreur : tout ce que la surface n’explique pas (étage, état, luminosité…).

Les hypothèses portent sur les erreurs :

1. $E[\varepsilon_i] = 0$ : le modèle est correct en moyenne (la relation est bien linéaire) ;
2. $V(\varepsilon_i) = \sigma^2$ pour tout $i$ : la variance est constante, on parle d’homoscédasticité ;
3. les $\varepsilon_i$ sont non corrélées (indépendantes) ;
4. pour les tests et intervalles, les $\varepsilon_i$ sont gaussiennes : $\varepsilon_i \sim \mathcal{N}(0, \sigma^2)$.

Estimer la droite par les moindres carrés

Quelle droite choisir ? On cherche $\beta_0$ et $\beta_1$ qui minimisent la somme des carrés des écarts verticaux entre les points et la droite :

$$\min_{\beta_0, \beta_1} \sum_{i=1}^{n} \left( y_i - \beta_0 - \beta_1 x_i \right)^2$$

En annulant les deux dérivées partielles, on obtient les estimateurs des moindres carrés :

$$\hat\beta_1 = \frac{\sum_{i}(x_i - \bar x)(Y_i - \bar Y)}{\sum_{i}(x_i - \bar x)^2}, \qquad \hat\beta_0 = \bar Y - \hat\beta_1 \bar x$$

La droite passe donc toujours par le point moyen $(\bar x, \bar y)$. Sous l’hypothèse gaussienne, ces estimateurs coïncident avec ceux du maximum de vraisemblance (chapitre 3) : maximiser la vraisemblance gaussienne revient à minimiser la somme des carrés.

Exemple

Avec $\bar x = 44{,}7$ m², $\bar y = 702$ €, $\sum (x_i - \bar x)^2 = 3\,202{,}1$ et $\sum (x_i - \bar x)(y_i - \bar y) = 32\,526$ :

$$\hat\beta_1 = \frac{32\,526}{3\,202{,}1} \approx 10{,}16 \text{ euros par m}^2, \qquad \hat\beta_0 = 702 - 10{,}16 \times 44{,}7 \approx 247{,}95 \text{ euros}$$

Chaque mètre carré supplémentaire augmente le loyer de 10,16 € en moyenne. L’ordonnée à l’origine (248 € pour 0 m²) n’a pas de sens physique : elle sert seulement à placer la droite dans la zone observée.

La droite des moindres carrés rend la somme des carrés des segments roses la plus petite possible, et passe toujours par le point moyen.

Les estimateurs $\hat\beta_0$ et $\hat\beta_1$ sont sans biais, et la variance de la pente vaut

$$V(\hat\beta_1) = \frac{\sigma^2}{\sum_i (x_i - \bar x)^2}$$

La pente est d’autant mieux estimée que les $x_i$ sont étalés : pour étudier l’effet de la surface, mieux vaut des studios et des grands T2 que dix appartements de 40 m².

Même bruit, même modèle : avec des surfaces étalées, la pente varie peu d'un échantillon à l'autre (écart-type 0,72 €/m²) ; avec des surfaces groupées autour de 40 m², elle devient presque indéterminée.

Résidus et variance résiduelle

Les valeurs ajustées sont $\hat y_i = \hat\beta_0 + \hat\beta_1 x_i$, et les résidus $\hat\varepsilon_i = y_i - \hat y_i$ (leur somme est nulle). La variance $\sigma^2$ des erreurs s’estime sans biais par

$$\hat\sigma^2 = \frac{1}{n - 2} \sum_{i=1}^{n} \hat\varepsilon_i^2 = \frac{SCR}{n - 2}$$

où $SCR$ est la somme des carrés résiduelle. On divise par $n - 2$ car on a estimé deux paramètres : c’est le même phénomène que le $n - 1$ de la variance corrigée. De plus, $(n - 2)\hat\sigma^2/\sigma^2 \sim \chi^2_{n-2}$.

Ici, $SCR \approx 13\,170$, donc $\hat\sigma^2 \approx 1\,646$ et $\hat\sigma \approx 40{,}6$ € : à surface égale, les loyers fluctuent d’une quarantaine d’euros autour de la droite.

Tester l’effet de la variable explicative

La surface a-t-elle vraiment un effet sur le loyer ? On teste $H_0 : \beta_1 = 0$ (pas de lien linéaire) contre $H_1 : \beta_1 \ne 0$. En remplaçant $\sigma$ par $\hat\sigma$ dans l’écart-type de $\hat\beta_1$, on obtient l’erreur-type $\hat\sigma_{\hat\beta_1} = \hat\sigma/\sqrt{\sum (x_i - \bar x)^2}$, et

$$T_0 = \frac{\hat\beta_1}{\hat\sigma_{\hat\beta_1}} \sim \mathcal{T}_{n-2} \text{ sous } H_0$$

Ici, $\hat\sigma_{\hat\beta_1} = 40{,}57/\sqrt{3\,202{,}1} \approx 0{,}717$ et $t_0 = 10{,}16/0{,}717 \approx 14{,}17$, à comparer à $t_{8;0,975} \approx 2{,}306$. On rejette très largement $H_0$ (p-valeur de l’ordre de $6 \times 10^{-7}$). L’IC à 95 % de la pente vaut $10{,}16 \pm 2{,}306 \times 0{,}717$, soit environ $[8{,}50 \ ;\ 11{,}81]$ €/m².

La décomposition de la variabilité et le R²

La variabilité totale des loyers se décompose en une part expliquée par la droite et une part résiduelle :

$$\underbrace{\sum_i (y_i - \bar y)^2}_{SCT} = \underbrace{\sum_i (\hat y_i - \bar y)^2}_{SCM} + \underbrace{\sum_i (y_i - \hat y_i)^2}_{SCR}$$

$SCT$ est la somme des carrés totale, $SCM$ celle du modèle (expliquée), $SCR$ la résiduelle. Le coefficient de détermination mesure la part expliquée :

$$R^2 = \frac{SCM}{SCT} = 1 - \frac{SCR}{SCT}$$

Il est compris entre 0 et 1 et, en régression simple, il vaut le carré du coefficient de corrélation linéaire $r$ entre $x$ et $y$. Ici, $R^2 = 1 - 13\,170/343\,560 \approx 0{,}962$ : la surface explique 96 % de la variabilité des loyers de l’échantillon.

Pour chaque appartement, l'écart à ȳ se coupe en une part expliquée par la droite et un résidu ; en sommant les carrés, SCT = SCM + SCR, et la part expliquée donne R² = 0,962.

Le test de Fisher global compare la part expliquée à la part résiduelle :

$$F_0 = \frac{SCM / 1}{SCR / (n - 2)} \sim \mathcal{F}_{1, n-2} \text{ sous } H_0 : \beta_1 = 0$$

Ici $F_0 \approx 200{,}7 = 14{,}17^2$ : en régression simple, ce test est identique au test de Student sur la pente (le carré d’une $\mathcal{T}_{n-2}$ est une $\mathcal{F}_{1,n-2}$).

logements <- data.frame(
  surface = c(18, 25, 30, 34, 40, 45, 52, 60, 68, 75),
  loyer   = c(430, 480, 590, 560, 700, 660, 810, 820, 990, 980)
)
reg <- lm(loyer ~ surface, data = logements)
summary(reg)    # coefficients, erreurs-types, t, p-valeurs, sigma, R², F
confint(reg)    # IC à 95 % des coefficients

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

logements = pd.DataFrame({
    "surface": [18, 25, 30, 34, 40, 45, 52, 60, 68, 75],
    "loyer":   [430, 480, 590, 560, 700, 660, 810, 820, 990, 980],
})
reg = smf.ols("loyer ~ surface", data=logements).fit()
print(reg.summary())    # même tableau que summary() en R

Piège

Un $R^2$ élevé ne prouve pas que le modèle est bon, et une pente significative ne prouve pas une causalité. Le célèbre quartet d’Anscombe contient quatre jeux de 11 points ayant tous la même droite ($\hat y \approx 3{,}00 + 0{,}500\,x$) et le même $R^2 \approx 0{,}67$ : l’un est bien linéaire, un autre est une parabole, un troisième est une droite parfaite perturbée par un point aberrant, le dernier tient entièrement à un seul point extrême. Seul un graphique permet de les distinguer. En TP, un exercice reposait précisément sur ce piège.

Quatre jeux, une seule droite et un même R² ≈ 0,67. Seul le premier est bien décrit par le modèle : sans graphique, impossible de le savoir.

Prévoir : intervalle de confiance ou intervalle de prévision ?

Pour un appartement de $x_0 = 40$ m², la droite prévoit $\hat y_0 = 247{,}95 + 10{,}158 \times 40 \approx 654{,}3$ €. Mais cette prévision peut répondre à deux questions différentes.

Le loyer moyen de tous les 40 m² du quartier, $\beta_0 + \beta_1 x_0$ : c’est un paramètre, et son IC tient compte de l’incertitude sur la droite seulement :

$$\hat y_0 \pm t_{n-2;1-\alpha/2}\ \hat\sigma \sqrt{\frac{1}{n} + \frac{(x_0 - \bar x)^2}{\sum_j (x_j - \bar x)^2}}$$

Le loyer d’un appartement précis de 40 m², $Y_0 = \beta_0 + \beta_1 x_0 + \varepsilon_0$ : il faut ajouter la variabilité propre de cet appartement, d’où un « $1 +$ » sous la racine. C’est l’intervalle de prévision :

$$\hat y_0 \pm t_{n-2;1-\alpha/2}\ \hat\sigma \sqrt{1 + \frac{1}{n} + \frac{(x_0 - \bar x)^2}{\sum_j (x_j - \bar x)^2}}$$

Exemple

Pour $x_0 = 40$ m², à 95 % :

• IC du loyer moyen : $654{,}3 \pm 30{,}6$, soit $[623{,}7 \ ;\ 684{,}9]$ € ;
• intervalle de prévision d’un loyer : $654{,}3 \pm 98{,}4$, soit $[555{,}8 \ ;\ 752{,}7]$ €.

L’intervalle de prévision est trois fois plus large : même si l’on connaissait parfaitement la droite, un appartement particulier s’en écarterait. Et pour un 100 m², hors de la zone observée, la demi-largeur de prévision monte à 134 € : le terme $(x_0 - \bar x)^2$ pénalise l’extrapolation, qui reste de toute façon hasardeuse (rien ne dit que la relation reste linéaire au-delà de 75 m²).

Les deux bandes s'évasent loin du point moyen. La bande de prévision reste large même là où la droite est bien connue : un appartement précis s'écarte toujours de la moyenne.

nouveau <- data.frame(surface = 40)
predict(reg, nouveau, interval = "confidence")   # loyer moyen
predict(reg, nouveau, interval = "prediction")   # loyer d'un appartement

Valider le modèle : regarder les résidus

Toutes ces conclusions reposent sur les hypothèses du modèle. On les vérifie graphiquement sur les résidus.

• Résidus en fonction des valeurs ajustées $\hat y_i$ : le nuage doit être sans structure, centré sur 0, d’épaisseur constante. On trace souvent les droites $\pm 2\hat\sigma$ : environ 95 % des résidus doivent être entre les deux.
• Une courbure (résidus positifs aux extrémités, négatifs au milieu) signale que la relation n’est pas linéaire.
• Un entonnoir (résidus de plus en plus dispersés) signale une variance non constante.
• Une droite de Henry (QQ-plot) des résidus vérifie la normalité.

Pour nos loyers, les résidus vont de $-45$ à $+51$ €, tous dans la bande $\pm 2\hat\sigma \approx \pm 81$ €, sans forme particulière.

En vert, les loyers : résidus sans forme et alignés sur la droite de Henry. En rose, deux cas simulés à reconnaître : une courbure (relation non linéaire) et un entonnoir (variance non constante).

Quand un graphique révèle un problème, on transforme les variables : $\ln Y$ ou $\sqrt{Y}$ stabilisent une variance qui croît avec la moyenne ; $x^2$ ou $\ln x$ capturent une courbure. Le modèle reste « linéaire » au sens où il est linéaire en les paramètres. En TP, nous avons ainsi étudié la distance de freinage d’un véhicule en fonction de sa vitesse : la physique (l’énergie cinétique) suggère une distance proportionnelle au carré de la vitesse, et c’est le graphe des résidus qui permet de juger si une simple droite suffit ou s’il faut régresser sur ce carré.

À retenir

Régression simple : $Y_i = \beta_0 + \beta_1 x_i + \varepsilon_i$, estimée par moindres carrés. On teste $\beta_1 = 0$ par Student ($\mathcal{T}_{n-2}$), on mesure l’ajustement par $R^2 = 1 - SCR/SCT$, on distingue IC de la moyenne et intervalle de prévision (avec « $1 +$ »), et on valide toujours par le graphe des résidus.

En entretien

« Les moindres carrés minimisent la somme des carrés des résidus ; sous hypothèse gaussienne, c’est aussi le maximum de vraisemblance. Avant de lire les p-valeurs, je vérifie les résidus : linéarité, variance constante, indépendance et normalité. »

En résumé

• Modèle : $Y_i = \beta_0 + \beta_1 x_i + \varepsilon_i$, erreurs centrées, de variance constante, non corrélées, gaussiennes pour l’inférence.
• Moindres carrés : $\hat\beta_1 = \sum(x_i - \bar x)(y_i - \bar y)/\sum(x_i - \bar x)^2$ et $\hat\beta_0 = \bar y - \hat\beta_1 \bar x$ ; $\hat\sigma^2 = SCR/(n-2)$.
• Test de la pente : $\hat\beta_1/\hat\sigma_{\hat\beta_1} \sim \mathcal{T}_{n-2}$ ; en simple, il équivaut au test de Fisher global.
• $SCT = SCM + SCR$ et $R^2 = SCM/SCT$, égal à $r^2$.
• Intervalle de prévision plus large que l’IC de la moyenne ; les deux s’élargissent loin de $\bar x$.
• Toujours tracer les résidus ; en cas de courbure ou d’entonnoir, transformer les variables.

Exercices

Exercice · Calculer une droite de régression

facile

Sur 5 jours, un glacier relève la température maximale $x$ (°C) et ses ventes $y$ (centaines de cornets) : (20 ; 3), (22 ; 4), (25 ; 6), (28 ; 7), (30 ; 10). Calculez $\hat\beta_1$ et $\hat\beta_0$, puis la vente prévue à 26 °C.

Voir la solution

$\bar x = 125/5 = 25$ et $\bar y = 30/5 = 6$. Les écarts $x_i - \bar x$ valent $-5, -3, 0, 3, 5$ et les écarts $y_i - \bar y$ valent $-3, -2, 0, 1, 4$.

• $\sum (x_i - \bar x)(y_i - \bar y) = 15 + 6 + 0 + 3 + 20 = 44$ ;
• $\sum (x_i - \bar x)^2 = 25 + 9 + 0 + 9 + 25 = 68$.

Donc $\hat\beta_1 = 44/68 \approx 0{,}647$ centaine de cornets par degré, et $\hat\beta_0 = 6 - 0{,}647 \times 25 \approx -10{,}18$. À 26 °C : $\hat y = -10{,}18 + 0{,}647 \times 26 \approx 6{,}65$, soit environ 665 cornets.

En code

Python

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

glaces = pd.DataFrame({"temperature": [20, 22, 25, 28, 30],
                       "ventes": [3, 4, 6, 7, 10]})
reg = smf.ols("ventes ~ temperature", data=glaces).fit()
print(reg.params)                                          # Intercept -10.18, temperature 0.647
print(reg.predict(pd.DataFrame({"temperature": [26]})))    # 6.65 centaines de cornets

R

glaces <- data.frame(temperature = c(20, 22, 25, 28, 30),
                     ventes      = c(3, 4, 6, 7, 10))
reg <- lm(ventes ~ temperature, data = glaces)
coef(reg)                                            # (Intercept) -10.18, temperature 0.647
predict(reg, newdata = data.frame(temperature = 26)) # 6.65 centaines de cornets

Exercice · Lire une sortie de R

moyen

Une régression du temps de trajet (minutes) sur la distance (km), sur $n = 27$ trajets, donne : pente estimée 1,85 avec une erreur-type de 0,40, et $R^2 = 0{,}46$.

1. La distance a-t-elle un effet significatif au niveau 5 % ($t_{25;0,975} \approx 2{,}060$) ?
2. Que vaut la statistique du test de Fisher global ?
3. Que penser du modèle ?

Voir la solution

1. $t_0 = 1{,}85/0{,}40 = 4{,}625 > 2{,}060$ : on rejette $\beta_1 = 0$, l’effet est significatif.
2. En régression simple, $F_0 = t_0^2 \approx 21{,}4$, à comparer à une loi $\mathcal{F}_{1,25}$. On peut le retrouver par $F_0 = \frac{R^2/1}{(1 - R^2)/(n-2)} = \frac{0{,}46}{0{,}54/25} \approx 21{,}3$ (l’écart vient des arrondis).
3. L’effet de la distance est net (chaque kilomètre ajoute environ 1,85 minute), mais la distance n’explique que 46 % de la variabilité : d’autres facteurs (trafic, heure, mode de transport) comptent beaucoup. Significatif ne veut pas dire suffisant pour prévoir précisément un trajet.

En code

Python

from scipy import stats

n, pente, erreur_type, R2 = 27, 1.85, 0.40, 0.46
t0 = pente / erreur_type
print(t0, stats.t.ppf(0.975, n - 2), 2 * stats.t.sf(t0, n - 2))   # 4.625 > 2.060, p = 0.0001

F_par_t = t0**2                                    # F = t^2 en régression simple
F_par_R2 = (R2 / 1) / ((1 - R2) / (n - 2))         # à partir du R²
print(F_par_t, F_par_R2, stats.f.sf(F_par_t, 1, n - 2))   # 21.4  21.3  p = 0.0001

R

n <- 27; pente <- 1.85; erreur_type <- 0.40; R2 <- 0.46
t0 <- pente / erreur_type
t0; qt(0.975, n - 2); 2 * pt(-t0, n - 2)    # 4.625 > 2.060, p = 0.0001

F_par_t  <- t0^2                            # F = t^2 en régression simple
F_par_R2 <- (R2 / 1) / ((1 - R2) / (n - 2)) # à partir du R²
F_par_t; F_par_R2; 1 - pf(F_par_t, 1, n - 2)   # 21.4  21.3  p = 0.0001

Exercice · Linéariser un modèle exponentiel

difficile

Le nombre d’utilisateurs d’une application semble croître exponentiellement : $U(t) = a\,e^{b t}$. Comment estimer $a$ et $b$ par une régression linéaire simple ? Quelle hypothèse cela suppose-t-il sur les erreurs ?

Voir la solution

En passant au logarithme : $\ln U(t) = \ln a + b\,t$. On régresse donc $Y_i = \ln u_i$ sur $t_i$ : le modèle $Y_i = \beta_0 + \beta_1 t_i + \varepsilon_i$ donne $\hat b = \hat\beta_1$ et $\hat a = e^{\hat\beta_0}$.

Cela suppose que les erreurs sont additives et de variance constante sur l’échelle logarithmique, donc multiplicatives sur l’échelle d’origine : $U_i = a\,e^{b t_i}\,e^{\varepsilon_i}$. Autrement dit, l’écart typique est un pourcentage constant du nombre d’utilisateurs, ce qui est souvent réaliste pour une croissance. On le vérifie sur le graphe des résidus du modèle en logarithme.

Attention : $e^{\hat\beta_0 + \hat\beta_1 t}$ estime la médiane de $U(t)$ plutôt que sa moyenne, car l’exponentielle d’une moyenne n’est pas la moyenne de l’exponentielle.

En code

Python

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

# Croissance simulée avec erreur multiplicative : U = a exp(b t) exp(eps)
rng = np.random.default_rng(0)
a, b = 500, 0.15
semaines = np.arange(1, 31)
U = a * np.exp(b * semaines) * np.exp(rng.normal(0, 0.2, size=30))
donnees = pd.DataFrame({"t": semaines, "U": U})

reg = smf.ols("np.log(U) ~ t", data=donnees).fit()     # régression sur le log
print(np.exp(reg.params["Intercept"]), reg.params["t"])   # a_hat = 493, b_hat = 0.149 (vrais : 500 et 0.15)

# exp(prévision) estime la médiane ; la moyenne vaut médiane * exp(sigma^2 / 2)
print(np.exp(0.2**2 / 2))                              # 1.020 : écart de 2 %

R

# Croissance simulée avec erreur multiplicative : U = a exp(b t) exp(eps)
set.seed(0)
a <- 500; b <- 0.15
semaines <- 1:30
U <- a * exp(b * semaines) * exp(rnorm(30, 0, 0.2))

reg <- lm(log(U) ~ semaines)                 # régression sur le log
exp(coef(reg)[1]); coef(reg)[2]              # a_hat = 531, b_hat = 0.146 (vrais : 500 et 0.15)
plot(fitted(reg), resid(reg))                # résidus : pas de structure attendue

# exp(prévision) estime la médiane ; la moyenne vaut médiane * exp(sigma^2 / 2)
exp(0.2^2 / 2)                               # 1.020 : écart de 2 %