Passez à la régression linéaire multiple

La surface explique bien le loyer d’un appartement, mais pas tout : à surface égale, un logement près du centre se loue plus cher. Et l’étage ? Et la présence d’un balcon ? Pour tenir compte de plusieurs facteurs à la fois, on passe à la régression linéaire multiple.

Le principe est le même qu’en régression simple, mais deux questions nouvelles apparaissent : comment interpréter un coefficient quand les variables sont liées entre elles, et comment choisir les variables à garder. Notre fil rouge : 12 appartements rennais décrits par trois variables.

Le fil rouge et le modèle

Surface (m²)	18	25	30	34	40	45	52	60	68	75	28	55
Distance au centre (km)	0,5	2,0	1,0	3,5	1,5	4,0	2,5	1,0	3,0	5,0	4,5	0,8
Étage	3	1	4	0	2	5	1	3	0	2	4	1
Loyer (€)	470	470	600	500	700	610	790	900	930	900	450	880

Le modèle de régression linéaire multiple à $p$ variables explicatives s’écrit

$$Y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \dots + \beta_p x_{ip} + \varepsilon_i, \qquad \varepsilon_i \text{ i.i.d. } \mathcal{N}(0, \sigma^2)$$

avec les mêmes hypothèses qu’en régression simple. Ici, $p = 3$ : la surface $x_1$, la distance $x_2$ et l’étage $x_3$.

L’écriture matricielle

Avec plusieurs variables, les formules deviennent lisibles sous forme matricielle :

$$Y = X\beta + \varepsilon$$

• $Y$ est le vecteur des $n$ réponses (les loyers) ;
• $X$ est la matrice du plan d’expérience, de taille $n \times (p+1)$ : une première colonne de 1 (pour la constante $\beta_0$), puis une colonne par variable explicative ;
• $\beta = (\beta_0, \dots, \beta_p)^\top$ est le vecteur des paramètres ;
• $\varepsilon \sim \mathcal{N}_n(0, \sigma^2 I_n)$ est le vecteur des erreurs, $I_n$ désignant la matrice identité.

Minimiser $\lVert Y - X\beta \rVert^2$ conduit aux équations normales $X^\top X \hat\beta = X^\top Y$, d’où, si $X^\top X$ est inversible :

$$\hat\beta = (X^\top X)^{-1} X^\top Y$$

Cet estimateur est sans biais et suit la loi $\mathcal{N}_{p+1}\left(\beta,\ \sigma^2 (X^\top X)^{-1}\right)$ : les éléments diagonaux de $\sigma^2 (X^\top X)^{-1}$ donnent les variances des coefficients. La variance des erreurs s’estime par

$$\hat\sigma^2 = \frac{SCR}{n - p - 1}$$

car on a estimé $p + 1$ paramètres. En TP, nous avons d’ailleurs calculé $(X^\top X)^{-1} X^\top y$ « à la main » avec solve() et t() en R, puis vérifié que lm() donnait les mêmes coefficients.

Les moindres carrés cherchent le point du plan des Xβ le plus proche de Y : c'est sa projection orthogonale. Les équations normales disent exactement que le résidu est orthogonal à chaque colonne de X.

logements <- data.frame(
  surface  = c(18, 25, 30, 34, 40, 45, 52, 60, 68, 75, 28, 55),
  distance = c(0.5, 2.0, 1.0, 3.5, 1.5, 4.0, 2.5, 1.0, 3.0, 5.0, 4.5, 0.8),
  etage    = c(3, 1, 4, 0, 2, 5, 1, 3, 0, 2, 4, 1),
  loyer    = c(470, 470, 600, 500, 700, 610, 790, 900, 930, 900, 450, 880)
)
X <- cbind(1, as.matrix(logements[, c("surface", "distance", "etage")]))
y <- logements$loyer
solve(t(X) %*% X) %*% t(X) %*% y        # beta chapeau, calculé à la main

complet <- lm(loyer ~ surface + distance + etage, data = logements)
summary(complet)                         # mêmes coefficients

Lire et tester les coefficients

La sortie de summary(complet) donne :

Coefficient	Estimation	Erreur-type	$t$	p-valeur
Constante	299,6	23,9	12,53	$1{,}5 \times 10^{-6}$
Surface	11,01	0,43	25,72	$5{,}6 \times 10^{-9}$
Distance	$-43{,}16$	4,78	$-9{,}03$	$1{,}8 \times 10^{-5}$
Étage	1,34	4,49	0,30	0,77

avec $\hat\sigma \approx 23{,}1$ €, $R^2 \approx 0{,}989$ et $F \approx 244{,}5$ sur $(3, 8)$ degrés de liberté.

Interpréter « toutes choses égales par ailleurs »

Le coefficient de la distance, $-43{,}16$, se lit ainsi : à surface et étage égaux, chaque kilomètre supplémentaire depuis le centre fait baisser le loyer de 43 € en moyenne. C’est un effet partiel, qui neutralise les autres variables du modèle.

Cette nuance est essentielle. Regardez la distance seule : sa corrélation avec le loyer vaut $-0{,}03$, et la régression simple du loyer sur la distance donne une pente non significative (p-valeur 0,91). Pourtant, dans le modèle multiple, son effet est net. Pourquoi ? Parce que dans cet échantillon, les grands appartements sont plutôt loin du centre (corrélation de 0,30 entre surface et distance) : l’effet positif de la surface masque l’effet négatif de la distance tant qu’on ne les sépare pas.

Seule, la distance semble sans effet : les grands appartements (gros points) sont souvent loin du centre. Une fois l'effet de la surface retiré, la baisse de 43 € par kilomètre apparaît nettement.

Piège

Un coefficient de régression multiple ne se lit jamais seul : il dépend des autres variables présentes dans le modèle. Ajouter ou retirer une variable corrélée peut changer sa valeur, voire son signe. Et une variable non significative seule peut devenir essentielle à côté des autres, comme la distance ici.

Le test de Student sur chaque coefficient

Pour chaque $j$, on teste $H_0 : \beta_j = 0$ (« la variable $x_j$ n’apporte rien, les autres étant dans le modèle ») par

$$T_j = \frac{\hat\beta_j}{\hat\sigma_{\hat\beta_j}} \sim \mathcal{T}_{n-p-1} \text{ sous } H_0$$

où $\hat\sigma_{\hat\beta_j}$ est l’erreur-type, racine du $j$-ème élément diagonal de $\hat\sigma^2 (X^\top X)^{-1}$. Ici, avec $n - p - 1 = 8$ degrés de liberté, la surface et la distance sont très significatives, l’étage non (p-valeur 0,77).

Le test de Fisher global

Le modèle sert-il à quelque chose ? On teste $H_0 : \beta_1 = \dots = \beta_p = 0$ (aucune variable n’a d’effet) contre $H_1$ : au moins un $\beta_j$ est non nul :

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

Ici $F_0 \approx 244{,}5$ pour un seuil $f_{3,8;0,95} \approx 4{,}07$ : on rejette très largement. C’est la dernière ligne de summary().

Attention

Le test global et les tests individuels ne répondent pas à la même question. Avec des variables très corrélées entre elles (on parle de colinéarité), le test global peut être très significatif alors qu’aucun test individuel ne l’est : chaque variable « prend la place » de l’autre. Les erreurs-types gonflent, les coefficients deviennent instables. En TP, sur des données de composition du lait (densité, taux protéique, taux de caséine…), il fallait justement repérer ces fortes corrélations avant d’interpréter les coefficients.

Comparer des modèles : R² ajusté et AIC

Faut-il garder l’étage ? Le $R^2$ ne peut pas trancher : il augmente toujours quand on ajoute une variable, même inutile, car la somme des carrés résiduelle ne peut que diminuer. Il faut des critères qui pénalisent la complexité.

Le $R^2$ ajusté corrige le $R^2$ par les degrés de liberté :

$$R^2_{aj} = 1 - \frac{SCR / (n - p - 1)}{SCT / (n - 1)}$$

Il peut baisser quand on ajoute une variable inutile.

Le critère d’information d’Akaike (AIC) vaut $-2 \ln \hat L + 2k$, où $\hat L$ est la vraisemblance maximale et $k$ le nombre de paramètres. Il récompense l’ajustement (premier terme) et pénalise la complexité (second terme) : plus il est petit, mieux c’est. Pour un modèle linéaire gaussien, R l’évalue à une constante près par $n \ln(SCR/n) + 2k$ (fonction extractAIC).

Modèle	$R^2$	$R^2_{aj}$	AIC
surface	0,878	0,866	103,65
surface + distance	0,9891	0,9867	76,65
surface + distance + étage	0,9892	0,9852	78,51

Le $R^2$ préfère le modèle complet, de très peu. Le $R^2$ ajusté et l’AIC préfèrent tous deux le modèle sans l’étage.

Ajouter l'étage fait encore monter le R² (de 0,0001), mais le R² ajusté baisse et l'AIC remonte : la variable coûte plus qu'elle n'apporte.

On peut aussi comparer deux modèles emboîtés par un test de Fisher partiel : anova(petit, grand) teste si les variables ajoutées réduisent significativement la $SCR$. Pour l’étage, on trouve $F = 0{,}089$, soit exactement le carré du $t = 0{,}30$ de la table : retirer une seule variable, c’est la même chose que son test de Student.

Sélectionner les variables avec step

Avec 3 variables, on peut comparer les $2^3 = 8$ sous-modèles (recherche exhaustive). Avec 30, il y en a plus d’un milliard. On utilise alors des procédures pas à pas :

• descendante (backward) : on part du modèle complet et on retire, à chaque étape, la variable dont le retrait diminue le plus l’AIC, tant que c’est possible ;
• ascendante (forward) : on part du modèle vide et on ajoute la variable qui diminue le plus l’AIC ;
• pas à pas (both) : à chaque étape, on autorise ajouts et retraits.

Exemple

Procédure descendante sur nos logements :

1. Modèle complet : AIC = 78,51. Retirer l’étage donne 76,65 ; retirer la distance, 105,50 ; retirer la surface, 129,64. On retire l’étage.
2. Modèle surface + distance : AIC = 76,65. Retirer la distance donne 103,65 ; retirer la surface, 128,85. Aucun retrait n’améliore l’AIC : on s’arrête.

Le modèle retenu est $\widehat{\text{loyer}} = 303{,}8 + 10{,}97 \times \text{surface} - 42{,}95 \times \text{distance}$. Pour un 40 m² à 2 km du centre, il prévoit $656{,}6$ €, avec un intervalle de prévision à 95 % de $\pm 51{,}8$ €.

À chaque étape, step essaie tous les retraits possibles et garde celui qui fait le plus baisser l'AIC (en vert) ; il s'arrête quand tous les retraits le font monter (en rose).

retenu <- step(complet, direction = "backward")   # affiche chaque étape et l'AIC
summary(retenu)
extractAIC(complet); extractAIC(retenu)
predict(retenu, data.frame(surface = 40, distance = 2), interval = "prediction")

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

logements = pd.DataFrame({
    "surface":  [18, 25, 30, 34, 40, 45, 52, 60, 68, 75, 28, 55],
    "distance": [0.5, 2.0, 1.0, 3.5, 1.5, 4.0, 2.5, 1.0, 3.0, 5.0, 4.5, 0.8],
    "etage":    [3, 1, 4, 0, 2, 5, 1, 3, 0, 2, 4, 1],
    "loyer":    [470, 470, 600, 500, 700, 610, 790, 900, 930, 900, 450, 880],
})
for formule in ["loyer ~ surface + distance + etage", "loyer ~ surface + distance"]:
    m = smf.ols(formule, data=logements).fit()
    print(formule, round(m.rsquared_adj, 4), round(m.aic, 2))   # AIC complet (constante différente de R)

Attention

L’AIC de statsmodels (et de la fonction AIC() de R) inclut toutes les constantes de la vraisemblance, contrairement à extractAIC, qui est utilisé par step. Les valeurs diffèrent, mais les écarts entre modèles, et donc le classement, sont identiques.

En TP, nous avons suivi cette démarche sur un jeu de données marketing : les ventes d’un produit en fonction des budgets publicitaires à la télévision, à la radio et dans les journaux. Le travail consistait à tester chaque média, à ne garder que ceux qui ont un effet, puis à conseiller le client sur le média où investir, avec un intervalle de prévision. Sur les données laitières, step en mode descendant servait à simplifier le modèle complet, et l’AIC à départager les candidats.

À retenir

En régression multiple, $\hat\beta = (X^\top X)^{-1}X^\top Y$, et chaque coefficient est un effet à autres variables fixées. Student teste chaque variable en présence des autres ; Fisher teste le modèle entier. Pour choisir les variables, on ne regarde pas le $R^2$ (qui augmente toujours) mais le $R^2$ ajusté ou l’AIC, par exemple avec step.

En entretien

« Le R² augmente mécaniquement avec le nombre de variables, donc je compare les modèles avec un critère pénalisé comme l’AIC, ou mieux, sur des données de validation. Et j’interprète chaque coefficient toutes choses égales par ailleurs, en me méfiant de la colinéarité. »

En résumé

• Modèle $Y = X\beta + \varepsilon$, $X$ de taille $n \times (p+1)$ avec une colonne de 1 ; $\hat\beta = (X^\top X)^{-1}X^\top Y$, $\hat\sigma^2 = SCR/(n - p - 1)$.
• Un coefficient mesure un effet partiel ; une variable peut sembler inutile seule et être essentielle avec les autres, ou l’inverse.
• Test de Student par coefficient ($\mathcal{T}_{n-p-1}$), test de Fisher global ($\mathcal{F}_{p, n-p-1}$), test de Fisher partiel pour des modèles emboîtés.
• Le $R^2$ augmente toujours avec $p$ : on compare avec le $R^2$ ajusté ou l’AIC (plus petit = meilleur).
• Sélection exhaustive, descendante, ascendante ou pas à pas (step en R).

Exercices

Exercice · Construire la matrice X

facile

On modélise la consommation $y$ d’une voiture (L/100 km) par sa masse $x_1$ (tonnes) et sa puissance $x_2$ (ch) sur 4 voitures : (1,2 ; 90 ; 6,1), (1,5 ; 120 ; 7,0), (1,1 ; 75 ; 5,4), (1,8 ; 150 ; 8,3). Écrivez $X$ et $Y$, et donnez la taille de $X^\top X$. Combien de degrés de liberté reste-t-il pour estimer $\sigma^2$ ?

Voir la solution

$$X = \begin{pmatrix} 1 & 1{,}2 & 90 \\ 1 & 1{,}5 & 120 \\ 1 & 1{,}1 & 75 \\ 1 & 1{,}8 & 150 \end{pmatrix}, \qquad Y = \begin{pmatrix} 6{,}1 \\ 7{,}0 \\ 5{,}4 \\ 8{,}3 \end{pmatrix}$$

$X$ est de taille $4 \times 3$, donc $X^\top X$ est une matrice carrée $3 \times 3$. Il reste $n - p - 1 = 4 - 2 - 1 = 1$ degré de liberté pour estimer $\sigma^2$ : c’est bien trop peu pour une inférence sérieuse. Il faudrait beaucoup plus de voitures.

En code

Python

import numpy as np

masse = np.array([1.2, 1.5, 1.1, 1.8])
puissance = np.array([90, 120, 75, 150])
Y = np.array([6.1, 7.0, 5.4, 8.3])

X = np.column_stack([np.ones(4), masse, puissance])   # colonne de 1, puis x1 et x2
print(X.shape, (X.T @ X).shape)                        # (4, 3) (3, 3)
n, p = X.shape[0], X.shape[1] - 1
print(n - p - 1)                                       # 1 degré de liberté résiduel
print(np.linalg.solve(X.T @ X, X.T @ Y))               # 3.23  -1.67  0.053 : masse négative, estimation très fragile

R

voitures <- data.frame(masse     = c(1.2, 1.5, 1.1, 1.8),
                       puissance = c(90, 120, 75, 150),
                       conso     = c(6.1, 7.0, 5.4, 8.3))
X <- model.matrix(conso ~ masse + puissance, data = voitures)   # colonne de 1 incluse
Y <- voitures$conso
X                                          # matrice 4 x 3
dim(X); dim(t(X) %*% X)                    # 4 3  puis  3 3
reg <- lm(conso ~ masse + puissance, data = voitures)
df.residual(reg)                           # 1 degré de liberté résiduel
coef(reg)                                  # 3.23  -1.67  0.053 : masse négative, estimation très fragile

Exercice · Choisir entre deux modèles

moyen

Sur $n = 50$ observations, avec $SCT = 1\,000$ : le modèle A (2 variables) a $SCR = 400$ ; le modèle B (5 variables, dont celles de A) a $SCR = 380$. Calculez le $R^2$ et le $R^2$ ajusté des deux modèles. Lequel préférer ?

Voir la solution

• Modèle A : $R^2 = 1 - 400/1\,000 = 0{,}60$ ; $R^2_{aj} = 1 - \frac{400/47}{1\,000/49} = 1 - \frac{8{,}511}{20{,}408} \approx 0{,}583$.
• Modèle B : $R^2 = 0{,}62$ ; $R^2_{aj} = 1 - \frac{380/44}{20{,}408} = 1 - \frac{8{,}636}{20{,}408} \approx 0{,}577$.

B a un meilleur $R^2$, comme toujours avec plus de variables, mais un moins bon $R^2$ ajusté : les trois variables supplémentaires réduisent trop peu la $SCR$ pour justifier leur coût. On préfère A. Un test de Fisher partiel le confirmerait : $F = \frac{(400 - 380)/3}{380/44} \approx 0{,}77$, bien en dessous du seuil de $\mathcal{F}_{3,44}$ (environ 2,8).

En code

Python

from scipy import stats

n, SCT = 50, 1000
modeles = {"A": (2, 400), "B": (5, 380)}          # (p, SCR)
for nom, (p, SCR) in modeles.items():
    R2 = 1 - SCR / SCT
    R2_aj = 1 - (SCR / (n - p - 1)) / (SCT / (n - 1))
    print(nom, R2, round(R2_aj, 3))                # A 0.6 0.583 ; B 0.62 0.577

# Test de Fisher partiel : les 3 variables de B en plus de A
F = ((400 - 380) / 3) / (380 / 44)
print(F, stats.f.ppf(0.95, 3, 44), stats.f.sf(F, 3, 44))   # 0.77  2.82  p = 0.52

R

n <- 50; SCT <- 1000
p   <- c(A = 2, B = 5)
SCR <- c(A = 400, B = 380)
R2    <- 1 - SCR / SCT
R2_aj <- 1 - (SCR / (n - p - 1)) / (SCT / (n - 1))
R2; round(R2_aj, 3)                         # A 0.60 0.583 ; B 0.62 0.577

# Test de Fisher partiel : les 3 variables de B en plus de A
F <- ((400 - 380) / 3) / (380 / 44)
F; qf(0.95, 3, 44); 1 - pf(F, 3, 44)        # 0.77  2.82  p = 0.52

Exercice · Un modèle non linéaire rendu linéaire

difficile

Une loi physique relie la pression $P$ à la température $T$ par $P = \alpha_1 T^{\alpha_2} e^{\alpha_3/T}$. Montrez qu’on peut estimer $\alpha_1, \alpha_2, \alpha_3$ par une régression linéaire multiple, en précisant la réponse et les variables explicatives. Quelle précaution prendre si une température vaut 0 ?

Voir la solution

En passant au logarithme : $\ln P = \ln\alpha_1 + \alpha_2 \ln T + \alpha_3 \cdot \frac{1}{T}$. C’est un modèle linéaire en les paramètres :

$$\ln P_i = \beta_0 + \beta_1 \ln T_i + \beta_2 \frac{1}{T_i} + \varepsilon_i$$

avec la réponse $\ln P$, les deux variables explicatives $\ln T$ et $1/T$, et $\beta_0 = \ln\alpha_1$, $\beta_1 = \alpha_2$, $\beta_2 = \alpha_3$. On estime par lm(log(P) ~ log(T) + I(1/T)), puis $\hat\alpha_1 = e^{\hat\beta_0}$.

Si une température vaut 0, ni $\ln T$ ni $1/T$ ne sont définis : il faut soit retirer ce point, soit (si c’est justifié physiquement) travailler en kelvins, où la température est strictement positive. C’est exactement la situation d’un exercice de TP sur la tension de vapeur du mercure, où l’on contournait le problème en remplaçant la température nulle par une petite valeur positive.

En code

Python

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

# Données simulées : P = 2 * T^0.8 * exp(-20 / T), bruit multiplicatif
rng = np.random.default_rng(0)
temp = np.linspace(5, 100, 25)                     # températures, toutes > 0
P = 2 * temp**0.8 * np.exp(-20 / temp) * np.exp(rng.normal(0, 0.02, size=25))
donnees = pd.DataFrame({"temp": temp, "P": P})

# Réponse log(P), variables explicatives log(T) et 1/T
reg = smf.ols("np.log(P) ~ np.log(temp) + I(1 / temp)", data=donnees).fit()
b0, b1, b2 = reg.params
print(np.exp(b0), b1, b2)          # 1.96  0.804  -19.90 (vrais : 2, 0.8, -20)

R

# Données simulées : P = 2 * T^0.8 * exp(-20 / T), bruit multiplicatif
set.seed(0)
temp <- seq(5, 100, length.out = 25)        # températures, toutes > 0
P <- 2 * temp^0.8 * exp(-20 / temp) * exp(rnorm(25, 0, 0.02))

# Réponse log(P), variables explicatives log(T) et 1/T
reg <- lm(log(P) ~ log(temp) + I(1 / temp))
b <- coef(reg)
c(alpha1 = exp(b[[1]]), alpha2 = b[[2]], alpha3 = b[[3]])
## alpha1 1.98, alpha2 0.801, alpha3 -19.84 (vrais : 2, 0.8, -20)