M2SP — Régression linéaire. Ce TD suit pas à pas le cours LIN_REG (diapositives 1 à 105). Chaque section rappelle les notions statistiques, puis propose le code R qui les met en œuvre, avec une explication du code. Les notions s’appuient aussi sur Cornillon et al., Régression avec R ; James et al., An Introduction to Statistical Learning (ISLR, chap. 3) ; Lilja, Linear Regression Using R.
Durée indicative : 3 h. Les sections marquées ★ sont facultatives (à traiter si le temps le permet, ou en autonomie).
| Partie | Contenu | Diapos | Durée |
|---|---|---|---|
| 0 | Mise en route, données | — | 10 min |
| 1 | Régression linéaire simple | 3–43 | 70 min |
| 2 | Régression multiple, propriétés des estimateurs, EQM, moments, MV | 44–67 | 35 min |
| 3 | Validation du modèle | 68–91 | 40 min |
| 4 | Variables qualitatives : ANOVA, ANCOVA | 92–105 | 20 min |
| 5 | Synthèse : modèles complets et cas à diagnostiquer | tout | 25 min |
Philosophie du TD. Pour chaque notion, on calcule
d’abord à la main (avec des vecteurs et des matrices), puis on
vérifie avec les fonctions de R (lm, summary,
predict, …). Quand une propriété est « en moyenne sur les
échantillons » (biais, variance, couverture d’un IC), on la vérifie
par simulation : c’est le moyen le plus sûr de
comprendre ce que dit un théorème.
Tout le TD fonctionne en R de base. Deux packages
sont utilisés ponctuellement pour vérifier des calculs faits à
la main : car (VIF) et lmtest (test de
Durbin-Watson).
# install.packages(c("car", "lmtest")) # une seule fois
set.seed(2026) # graine : rend les simulations reproductiblesLes quatre fichiers proviennent de l’association Air Breizh (Rennes)
et sont ceux utilisés dans Cornillon et al. Placez-les dans un
dossier data/ à côté de votre script.
| Fichier | n | Contenu |
|---|---|---|
ozone_simple.txt |
50 | O3 (maximum journalier d’ozone, µg/m³) et
T12 (température à 12 h) |
ozone.txt |
50 | mêmes jours, avec T15, nébulosité Ne12,
vents N12 S12 E12 W12, Vx, ozone de la veille
O3v, et deux facteurs nebulosite,
vent |
ozone_long.txt |
1014 | O3, T6, T12,
Ne12, Ne15, Vx,
O3v |
ozone_complet.txt |
1464 | toutes les mesures (températures, nébulosité, direction et vitesse du vent à 6/9/12/15/18 h), 1995–2002, avec valeurs manquantes |
dossier <- "data"
ozone_simple <- read.table(file.path(dossier, "ozone_simple.txt"), header = TRUE, sep = ";")
ozone <- read.table(file.path(dossier, "ozone.txt"), header = TRUE, sep = ";",
stringsAsFactors = TRUE)
ozone_long <- read.table(file.path(dossier, "ozone_long.txt"), header = TRUE, sep = ";")
ozone_complet <- read.table(file.path(dossier, "ozone_complet.txt"), header = TRUE, sep = ";")
str(ozone_simple)
str(ozone)
dim(ozone_long)
dim(ozone_complet)
head(rownames(ozone_complet))Explication du code. -
read.table(..., header = TRUE, sep = ";") lit un fichier
texte dont la première ligne contient les noms de colonnes et dont le
séparateur est le point-virgule. - stringsAsFactors = TRUE
transforme les colonnes de texte (nebulosite,
vent) en facteurs : c’est ce type qui
permet à R de les traiter comme variables qualitatives (partie 4). -
Dans ozone_complet.txt, la ligne d’en-tête contient une
colonne de moins que les lignes de données : R utilise alors
automatiquement la première colonne (la date AAAAMMJJ)
comme noms de lignes. On la récupère comme vraie date
:
ozone_complet$date <- as.Date(rownames(ozone_complet), format = "%Y%m%d")
colSums(is.na(ozone_complet)) # valeurs manquantes par variable
oc <- na.omit(ozone_complet) # jeu sans valeur manquante, utilisé en partie 3 et 5
nrow(oc)Notions. On cherche à expliquer le maximum journalier d’ozone \(y_i\) par la température à midi \(x_i\), avec deux objectifs : ajuster un modèle (comprendre la liaison) et prédire \(y\) pour de nouvelles valeurs de \(x\). Toute étude de régression commence par le nuage de points \((x_i, y_i)\).
x <- ozone_simple$T12
y <- ozone_simple$O3
n <- length(y)
head(ozone_simple, 5) # les 5 premières mesures de la diapo 4
summary(ozone_simple)
plot(x, y, pch = 19, col = "steelblue",
xlab = "T12 (°C)", ylab = "max O3 (µg/m³)", main = "Ozone en fonction de la température")
cor(x, y)Explication. On stocke les deux variables dans des
vecteurs courts x et y pour alléger les
formules qui suivent. cor() donne le coefficient de
corrélation linéaire de Pearson \(r\) ;
on verra en 1.11 que \(R^2 = r^2\) en
régression simple.
Le mot régression vient de Galton : les fils de pères grands sont en moyenne plus petits que leur père (et inversement). Ce n’est pas un phénomène biologique mais une conséquence d’une corrélation \(|\rho| < 1\) :
pere <- rnorm(1000, mean = 175, sd = 7)
fils <- 175 + 0.5 * (pere - 175) + rnorm(1000, 0, 7 * sqrt(1 - 0.5^2)) # corrélation 0.5
plot(pere, fils, pch = 20, col = "grey50", asp = 1)
abline(0, 1, lty = 2) # « le fils a la taille du père »
abline(lm(fils ~ pere), col = "red", lwd = 2) # droite de régression : pente 0.5 < 1
tapply(fils - pere, pere > 175, mean) # grands pères : fils plus petits qu'euxLa simulation construit fils avec une pente 0.5 et un
bruit choisi pour que fils ait la même variance que
pere. La droite de régression (rouge) est « aplatie » par
rapport à la bissectrice : c’est le taller than mediocrity de
Galton.
Notions. Donner un sens à \(y_i \approx f(x_i)\) demande de choisir une fonction de coût \(l\) et une classe de fonctions \(\mathcal{G}\), puis de résoudre
\hat f = \operatorname*{argmin}_{f \in \mathcal{G}} \sum_{i=1}^n l\big(y_i - f(x_i)\big).Deux choix classiques : le coût quadratique \(l(u)=u^2\) (moindres carrés) et le coût absolu \(l(u)=|u|\) (moindres écarts absolus). Le coût quadratique pénalise beaucoup plus les grands écarts : il est plus simple à calculer (solution explicite) mais moins robuste aux points éloignés.
# Les deux fonctions de coût
curve(x^2, from = -3, to = 3, lwd = 2, xlab = "u", ylab = "l(u)")
curve(abs(x), add = TRUE, lty = 2, lwd = 2)
legend("top", c("quadratique", "absolu"), lty = 1:2, lwd = 2, bty = "n")
# Critères à minimiser pour f(x) = b1 + b2 x
S_quad <- function(b, x, y) sum((y - b[1] - b[2] * x)^2)
S_abs <- function(b, x, y) sum(abs(y - b[1] - b[2] * x))
b_quad <- optim(c(mean(y), 0), S_quad, x = x, y = y, method = "BFGS")$par
b_abs <- optim(b_quad, S_abs, x = x, y = y, control = list(maxit = 5000))$par
rbind(quadratique = b_quad, absolu = b_abs)Explication. optim() minimise
numériquement une fonction : on lui donne un point de départ
(c(mean(y), 0) : droite horizontale) et la fonction à
minimiser ; les arguments supplémentaires (x = x, y = y)
sont transmis à la fonction. Pour le coût absolu, non dérivable, on
garde la méthode par défaut (Nelder-Mead), qui n’utilise pas de
gradient.
Robustesse. Ajoutons un point très éloigné (journée chaude avec peu d’ozone) :
x_pert <- c(x, 32)
y_pert <- c(y, 20)
b_quad_p <- optim(c(mean(y_pert), 0), S_quad, x = x_pert, y = y_pert, method = "BFGS")$par
b_abs_p <- optim(b_quad_p, S_abs, x = x_pert, y = y_pert, control = list(maxit = 5000))$par
plot(x_pert, y_pert, pch = 19, col = c(rep("grey40", n), "red"), xlab = "T12", ylab = "O3")
abline(b_quad, col = "blue"); abline(b_quad_p, col = "blue", lty = 2)
abline(b_abs, col = "darkgreen"); abline(b_abs_p, col = "darkgreen", lty = 2)
legend("topleft", c("quadratique", "quadratique + point", "absolu", "absolu + point"),
col = c("blue", "blue", "darkgreen", "darkgreen"), lty = c(1, 2, 1, 2), bty = "n")
rbind(quad_sans = b_quad, quad_avec = b_quad_p, abs_sans = b_abs, abs_avec = b_abs_p)À vous. Quelle pente bouge le plus ? Relier à la phrase de la diapo 12 : coût quadratique peu robuste. Gauss choisit néanmoins les moindres carrés pour leur simplicité de calcul (et on verra en 2.7 qu’ils ont une justification probabiliste).
Notions. Si \(\mathcal{G}\) est trop large, une infinité de fonctions interpolent les données et annulent le critère : le critère seul ne suffit pas. On restreint \(\mathcal{G}\) par connaissance a priori ou par l’observation du nuage, par exemple aux droites : \(\mathcal{G}=\{f : f(x)=ax+b\}\).
set.seed(1)
petit <- data.frame(x = 1:6)
petit$y <- 2 + 0.5 * petit$x + rnorm(6, sd = 0.5)
interp <- lm(y ~ poly(x, 5, raw = TRUE), data = petit) # polynôme de degré 5 : passe par les 6 points
droite <- lm(y ~ x, data = petit)
c(SCR_interpolation = sum(resid(interp)^2), SCR_droite = sum(resid(droite)^2))
grille <- data.frame(x = seq(0.8, 6.2, length.out = 300))
plot(petit, pch = 19, ylim = range(c(petit$y, predict(interp, grille))))
lines(grille$x, predict(interp, grille), col = "red", lwd = 2)
abline(droite, col = "blue", lwd = 2)Explication. poly(x, 5, raw = TRUE)
crée les colonnes \(x, x^2, \dots,
x^5\) : avec 6 coefficients et 6 points, le polynôme passe
exactement par tous les points (SCR ≈ 0), mais il oscille entre eux et
serait catastrophique en prédiction.
Le choix de \(\mathcal{G}\) se fait en regardant le nuage (diapo 15, graphiques a, b, c) :
xs <- runif(60, 0, 10)
par(mfrow = c(1, 3))
plot(xs, 3 + rnorm(60), main = "(a) pas de liaison")
plot(xs, (xs - 5)^2 + rnorm(60, 0, 2), main = "(b) liaison non linéaire")
plot(xs, 2 + 1.5 * xs + rnorm(60, 0, 2), main = "(c) liaison linéaire")
par(mfrow = c(1, 1))Notions. Le modèle de régression linéaire simple s’écrit
Y_i = \beta_1 + \beta_2 x_i + \varepsilon_i, \qquad i = 1,\dots,n,où \(\beta_1, \beta_2\) sont fixes mais inconnus, les \(x_i\) sont connus (non aléatoires) et \(\varepsilon_i\) est une erreur aléatoire inconnue. Hypothèses :
H1 ⇒ existence/unicité de \(\hat\beta\) ; H1+H2 ⇒ Gauss–Markov ; H1+H3 ⇒ lois exactes (tests, IC).
Pour s’entraîner, on écrit une fonction qui simule des données sous le modèle. C’est l’outil central de tout le TD : comme on connaît les vrais paramètres, on peut vérifier ce que font les estimateurs.
simule_reg <- function(x, beta1 = 10, beta2 = 2, sigma = 3) {
data.frame(x = x, y = beta1 + beta2 * x + rnorm(length(x), mean = 0, sd = sigma))
}
x_plan <- seq(1, 10, length.out = 10) # plan d'expérience fixe (n = 10)
d1 <- simule_reg(x_plan)
coef(lm(y ~ x, data = d1)) # proche de (10, 2), mais pas égal
# Violation de H1 : tous les x identiques
d_H1 <- data.frame(x = rep(5, 10), y = rnorm(10))
coef(lm(y ~ x, data = d_H1)) # pente = NA : non identifiableExplication. simule_reg() prend un
vecteur de \(x\) fixés et fabrique
\(y\) selon le modèle, avec des erreurs
gaussiennes (H3, donc H2). Lorsque tous les \(x_i\) sont égaux, la colonne \(x\) est proportionnelle à la constante : R
ne peut pas séparer \(\beta_1\) de
\(\beta_2\) et renvoie
NA.
Notions. Les estimateurs MC minimisent
S(\beta_1,\beta_2) = \sum_{i=1}^n (y_i - \beta_1 - \beta_2 x_i)^2 = \|Y - \beta_1\mathbf{1} - \beta_2 X\|^2 .\(S\) est strictement convexe (sous H1). En annulant les dérivées partielles on obtient les équations normales
\sum_i (y_i - \hat\beta_1 - \hat\beta_2 x_i) = 0, \qquad \sum_i x_i\,(y_i - \hat\beta_1 - \hat\beta_2 x_i) = 0,d’où
\hat\beta_2 = \frac{\sum_i (x_i-\bar x)\,y_i}{\sum_i (x_i-\bar x)^2}, \qquad \hat\beta_1 = \bar y - \hat\beta_2\,\bar x .La droite passe par le centre de gravité \((\bar x, \bar y)\).
xb <- mean(x)
yb <- mean(y)
sxx <- sum((x - xb)^2)
b2_chap <- sum((x - xb) * y) / sxx
b1_chap <- yb - b2_chap * xb
c(b1_chap, b2_chap)
reg <- lm(O3 ~ T12, data = ozone_simple)
coef(reg) # identique
# Les équations normales sont vérifiées
e_chap <- y - b1_chap - b2_chap * x
c(somme_residus = sum(e_chap), somme_x_residus = sum(x * e_chap)) # ≈ 0 (erreurs d'arrondi)
# Centre de gravité
plot(x, y, pch = 19, col = "grey50", xlab = "T12", ylab = "O3")
abline(reg, col = "blue", lwd = 2)
points(xb, yb, pch = 4, cex = 3, lwd = 3, col = "red")
predict(reg, newdata = data.frame(T12 = xb)); ybExplication. On code directement les formules du
cours. lm(O3 ~ T12, data = ozone_simple) ajuste le même
modèle ; la formule O3 ~ T12 inclut la constante par
défaut. predict() au point \(\bar
x\) renvoie bien \(\bar y\).
g1 <- seq(b1_chap - 40, b1_chap + 40, length.out = 100)
g2 <- seq(b2_chap - 2, b2_chap + 2, length.out = 100)
S_grille <- outer(g1, g2, Vectorize(function(u, v) S_quad(c(u, v), x, y)))
contour(g1, g2, S_grille, nlevels = 25, xlab = "beta1", ylab = "beta2")
points(b1_chap, b2_chap, pch = 19, col = "red")outer() évalue \(S\)
sur toute une grille de valeurs \((\beta_1,\beta_2)\) ; les courbes de niveau
sont des ellipses emboîtées autour d’un minimum unique. Leur forme
allongée et inclinée annonce la corrélation entre \(\hat\beta_1\) et \(\hat\beta_2\) (1.6).
Notions. \(\hat\beta_1,\hat\beta_2\) sont des variables aléatoires : ils changent d’un échantillon à l’autre, alors que \(\beta_1, \beta_2\) restent fixes. Sous H1–H2 :
E(\hat\beta_j)=\beta_j,\qquad
V(\hat\beta_2)=\frac{\sigma^2}{\sum (x_i-\bar x)^2},\qquad
V(\hat\beta_1)=\frac{\sigma^2\sum x_i^2}{n\sum (x_i-\bar x)^2},\qquad
\operatorname{Cov}(\hat\beta_1,\hat\beta_2)=-\frac{\sigma^2\,\bar x}{\sum (x_i-\bar x)^2}.Les variances sont faibles si \(\sigma^2\) est petit, si les \(x_i\) sont dispersés, et (pour \(\hat\beta_1\)) si les \(x_i\) sont proches de 0. Gauss–Markov : parmi les estimateurs linéaires en \(Y\) et sans biais, les MC sont de variance minimale.
# Trois échantillons de taille 10 (diapo 24)
plot(NULL, xlim = c(0, 11), ylim = c(0, 40), xlab = "x", ylab = "y")
for (k in 1:3) {
d <- simule_reg(x_plan)
points(d, pch = 19, col = k + 1)
abline(lm(y ~ x, data = d), col = k + 1)
}
abline(10, 2, lwd = 3, lty = 2) # la vraie droite, inconnue en pratiqueMonte-Carlo. On répète l’expérience B = 5000 fois et on compare moyennes/variances empiriques aux formules :
B <- 5000
sigma <- 3
est <- t(replicate(B, coef(lm(y ~ x, data = simule_reg(x_plan, sigma = sigma)))))
colnames(est) <- c("b1", "b2")
sxx_plan <- sum((x_plan - mean(x_plan))^2)
theorie <- c(V_b1 = sigma^2 * sum(x_plan^2) / (10 * sxx_plan),
V_b2 = sigma^2 / sxx_plan,
Cov = -sigma^2 * mean(x_plan) / sxx_plan)
empirique <- c(var(est[, 1]), var(est[, 2]), cov(est[, 1], est[, 2]))
colMeans(est) # ≈ (10, 2) : sans biais
round(rbind(theorie, empirique), 4)
hist(est[, 2], breaks = 50, freq = FALSE, main = "Distribution de la pente estimée", xlab = "b2")
abline(v = 2, col = "red", lwd = 2)Explication. replicate(B, expr) évalue
expr B fois ; chaque évaluation simule un nouvel
échantillon et renvoie les deux coefficients. t() met le
résultat sous forme d’une matrice B × 2. Les moyennes empiriques
approchent les vraies valeurs (absence de biais) et les variances
empiriques les formules théoriques.
Influence du plan d’expérience.
var_pente <- function(xp) var(replicate(2000, coef(lm(y ~ x, data = simule_reg(xp)))[2]))
var_const <- function(xp) var(replicate(2000, coef(lm(y ~ x, data = simule_reg(xp)))[1]))
plans <- list(resserre = seq(5, 6, length.out = 10),
disperse = seq(1, 10, length.out = 10),
loin_de_0 = seq(101, 110, length.out = 10))
sapply(plans, function(p) c(V_b2 = var_pente(p), V_b1 = var_const(p)))À vous. Interpréter les trois colonnes à l’aide de la diapo 29.
Gauss–Markov illustré. La « pente des extrêmes » \((y_{10}-y_1)/(x_{10}-x_1)\) est aussi linéaire en \(Y\) et sans biais, mais de variance \(2\sigma^2/(x_{10}-x_1)^2\) plus grande :
pente_ext <- replicate(B, { d <- simule_reg(x_plan); (d$y[10] - d$y[1]) / (d$x[10] - d$x[1]) })
rbind(moyenne = c(MC = mean(est[, 2]), extremes = mean(pente_ext)),
variance = c(MC = var(est[, 2]), extremes = var(pente_ext)))Notions. Les résidus \(\hat\varepsilon_i = y_i - \hat y_i\) estiment les erreurs inconnues. On estime \(\sigma^2\) par
\hat\sigma^2 = \frac{1}{n-p}\sum_{i=1}^n \hat\varepsilon_i^2 = \frac{SCR}{n-p}, \qquad p = 2 \text{ en régression simple.}Diviser par \(n-p\) (et non \(n\)) rend l’estimateur sans biais : les résidus sont liés par \(p\) contraintes (\(\sum\hat\varepsilon_i=0\), \(\sum x_i\hat\varepsilon_i=0\)), il ne reste que \(n-p\) degrés de liberté. \(\hat\sigma\) est le Residual standard error de R.
y_chap <- b1_chap + b2_chap * x
eps_chap <- y - y_chap
all.equal(unname(eps_chap), unname(resid(reg)))
SCR <- sum(eps_chap^2)
sigma2_chap <- SCR / (n - 2)
c(sigma_chap = sqrt(sigma2_chap), sigma_R = sigma(reg))
# Biais de SCR/n : simulation avec sigma^2 = 9
s2 <- t(replicate(B, {
r <- resid(lm(y ~ x, data = simule_reg(x_plan, sigma = 3)))
c(div_n = sum(r^2) / 10, div_n_moins_2 = sum(r^2) / 8)
}))
colMeans(s2) # la première sous-estime 9 d'un facteur (n-2)/n = 0.8Explication. resid() et
sigma() sont les extracteurs de R ; on vérifie qu’ils
coïncident avec nos calculs. La simulation montre que \(SCR/n\) vaut en moyenne \(\sigma^2 (n-2)/n\).
Notions. Espace des individus (\(\mathbb{R}^2\)) : la droite minimise la somme des carrés des distances verticales. Espace des variables (\(\mathbb{R}^n\)) : \(\hat Y\) est la projection orthogonale de \(Y\) sur le sous-espace engendré par \(\mathbf 1\) et \(X\) :
\hat Y = P_X Y, \qquad P_X = X(X'X)^{-1}X', \qquad P_X^2 = P_X,\quad P_X' = P_X,\quad \operatorname{tr}(P_X)=p .Le vecteur des résidus est orthogonal à \(\mathbf 1\) et à \(X\).
# Espace des individus : résidus = segments verticaux (diapo 34)
o <- order(x)
i9 <- o[9] # 9e plus petite valeur de x
plot(x, y, pch = 19, col = "grey50", xlab = "T12", ylab = "O3")
abline(reg, col = "blue", lwd = 2)
segments(x, y, x, fitted(reg), col = "grey70")
segments(x[i9], y[i9], x[i9], fitted(reg)[i9], col = "red", lwd = 3)
points(x[i9], fitted(reg)[i9], pch = 4, col = "red", cex = 2)
# Espace des variables : matrice de projection
X <- cbind(1, x) # matrice n x 2 : constante et T12
P <- X %*% solve(t(X) %*% X) %*% t(X)
Y_chap <- P %*% y
all.equal(c(Y_chap), unname(fitted(reg)))
all.equal(P %*% P, P) # idempotente
isSymmetric(P) # symétrique
sum(diag(P)) # trace = p = 2
t(X) %*% (y - Y_chap) # ≈ 0 : résidus orthogonaux à 1 et à XExplication. %*% est le produit
matriciel, t() la transposée, solve(A)
l’inverse de A. order(x) donne les indices qui trient
x : o[9] est l’individu \(x_{(9)}\) de la diapo 34.
Notions. Sous H1 + H3 : - \(\hat\beta_j \sim \mathcal{N}(\beta_j, \sigma^2_{\hat\beta_j})\) ; \((n-2)\hat\sigma^2/\sigma^2 \sim \chi^2_{n-2}\), indépendant de \(\hat\beta\) ; - \(\sigma\) inconnu : \(\dfrac{\hat\beta_j-\beta_j}{\hat\sigma_{\hat\beta_j}} \sim \mathcal{T}_{n-2}\), avec \(\hat\sigma_{\hat\beta_2} = \hat\sigma / \sqrt{\sum (x_i-\bar x)^2}\) ; - région de confiance jointe : \(\dfrac{1}{2\hat\sigma^2}(\hat\beta-\beta)'V^{-1}(\hat\beta-\beta) \sim \mathcal{F}_{2,n-2}\), avec \(V^{-1} = X'X\).
Test de liaison : \(H_0:\beta_2=0\) contre \(H_1:\beta_2\neq 0\), statistique \(T = \hat\beta_2/\hat\sigma_{\hat\beta_2}\) ; on rejette si \(|T| > t_{n-2}(1-\alpha/2)\). IC : \(\hat\beta_2 \pm t_{n-2}(1-\alpha/2)\,\hat\sigma_{\hat\beta_2}\).
se_b2 <- sqrt(sigma2_chap / sxx)
se_b1 <- sqrt(sigma2_chap * sum(x^2) / (n * sxx))
T_b2 <- b2_chap / se_b2
p_b2 <- 2 * pt(-abs(T_b2), df = n - 2)
t_q <- qt(0.975, df = n - 2)
IC_b2 <- b2_chap + c(-1, 1) * t_q * se_b2
IC_b1 <- b1_chap + c(-1, 1) * t_q * se_b1
c(se_b2 = se_b2, T = T_b2, p = p_b2, quantile = t_q)
summary(reg)$coefficients
rbind(IC_b1, IC_b2); confint(reg)Explication. pt() est la fonction de
répartition de Student, qt() sa fonction quantile. La
p-valeur bilatérale vaut \(2\,P(\mathcal{T}_{n-2} > |T|)\).
Vérification par simulation : loi de Student et couverture de l’IC.
sim_T <- t(replicate(B, {
m <- lm(y ~ x, data = simule_reg(x_plan))
cf <- summary(m)$coefficients
c(T = (cf[2, 1] - 2) / cf[2, 2], # (b2 - beta2) / se
couvre = abs(cf[2, 1] - 2) <= qt(0.975, 8) * cf[2, 2])
}))
hist(sim_T[, "T"], breaks = 60, freq = FALSE, xlim = c(-6, 6), main = "(b2 - beta2)/se", xlab = "")
curve(dt(x, df = 8), add = TRUE, col = "red", lwd = 2)
curve(dnorm(x), add = TRUE, col = "blue", lty = 2)
mean(sim_T[, "couvre"]) # ≈ 0.95Région de confiance jointe vs rectangle des deux IC.
XtX <- t(X) %*% X # V^{-1}
bhat <- c(b1_chap, b2_chap)
g1 <- seq(b1_chap - 4 * se_b1, b1_chap + 4 * se_b1, length.out = 150)
g2 <- seq(b2_chap - 4 * se_b2, b2_chap + 4 * se_b2, length.out = 150)
Fq <- outer(g1, g2, Vectorize(function(u, v) {
d <- bhat - c(u, v)
drop(t(d) %*% XtX %*% d) / (2 * sigma2_chap)
}))
contour(g1, g2, Fq, levels = qf(0.95, 2, n - 2), drawlabels = FALSE,
xlab = "beta1", ylab = "beta2", main = "Région de confiance à 95 %")
rect(IC_b1[1], IC_b2[1], IC_b1[2], IC_b2[2], border = "red", lty = 2)
points(b1_chap, b2_chap, pch = 19)L’ellipse est inclinée (covariance négative car \(\bar x>0\)) : certains couples dans le rectangle sont exclus de la région jointe, et inversement.
Significatif ≠ fort ≠ causal (diapo 40).
n_grand <- 1e5
xg <- rnorm(n_grand, 20, 5)
yg <- 50 + 0.05 * xg + rnorm(n_grand, 0, 20)
m_grand <- lm(yg ~ xg)
summary(m_grand)$coefficients["xg", ]
summary(m_grand)$r.squared # liaison dérisoire...
confint(m_grand)["xg", ] # ... mais « significative » : rapporter l'IC !Notions. Pour une nouvelle valeur \(x_{n+1}\), la prévision ponctuelle est \(\hat y_{n+1}=\hat\beta_1+\hat\beta_2x_{n+1}\).
\text{IC (moyenne)}:\ \hat y_{n+1} \pm t_{n-2}(1-\tfrac{\alpha}{2})\,\hat\sigma\sqrt{\tfrac1n + \tfrac{(x_{n+1}-\bar x)^2}{\sum (x_i-\bar x)^2}}
\qquad
\text{IP (individu)}:\ \hat y_{n+1} \pm t_{n-2}(1-\tfrac{\alpha}{2})\,\hat\sigma\sqrt{1+\tfrac1n + \tfrac{(x_{n+1}-\bar x)^2}{\sum (x_i-\bar x)^2}}x0 <- 25
y0 <- b1_chap + b2_chap * x0
se_moy <- sqrt(sigma2_chap * (1 / n + (x0 - xb)^2 / sxx))
se_pred <- sqrt(sigma2_chap * (1 + 1 / n + (x0 - xb)^2 / sxx))
rbind(IC = c(y0, y0 + c(-1, 1) * t_q * se_moy),
IP = c(y0, y0 + c(-1, 1) * t_q * se_pred))
nouveau <- data.frame(T12 = x0)
predict(reg, newdata = nouveau, interval = "confidence")
predict(reg, newdata = nouveau, interval = "prediction")
# Forme en sablier
grille <- data.frame(T12 = seq(min(x) - 8, max(x) + 8, length.out = 200))
ic <- predict(reg, grille, interval = "confidence")
ip <- predict(reg, grille, interval = "prediction")
plot(x, y, pch = 19, col = "grey50", xlim = range(grille$T12), ylim = range(ip),
xlab = "T12", ylab = "O3")
matlines(grille$T12, ic, lty = c(1, 2, 2), col = "blue")
matlines(grille$T12, ip[, 2:3], lty = 3, col = "red", lwd = 2)
abline(v = range(x), col = "grey", lty = 4) # domaine observé : ne pas extrapoler au-delàExplication. Dans predict(),
newdata doit être un data.frame dont les colonnes portent
les mêmes noms que dans la formule (T12).
matlines() trace plusieurs colonnes d’une matrice d’un
coup.
★ L’IP ne se resserre pas quand n augmente.
largeurs <- sapply(c(10, 100, 10000), function(nn) {
d <- simule_reg(runif(nn, 1, 10))
pr <- predict(lm(y ~ x, d), data.frame(x = 5), interval = "prediction")
pc <- predict(lm(y ~ x, d), data.frame(x = 5), interval = "confidence")
c(n = nn, largeur_IC = pc[3] - pc[2], largeur_IP = pr[3] - pr[2])
})
round(largeurs, 2) # IC -> 0 ; IP -> 2 x 1.96 x sigma ≈ 11.8Notions. Pythagore dans \(\mathbb{R}^n\) : \(SCT = SCE + SCR\), d’où
R^2 = \frac{SCE}{SCT} = 1 - \frac{SCR}{SCT}, \qquad
R^2_{aj} = 1 - \frac{SCR/(n-p)}{SCT/(n-1)} .En régression simple \(R^2 = r^2\). Le \(R^2\) croît mécaniquement avec le nombre de variables.
SCT <- sum((y - mean(y))^2)
SCE <- sum((fitted(reg) - mean(y))^2)
c(SCT = SCT, SCE_plus_SCR = SCE + SCR)
R2 <- SCE / SCT
R2aj <- 1 - (SCR / (n - 2)) / (SCT / (n - 1))
c(R2 = R2, r2 = cor(x, y)^2, R2_R = summary(reg)$r.squared,
R2aj = R2aj, R2aj_R = summary(reg)$adj.r.squared)
# Ajout de variables de pur bruit
bruit <- matrix(rnorm(n * 20), n, 20)
r2 <- sapply(0:20, function(k) {
m <- if (k == 0) reg else lm(y ~ x + bruit[, 1:k, drop = FALSE])
c(R2 = summary(m)$r.squared, R2aj = summary(m)$adj.r.squared)
})
matplot(0:20, t(r2), type = "b", pch = 19, lty = 1, col = c("black", "red"),
xlab = "nombre de variables de bruit ajoutées", ylab = "")
legend("topleft", c("R2", "R2 ajusté"), col = c("black", "red"), lty = 1, bty = "n")summary() (diapo
43)summary(reg)
anova(reg)
c(T_carre = summary(reg)$coefficients["T12", "t value"]^2,
F = summary(reg)$fstatistic["value"])Avec ozone_simple (n = 50) : \(\hat\beta_1 = 31.41\), \(\hat\beta_2 = 2.70\) (erreur-type 0.63),
\(\hat\sigma = 20.5\) sur 48 ddl, \(R^2 = 0.279\), \(F = 18.58 = 4.31^2\).
| Élément de la sortie | Où l’avons-nous calculé ? |
|---|---|
Estimate |
1.5 : b1_chap, b2_chap |
Std. Error |
1.9 : se_b1, se_b2 |
t value, Pr(>|t|) |
1.9 : T_b2, p_b2 |
Residual standard error |
1.7 : sqrt(sigma2_chap) |
Multiple / Adjusted R-squared |
1.11 : R2, R2aj |
F-statistic |
1.12 : \(F = T^2\) en régression simple |
Notions. On ajoute le vent Vx et la
nébulosité Ne12 :
y_i = \beta_1 + \beta_2 x_{i2} + \dots + \beta_p x_{ip} + \varepsilon_i, \qquad i=1,\dots,n,les \(x_{ij}\) sont connus, les \(\beta_j\) inconnus, les \(\varepsilon_i\) aléatoires. Le cours note \(\beta_1\) la constante (la colonne \(x_{i1}=1\)).
head(ozone[, c("O3", "T12", "Vx", "Ne12")], 5) # diapo 44
pairs(ozone[, c("O3", "T12", "Vx", "Ne12")], pch = 20)
round(cor(ozone[, c("O3", "T12", "Vx", "Ne12")]), 2)Explication. pairs() trace tous les
nuages de points deux à deux : on y voit les liaisons de O3
avec chaque variable, mais aussi les liaisons entre variables
explicatives (ici T12 et Ne12 sont
corrélées négativement), qui sont au cœur de l’interprétation « toutes
choses égales par ailleurs ».
Notions. \(Y_{n\times1} = X_{n\times p}\,\beta_{p\times1} + \varepsilon_{n\times1}\). Si \(\operatorname{rg}(X)=p\) (H1),
\hat\beta = \operatorname*{argmin}_{\beta\in\mathbb R^p}\|Y-X\beta\|^2 = (X'X)^{-1}X'Y .Y <- ozone$O3
X <- model.matrix(~ T12 + Vx + Ne12, data = ozone) # ajoute la colonne de 1
head(X)
dim(X); qr(X)$rank # rang = p = 4 : H1 vérifiée
beta_chap <- solve(t(X) %*% X) %*% t(X) %*% Y
reg_mult <- lm(O3 ~ T12 + Vx + Ne12, data = ozone)
cbind(a_la_main = drop(beta_chap), lm = coef(reg_mult))
# Numériquement préférable : résoudre le système X'X b = X'Y sans inverser
solve(crossprod(X), crossprod(X, Y))Explication. model.matrix() construit
la matrice du plan d’expérience à partir d’une formule, exactement comme
le fait lm() en interne. crossprod(X) calcule
\(X'X\) et
crossprod(X, Y) calcule \(X'Y\). En pratique lm()
n’inverse pas \(X'X\) mais utilise
une décomposition QR, plus stable quand les colonnes sont presque
colinéaires (3.8).
Le modèle est un hyperplan (diapos 49, 51). Simulons \(Y = 3X_1 + 4X_2 + \varepsilon\) :
n_s <- 100
x1 <- runif(n_s, 0, 10)
x2 <- runif(n_s, 0, 10)
y_plan <- 3 * x1 + 4 * x2 + rnorm(n_s, 0, 2)
m_plan <- lm(y_plan ~ x1 + x2)
coef(m_plan)
g <- seq(0, 10, length.out = 20)
z <- outer(g, g, function(a, b) predict(m_plan, data.frame(x1 = a, x2 = b)))
persp(g, g, z, theta = 30, phi = 20, col = "lightblue",
xlab = "X1", ylab = "X2", zlab = "Y", ticktype = "detailed")La surface est parfaitement plane : l’effet de \(X_1\) (+3 par unité) est le même quelle que soit la valeur de \(X_2\).
Notions. \(\hat\beta_j\) est la variation moyenne de \(Y\) associée à +1 unité de \(X_j\), les autres variables du modèle étant fixées. L’effet est conditionnel au modèle (il change si l’on ajoute/retire une variable), dépend de l’unité de \(X_j\), et n’est pas en soi causal.
coef(lm(O3 ~ T12, data = ozone))["T12"] # régression simple
coef(lm(O3 ~ T12 + Ne12, data = ozone))["T12"] # ajustée sur la nébulosité
cor(ozone$T12, ozone$Ne12)La pente de T12 passe d’environ 2.7 à environ 1.0 : dans
la régression simple, T12 « récupère » une partie de
l’effet de la nébulosité (jours chauds = jours ensoleillés).
# Changement d'unité : °C -> °F
T12F <- ozone$T12 * 9 / 5 + 32
c(celsius = coef(lm(O3 ~ T12 + Ne12, data = ozone))["T12"],
fahrenheit = coef(lm(O3 ~ T12F + Ne12, data = ozone))["T12F"],
ratio = 1.8)
# Confusion simulée : Z agit sur X et sur Y, X n'a AUCUN effet propre
z_c <- rnorm(200)
x_c <- z_c + rnorm(200, 0, 0.5)
y_c <- 3 * z_c + 0 * x_c + rnorm(200)
rbind(sans_Z = coef(lm(y_c ~ x_c))[2], avec_Z = coef(lm(y_c ~ x_c + z_c))[2])À vous. Dans la simulation, quelle est la vraie
valeur du coefficient de x_c ? Laquelle des deux
régressions l’estime ?
Notions. Le modèle \(y_i = \beta_1x_{i1} + \beta_2x_{i2} + \beta_3x_{i1}x_{i2} + \varepsilon_i\) permet à l’effet de \(X_1\) de dépendre de \(X_2\) : \(\partial E(Y)/\partial x_1 = \beta_1 + \beta_3 x_2\). Il n’existe plus « un » effet de \(X_1\).
y_int <- x1 + 3 * x2 + 6 * x1 * x2 + rnorm(n_s, 0, 5)
m_int <- lm(y_int ~ x1 * x2) # équivaut à x1 + x2 + x1:x2
coef(m_int)
cf <- coef(m_int)
sapply(c(x2_0 = 0, x2_5 = 5, x2_10 = 10), function(v) cf["x1"] + cf["x1:x2"] * v) # effet de X1
z2 <- outer(g, g, function(a, b) predict(m_int, data.frame(x1 = a, x2 = b)))
persp(g, g, z2, theta = 30, phi = 20, col = "salmon", xlab = "X1", ylab = "X2", zlab = "Y")
# Sur l'ozone : l'effet de la température dépend-il de la nébulosité ?
summary(lm(O3 ~ T12 * Ne12, data = ozone))$coefficientsExplication. Dans une formule R, a * b
développe en a + b + a:b, où a:b est le
produit. La surface n’est plus un plan mais une surface « vrillée ».
Notions. Sous H2 : \(E(\hat\beta)=(X'X)^{-1}X'X\beta=\beta\) (sans biais), et \(V(\hat\beta)=\sigma^2(X'X)^{-1}\), estimée par \(\hat\sigma^2(X'X)^{-1}\).
p_m <- ncol(X)
sigma2_m <- sum(resid(reg_mult)^2) / (nrow(X) - p_m)
V_chap <- sigma2_m * solve(crossprod(X))
all.equal(V_chap, vcov(reg_mult), check.attributes = FALSE)
sqrt(diag(V_chap)) # = colonne "Std. Error" de summary()
# Simulation avec le plan X de l'ozone et des paramètres connus
beta_vrai <- c(80, 1.5, 0.5, -5)
sim_b <- t(replicate(3000, {
Ys <- X %*% beta_vrai + rnorm(nrow(X), 0, 14)
drop(solve(crossprod(X), crossprod(X, Ys)))
}))
rbind(vrai = beta_vrai, moyenne_sim = colMeans(sim_b))
rbind(var_theorique = diag(14^2 * solve(crossprod(X))), var_sim = apply(sim_b, 2, var))Explication. vcov() renvoie la matrice
de variance-covariance estimée des coefficients. La simulation garde le
plan \(X\) fixe (hypothèse du cours :
\(x_{ij}\) non aléatoires) et ne tire
que les erreurs.
Notions. Pour un estimateur \(\hat\theta\) de \(\theta\) :
EQM(\hat\theta) = E\big[(\hat\theta-\theta)^2\big] = V(\hat\theta) + \big(E(\hat\theta)-\theta\big)^2 = \text{variance} + \text{biais}^2 .Un estimateur biaisé peut avoir une EQM plus faible qu’un estimateur sans biais si sa variance est plus petite. Parmi les estimateurs sans biais, celui de variance minimale est le MVUE ; pour un échantillon gaussien, \(\bar X\) est le MVUE de \(\mu\).
On compare trois estimateurs de \(\mu\) à partir de \(n=10\) observations \(\mathcal N(\mu=1,\sigma=2)\) : la moyenne \(\bar X\), la moyenne « rétrécie » \(0.7\,\bar X\), et la médiane.
mu <- 1; sig <- 2; n_e <- 10; c_r <- 0.7
ech <- replicate(10000, {
xs <- rnorm(n_e, mu, sig)
c(moyenne = mean(xs), retrecie = c_r * mean(xs), mediane = median(xs))
})
bilan <- apply(ech, 1, function(e) c(biais = mean(e) - mu, variance = var(e), EQM = mean((e - mu)^2)))
round(rbind(bilan, biais2_plus_var = bilan["biais", ]^2 + bilan["variance", ]), 3)
# Théorie : EQM(moyenne) = sig^2/n ; EQM(c * moyenne) = (1 - c)^2 mu^2 + c^2 sig^2 / n
curve(sig^2 / n_e + 0 * x, from = -3, to = 3, ylim = c(0, 1.5), lwd = 2,
xlab = "vraie valeur mu", ylab = "EQM")
curve((1 - c_r)^2 * x^2 + c_r^2 * sig^2 / n_e, add = TRUE, col = "red", lwd = 2)
legend("top", c("moyenne (sans biais)", "0.7 x moyenne (biaisée)"), col = c("black", "red"),
lwd = 2, bty = "n")Explication. apply(ech, 1, f) applique
f à chaque ligne de ech (un estimateur par
ligne). La décomposition biais² + variance retrouve l’EQM. L’estimateur
rétréci gagne quand \(\mu\) est proche
de 0, perd sinon : le biais peut payer, mais pas
partout. La médiane est sans biais mais de variance plus grande
que la moyenne (\(\approx \pi/2\) fois
pour \(n\) grand) : la moyenne est le
MVUE. C’est l’idée qui sous-tend la régression ridge (évoquée en
3.8).
Méthode des moments. On égale moments empiriques et théoriques : \(\frac1n\sum X_i = E(X)=\mu\), \(\frac1n\sum X_i^2 = E(X^2)=\sigma^2+\mu^2\).
xs <- rnorm(50, mean = 10, sd = 3)
m1 <- mean(xs); m2 <- mean(xs^2)
c(mu_moments = m1, sigma2_moments = m2 - m1^2, var_R = var(xs)) # var() divise par n-1
# Exemple où la méthode est vraiment utile : loi Gamma(forme k, taux lambda)
# E(X) = k / lambda, V(X) = k / lambda^2 => k = m^2 / v, lambda = m / v
x_gam <- rgamma(500, shape = 3, rate = 0.5)
m <- mean(x_gam); v <- mean(x_gam^2) - m^2
c(k_chap = m^2 / v, lambda_chap = m / v) # vrais : 3 et 0.5Maximum de vraisemblance, retour à la régression (diapo 67). Sous H3, \(Y\sim\mathcal N(X\beta,\sigma^2I_n)\) et
\ell(\beta,\sigma^2) = -\frac n2\ln(2\pi\sigma^2) - \frac{1}{2\sigma^2}\|Y-X\beta\|^2 .Maximiser \(\ell\) en \(\beta\) revient à minimiser \(\|Y-X\beta\|^2\) : \(\hat\beta_{MV}=\hat\beta_{MC}\). En revanche \(\hat\sigma^2_{MV}=SCR/n\) est biaisé.
moins_logvrais <- function(theta, X, Y) {
p <- ncol(X)
beta <- theta[1:p]
s2 <- exp(theta[p + 1]) # paramétrage log : garantit s2 > 0
-sum(dnorm(Y, mean = X %*% beta, sd = sqrt(s2), log = TRUE))
}
init <- c(mean(Y), rep(0, ncol(X) - 1), log(var(Y)))
mv <- optim(init, moins_logvrais, X = X, Y = Y, method = "BFGS",
control = list(maxit = 10000, reltol = 1e-14))
SCR_m <- sum(resid(reg_mult)^2)
cbind(MV = mv$par[1:4], MC = coef(reg_mult))
c(sigma2_MV = exp(mv$par[5]), SCR_sur_n = SCR_m / nrow(X), SCR_sur_n_moins_p = SCR_m / (nrow(X) - p_m))
c(logvrais_optim = -mv$value, logLik_R = as.numeric(logLik(reg_mult)))Explication. On écrit l’opposé de la
log-vraisemblance (car optim() minimise) avec
dnorm(..., log = TRUE). On optimise \(\ln\sigma^2\) plutôt que \(\sigma^2\) pour rester dans le domaine
autorisé. logLik() de R utilise précisément \(\hat\sigma^2_{MV}=SCR/n\) : c’est cette
log-vraisemblance qui sert au calcul de l’AIC et du BIC (3.9).
Notions (diapo 68). Une régression, c’est (1) une modélisation (\(Y=X\beta+\varepsilon\)), (2) une estimation, (3) une validation : vérifier H2/H3, la linéarité, et repérer les individus atypiques. L’examen des résidus est essentiellement graphique.
On travaille sur le modèle reg_mult :
O3 ~ T12 + Vx + Ne12, n = 50, p = 4. R propose un résumé
graphique immédiat, que l’on va reconstruire pièce par pièce :
par(mfrow = c(2, 2)); plot(reg_mult); par(mfrow = c(1, 1))Notions. Même sous H2, les résidus n’ont pas la même variance : \(V(\hat\varepsilon) = \sigma^2(I-P_X)\), soit \(V(\hat\varepsilon_i)=\sigma^2(1-h_{ii})\) avec \(h_{ii}\) le terme diagonal de \(P_X\). On les normalise :
t_i = \frac{\hat\varepsilon_i}{\hat\sigma\sqrt{1-h_{ii}}}\ \ \text{(standardisés)},\qquad
t_i^* = \frac{\hat\varepsilon_i}{\hat\sigma_{(i)}\sqrt{1-h_{ii}}}\ \ \text{(studentisés par validation croisée)},où \(\hat\sigma_{(i)}\) est estimé sans l’individu \(i\). Sous H1, H3, \(t_i^*\sim\mathcal T_{n-p-1}\).
X <- model.matrix(reg_mult)
n_m <- nrow(X); p_m <- ncol(X)
H <- X %*% solve(crossprod(X)) %*% t(X)
h <- diag(H)
all.equal(h, hatvalues(reg_mult), check.attributes = FALSE)
e <- resid(reg_mult)
s <- sigma(reg_mult)
t_std <- e / (s * sqrt(1 - h))
all.equal(t_std, rstandard(reg_mult))
# Studentisés par VC : on réajuste le modèle n fois, sans l'individu i
s_moins_i <- sapply(1:n_m, function(i) sigma(lm(O3 ~ T12 + Vx + Ne12, data = ozone[-i, ])))
t_etoile <- e / (s_moins_i * sqrt(1 - h))
all.equal(t_etoile, rstudent(reg_mult))
# Inutile de réajuster : formule explicite de sigma_(i)
s2_i <- ((n_m - p_m) * s^2 - e^2 / (1 - h)) / (n_m - p_m - 1)
all.equal(sqrt(s2_i), s_moins_i, check.attributes = FALSE)Explication. ozone[-i, ] retire la
ligne i. sapply(1:n_m, ...) boucle sur les
individus et collecte les \(\hat\sigma_{(i)}\). Les extracteurs
hatvalues(), rstandard() et
rstudent() donnent directement \(h_{ii}\), \(t_i\) et \(t_i^*\).
★ Les résidus bruts n’ont pas tous la même variance.
sim_e <- replicate(3000, resid(lm(X %*% beta_vrai + rnorm(n_m, 0, 14) ~ X - 1)))
plot(14^2 * (1 - h), apply(sim_e, 1, var), pch = 19,
xlab = "sigma^2 (1 - h_ii) théorique", ylab = "variance simulée du résidu i")
abline(0, 1, col = "red")Notions. Un individu est aberrant (atypique en \(Y\)) si \(|t_i^*| > t_{n-p-1}(1-\alpha/2)\).
seuil_t <- qt(0.975, df = n_m - p_m - 1)
plot(t_etoile, pch = 19, ylab = "t*_i", xlab = "indice", ylim = range(c(t_etoile, -3, 3)))
abline(h = c(-seuil_t, 0, seuil_t), lty = c(2, 1, 2), col = c("red", "grey", "red"))
which(abs(t_etoile) > seuil_t)
# Attention à la multiplicité : avec alpha = 5 %, on attend ~ 0.05 * n « aberrants » par hasard
0.05 * n_m
seuil_bonf <- qt(1 - 0.05 / (2 * n_m), df = n_m - p_m - 1) # correction de Bonferroni
which(abs(t_etoile) > seuil_bonf)
car::outlierTest(reg_mult) # même idée, package carÀ vous. Retrouver dans ozone les dates
des jours signalés (ozone$Date[...]). Sont-ils aberrants ou
simplement dans la queue attendue de la loi de Student ?
Notions. On compare les résidus ordonnés \(t^*_{(i)}\) aux quantiles théoriques. L’espérance de la \(i\)-ème statistique d’ordre d’une \(\mathcal N(0,1)\) est approchée par \(\Phi^{-1}\!\big(\tfrac{i-3/8}{n+1/4}\big)\) si \(n\le10\), \(\Phi^{-1}\!\big(\tfrac{i-1/2}{n}\big)\) sinon.
ppoints(5); (1:5 - 3/8) / (5 + 1/4) # ppoints() applique exactement la règle de la diapo
ppoints(50)[1:3]; ((1:3) - 1/2) / 50
q_th <- qnorm(ppoints(n_m))
plot(q_th, sort(t_etoile), pch = 19, xlab = "quantiles théoriques N(0,1)", ylab = "t* ordonnés")
abline(0, 1, col = "red")
qqnorm(t_etoile); qqline(t_etoile) # version R, identique
shapiro.test(t_etoile) # complément (peu puissant si n petit, trop si n grand)
# Contre-exemple : erreurs asymétriques
x_q <- runif(100, 0, 10)
y_q <- 1 + 2 * x_q + 3 * (rexp(100) - 1)
qqnorm(rstudent(lm(y_q ~ x_q)), main = "Erreurs exponentielles centrées"); qqline(rstudent(lm(y_q ~ x_q)))Explication. ppoints(n) renvoie les
probabilités \((i-a)/(n+1-2a)\) avec
\(a=3/8\) si \(n\le 10\) et \(a=1/2\) sinon : c’est la formule du cours.
Un Q-Q plot courbé en « banane » signale une asymétrie ; des extrémités
qui s’écartent en « S » signalent des queues lourdes.
Notions. Pas de test universel ; on trace \(t_i^*\) en fonction de \(\hat y_i\), et \(|t_i^*|\) avec un lisseur (lowess). Une structure (cône, tendance, vague) suggère une hétéroscédasticité.
par(mfrow = c(1, 2))
plot(fitted(reg_mult), t_etoile, pch = 19, xlab = "valeurs ajustées", ylab = "t*")
abline(h = 0, lty = 2)
plot(fitted(reg_mult), abs(t_etoile), pch = 19, xlab = "valeurs ajustées", ylab = "|t*|")
lines(lowess(fitted(reg_mult), abs(t_etoile)), col = "red", lwd = 2)
# Contre-exemple simulé : écart-type proportionnel à x
x_h <- runif(200, 1, 10)
y_h <- 5 + 2 * x_h + rnorm(200, 0, 0.5 * x_h)
m_h <- lm(y_h ~ x_h)
plot(fitted(m_h), rstudent(m_h), pch = 20, main = "cône"); abline(h = 0, lty = 2)
plot(fitted(m_h), abs(rstudent(m_h)), pch = 20)
lines(lowess(fitted(m_h), abs(rstudent(m_h))), col = "red", lwd = 2)
par(mfrow = c(1, 1))Variable oubliée (diapos 81–82). Si le vrai modèle est \(Y=\alpha+\beta_1X_1+\beta_2X_2+\varepsilon\) et qu’on ajuste \(Y=\alpha+\beta_1X_1+\varepsilon\), le terme \(\beta_2X_2\) passe dans les résidus.
n_o <- 100
x1o <- runif(n_o, 0, 10); x2o <- runif(n_o, 0, 10)
y_o <- 1 + 2 * x1o + 3 * x2o + rnorm(n_o)
m_oubli <- lm(y_o ~ x1o)
t_o <- rstudent(m_oubli)
par(mfrow = c(1, 3))
plot(t_o, pch = 19, xlab = "indice"); abline(h = 0, lty = 2)
plot(fitted(m_oubli), t_o, pch = 19, xlab = "Y ajusté"); abline(h = 0, lty = 2)
plot(x2o, t_o, pch = 19, xlab = "X2 (oubliée)"); abline(h = 0, lty = 2)
par(mfrow = c(1, 1))Seul le graphique contre \(X_2\) révèle la structure : il faut tracer les résidus contre toute variable candidate.
Structure temporelle (diapo 83). Les données
ozone_complet sont journalières : l’ozone d’un jour
ressemble à celui de la veille. Test de Durbin-Watson :
\(DW=\sum_{i\ge2}(\hat\varepsilon_i-\hat\varepsilon_{i-1})^2/\sum\hat\varepsilon_i^2\approx2(1-\hat\rho_1)\)
; \(DW\approx2\) sans autocorrélation,
\(DW<2\) si autocorrélation
positive.
m_temps <- lm(maxO3 ~ T12 + Ne12 + Vx, data = oc)
e_t <- resid(m_temps)
plot(oc$date, e_t, type = "l", xlab = "date", ylab = "résidus"); abline(h = 0, col = "red")
acf(e_t, main = "Autocorrélation des résidus")
DW <- sum(diff(e_t)^2) / sum(e_t^2)
rho1 <- cor(e_t[-1], e_t[-length(e_t)])
c(DW = DW, deux_fois_1_moins_rho = 2 * (1 - rho1))
lmtest::dwtest(m_temps)
# Ajouter l'ozone de la veille capte l'essentiel de la dépendance temporelle
m_temps2 <- update(m_temps, . ~ . + maxO3v)
lmtest::dwtest(m_temps2)Explication. diff(e_t) calcule \(\hat\varepsilon_i-\hat\varepsilon_{i-1}\).
acf() trace les autocorrélations aux différents retards.
update(m, . ~ . + v) réajuste le modèle en ajoutant
v. Remarque : na.omit() a supprimé des jours,
donc quelques résidus « consécutifs » ne correspondent pas à des jours
consécutifs : c’est une approximation acceptable ici.
★ Structure spatiale (diapo 84). On simule 100 sites où un gradient est-ouest a été oublié, puis on cartographie les résidus :
sx <- runif(100); sy <- runif(100); x_sp <- rnorm(100)
y_sp <- 2 + x_sp + 3 * sx + rnorm(100, 0, 0.5)
r_sp <- rstudent(lm(y_sp ~ x_sp))
plot(sx, sy, pch = 19, cex = 0.5 + abs(r_sp), col = ifelse(r_sp > 0, "red", "blue"),
xlab = "longitude", ylab = "latitude", main = "résidus : rouge > 0, bleu < 0")Notions. \(\hat y_i=\sum_j h_{ij}y_j = h_{ii}y_i+\sum_{j\ne i}h_{ij}y_j\) : \(h_{ii}\) mesure le poids de l’observation \(i\) dans sa propre prévision. \(\sum_i h_{ii} = p\). Un point est levier (atypique en \(X\)) si \(h_{ii}>2p/n\), ou \(3p/n\) (pour \(p>6\) et \(n-p>12\)), ou \(h_{ii}>0.5\).
sum(h) # = p
seuils_h <- c(2 * p_m / n_m, 3 * p_m / n_m, 0.5)
plot(h, type = "h", lwd = 2, ylim = c(0, max(0.55, h)), ylab = "h_ii", xlab = "indice")
abline(h = seuils_h, lty = c(2, 3, 1), col = c("orange", "red", "darkred"))
which(h > 2 * p_m / n_m)
ozone[which(h > 2 * p_m / n_m), c("O3", "T12", "Vx", "Ne12")] # pourquoi sont-ils atypiques ?Notions. Trois notions distinctes : aberrant (en \(Y\), \(t_i^*\)), levier (en \(X\), \(h_{ii}\)), influent (sur les coefficients, distance de Cook) :
C_i = \frac{\|\hat Y-\hat Y_{(i)}\|^2}{p\,\hat\sigma^2} = \frac{\hat\varepsilon_i^2}{p\,\hat\sigma^2}\cdot\frac{h_{ii}}{(1-h_{ii})^2}.Seuils usuels : \(C_i>4/n\) (alerte), \(C_i>1\) (préoccupant). Conduite à tenir : analyse de sensibilité, pas de suppression automatique.
C_formule <- e^2 / (p_m * s^2) * h / (1 - h)^2
all.equal(C_formule, cooks.distance(reg_mult))
# La définition : on réajuste sans i et on compare les prévisions de TOUS les individus
C_def <- sapply(1:n_m, function(i) {
m_i <- lm(O3 ~ T12 + Vx + Ne12, data = ozone[-i, ])
sum((fitted(reg_mult) - predict(m_i, newdata = ozone))^2) / (p_m * s^2)
})
all.equal(C_def, C_formule, check.attributes = FALSE)
plot(C_formule, type = "h", lwd = 2, ylab = "distance de Cook", xlab = "indice")
abline(h = 4 / n_m, lty = 2, col = "red")
# Analyse de sensibilité : on retire le point le plus influent
i_max <- which.max(C_formule)
round(cbind(avec = coef(reg_mult), sans = coef(update(reg_mult, data = ozone[-i_max, ]))), 3)Les trois situations, sur un petit jeu simulé.
set.seed(3)
x_b <- runif(20, 0, 10)
y_b <- 2 + x_b + rnorm(20)
cas <- list(aberrant = c(5, 12), # x au centre, y très loin de la droite
levier = c(25, 27), # x extrême, mais sur la droite y = 2 + x
influent = c(25, 5)) # x extrême ET hors de la droite
par(mfrow = c(1, 3))
for (k in names(cas)) {
xx <- c(x_b, cas[[k]][1]); yy <- c(y_b, cas[[k]][2])
mm <- lm(yy ~ xx)
plot(xx, yy, pch = 19, col = c(rep("grey40", 20), "red"), main = k, xlab = "x", ylab = "y")
abline(lm(y_b ~ x_b), lty = 2) # droite sans le point
abline(mm, col = "red", lwd = 2) # droite avec le point
cat(sprintf("%-9s t* = %6.2f h = %4.2f Cook = %5.2f\n", k,
rstudent(mm)[21], hatvalues(mm)[21], cooks.distance(mm)[21]))
}
par(mfrow = c(1, 1))À vous. Associer chaque panneau à la ligne du tableau de la diapo 89.
Notions. Colinéarité parfaite :
\(X'X\) non inversible (R renvoie
NA). Colinéarité forte : \(V(\hat\beta)=\sigma^2(X'X)^{-1}\)
explose (IC larges, signes instables, tests individuels non
significatifs mais test F global significatif).
VIF_j = \frac{1}{1-R_j^2},\qquad R_j^2 : R^2 \text{ de la régression de } X_j \text{ sur les autres explicatives.}\(VIF_j=1\) : aucune colinéarité ; \(>5\) : à surveiller ; \(>10\) : problématique.
expl <- c("T9", "T12", "T15", "Ne9", "Ne12", "Ne15", "Vx", "maxO3v")
m_col <- lm(reformulate(expl, response = "maxO3"), data = oc)
vif_manuel <- sapply(expl, function(v) {
r2 <- summary(lm(reformulate(setdiff(expl, v), response = v), data = oc))$r.squared
1 / (1 - r2)
})
round(vif_manuel, 2)
car::vif(m_col)
round(cor(oc[, c("T9", "T12", "T15")]), 3)Explication.
reformulate(c("a","b"), response = "y") fabrique la formule
y ~ a + b à partir de chaînes de caractères : pratique pour
boucler sur des variables. setdiff(expl, v) retire
v de la liste.
Le paradoxe typique, simulé.
set.seed(4)
n_c <- 30
x1c <- rnorm(n_c)
x2c <- x1c + rnorm(n_c, 0, 0.05) # quasi copie de x1c
y_c2 <- 1 + x1c + x2c + rnorm(n_c)
m_c2 <- lm(y_c2 ~ x1c + x2c)
summary(m_c2) # aucun coefficient significatif... F très significatif
car::vif(m_c2)
# Instabilité : on retire un individu à la fois
round(t(sapply(1:5, function(i) coef(lm(y_c2[-i] ~ x1c[-i] + x2c[-i])))), 2)
# Colinéarité parfaite : violation de H1
x3c <- 2 * x1c
coef(lm(y_c2 ~ x1c + x3c)) # NAQue faire ? Retirer une variable redondante (ici garder une seule température), les combiner (score, ACP), ou utiliser une régression pénalisée (ridge, qui accepte un peu de biais pour beaucoup moins de variance : cf. 2.6).
Notions. Trop peu de variables ⇒ biais ; trop ⇒ variance, sur-ajustement. Critères pénalisés à minimiser :
AIC = n\ln\!\Big(\frac{SCR}{n}\Big) + 2p, \qquad BIC = n\ln\!\Big(\frac{SCR}{n}\Big) + p\ln n .vars <- c("T6", "T9", "T12", "T15", "T18", "Ne6", "Ne9", "Ne12", "Ne15", "Ne18", "Vx", "maxO3v")
m_complet <- lm(reformulate(vars, response = "maxO3"), data = oc)
n_oc <- nrow(oc)
aic_cours <- function(m) { nn <- length(resid(m)); p <- length(coef(m)); nn * log(sum(resid(m)^2) / nn) + 2 * p }
bic_cours <- function(m) { nn <- length(resid(m)); p <- length(coef(m)); nn * log(sum(resid(m)^2) / nn) + p * log(nn) }
c(AIC_cours = aic_cours(m_complet), extractAIC = extractAIC(m_complet)[2])
c(BIC_cours = bic_cours(m_complet), extractAIC_k = extractAIC(m_complet, k = log(n_oc))[2])
AIC(m_complet) - aic_cours(m_complet) # AIC() ajoute une constante : mêmes comparaisons
m_aic <- step(m_complet, direction = "both", trace = 0) # pénalité 2
m_bic <- step(m_complet, direction = "both", k = log(n_oc), trace = 0) # pénalité ln(n)
formula(m_aic)
formula(m_bic)Explication. extractAIC() utilise
exactement la formule du cours ; AIC() utilise \(-2\ell+2(p+1)\) (avec \(\sigma^2\) compté comme paramètre), ce qui
ajoute une constante sans changer le classement des modèles.
step() part du modèle complet et ajoute/retire une variable
à la fois tant que le critère diminue ; k = log(n)
transforme l’AIC en BIC.
Précaution : l’inférence après sélection est trop optimiste.
set.seed(5)
bruit_df <- as.data.frame(matrix(rnorm(100 * 30), 100, 30))
bruit_df$y <- rnorm(100) # y indépendant de tout
m_b <- step(lm(y ~ ., data = bruit_df), trace = 0)
round(summary(m_b)$coefficients, 3) # des « effets » apparaissent dans du pur bruitEn recherche clinique, les variables d’ajustement se choisissent d’abord a priori (facteurs de confusion connus), pas par un algorithme.
Notions. On explique O3 par le secteur
du vent (facteur à \(I=4\) modalités).
Le facteur \(A\) est remplacé par son
codage disjonctif complet \(A_c=(\mathbf 1_{NORD},\mathbf 1_{SUD},\mathbf
1_{EST},\mathbf 1_{OUEST})\) et
Y = \mu\mathbf 1 + A_c\alpha + \varepsilon \quad\Longleftrightarrow\quad y_{ij} = \mu + \alpha_i + \varepsilon_{ij},\quad j=1,\dots,n_i .anova_ex <- data.frame(
O3 = c(64, 90, 79, 81, 88, 68, 139, 78, 114, 42),
vent = factor(c("EST", "NORD", "EST", "NORD", "OUEST", "SUD", "EST", "NORD", "SUD", "OUEST"),
levels = c("NORD", "SUD", "EST", "OUEST"))
)
anova_ex
table(anova_ex$vent) # effectifs n_i des cellules
split(anova_ex$O3, anova_ex$vent) # tableau regroupé de la diapo 96
tapply(anova_ex$O3, anova_ex$vent, mean) # moyennes par cellule
stripchart(O3 ~ vent, data = anova_ex, vertical = TRUE, pch = 19, method = "jitter")
Ac <- model.matrix(~ vent - 1, data = anova_ex) # codage disjonctif complet
AcExplication. factor(..., levels = ...)
fixe l’ordre des modalités (sinon ordre alphabétique).
split() regroupe les valeurs par modalité. Dans
model.matrix(~ vent - 1), le - 1 supprime la
constante et fait apparaître une indicatrice par modalité.
Notions. \(\mathbf 1 =
\mathbf 1_{NORD}+\mathbf 1_{SUD}+\mathbf 1_{EST}+\mathbf
1_{OUEST}\) : la matrice \((\mathbf 1
\mid A_c)\) n’est pas de plein rang, le modèle n’est pas
identifiable (\(\mu+\alpha_i =
\tilde\mu+\tilde\alpha_i\) pour une infinité de choix). On impose
une contrainte identifiante : - \(\alpha_1=0\) : une cellule de
référence (défaut de R, contr.treatment) ;
- \(\sum_i\alpha_i=0\) : effets
différentiels (contr.sum) ; - pas de
constante : \(y_{ij}=m_i+\varepsilon_{ij}\), les \(m_i\) sont les moyennes de cellule.
X_sur <- cbind(1, Ac)
qr(X_sur)$rank # 4 < 5 colonnes : H1 violée
coef(lm(anova_ex$O3 ~ Ac)) # un coefficient NA
m_ref <- lm(O3 ~ vent, data = anova_ex) # référence = NORD
m_somm <- lm(O3 ~ vent, data = anova_ex, contrasts = list(vent = "contr.sum"))
m_cell <- lm(O3 ~ vent - 1, data = anova_ex)
coef(m_ref); coef(m_somm); coef(m_cell)
moy <- tapply(anova_ex$O3, anova_ex$vent, mean)
c(mu_somme_nulle = mean(moy), moyenne_generale = mean(anova_ex$O3)) # différentes : plan déséquilibré
-sum(coef(m_somm)[-1]) # alpha de la dernière modalité (OUEST) - mu
# Trois paramétrisations, un seul modèle : mêmes valeurs ajustées
all.equal(fitted(m_ref), fitted(m_somm)); all.equal(fitted(m_ref), fitted(m_cell))Explication. Avec la contrainte de référence, la
constante est la moyenne de NORD et chaque coefficient est un
écart à NORD. Avec la somme nulle, \(\mu\) est la moyenne non pondérée
des moyennes de cellules (et non la moyenne générale quand les
effectifs diffèrent) ; R n’affiche que \(I-1\) effets, le dernier se déduit par
\(\alpha_I=-\sum_{i<I}\alpha_i\).
Changer relevel(anova_ex$vent, ref = "SUD") change la
référence.
Notions. \(H_0:\alpha_1=\dots=\alpha_I\) (le facteur n’a pas d’effet) contre \(H_1\) : au moins deux diffèrent. Sous \(H_0\), \(y_{ij}=\mu+\varepsilon_{ij}\) et \(\hat Y_0=\bar y\,\mathbf 1\). Sous H3 :
F = \frac{\|\hat Y-\hat Y_0\|^2/(I-1)}{\|Y-\hat Y\|^2/(n-I)} \sim \mathcal F_{I-1,\,n-I}.Plus généralement, pour deux modèles emboîtés \(\mathcal M_0\subset\mathcal M_1\) de dimensions \(p_0<p_1\) :
F = \frac{(SCR_0-SCR_1)/(p_1-p_0)}{SCR_1/(n-p_1)} \sim \mathcal F_{p_1-p_0,\,n-p_1}.Ya <- anova_ex$O3
I <- nlevels(anova_ex$vent); n_a <- length(Ya)
Y_chap_a <- fitted(m_ref)
Y_chap_0 <- mean(Ya)
F_a <- (sum((Y_chap_a - Y_chap_0)^2) / (I - 1)) / (sum((Ya - Y_chap_a)^2) / (n_a - I))
c(F = F_a, p_valeur = pf(F_a, I - 1, n_a - I, lower.tail = FALSE))
anova(m_ref)
m_0 <- lm(O3 ~ 1, data = anova_ex) # modèle sous H0
anova(m_0, m_ref) # test entre modèles emboîtés : même F
# Sur les 50 jours du fichier ozone
m_vent <- lm(O3 ~ vent, data = ozone)
anova(m_vent)
tapply(ozone$O3, ozone$vent, mean)
boxplot(O3 ~ vent, data = ozone, col = "lightblue")Explication.
pf(..., lower.tail = FALSE) donne \(P(\mathcal F > F_{obs})\).
anova(m0, m1) calcule la statistique F des modèles
emboîtés. Avec 10 observations l’effet visible sur le graphique n’est
pas significatif (puissance très faible) ; avec 50 jours, il l’est.
À vous. Vérifier H2/H3 pour m_vent :
résidus studentisés par modalité
(boxplot(rstudent(m_vent) ~ ozone$vent)), Q-Q plot, et
bartlett.test(O3 ~ vent, data = ozone) pour l’égalité des
variances.
Notions. - Pentes parallèles :
\(y_{ij}=\mu+\alpha_i+\beta
x_{ij}+\varepsilon_{ij}\) — lm(O3 ~ vent + T12) :
\(I\) droites parallèles, \(\alpha_i\) = décalage vertical. -
Avec interaction : \(y_{ij}=\mu+\alpha_i+(\beta+\gamma_i)x_{ij}+\varepsilon_{ij}\)
— lm(O3 ~ vent * T12) : \(I\) pentes différentes. - Choix : test de
modèles emboîtés anova(parallèle, interaction).
m_par <- lm(O3 ~ vent + T12, data = ozone)
m_inter <- lm(O3 ~ vent * T12, data = ozone)
summary(m_par)$coefficients
anova(m_par, m_inter) # les pentes diffèrent-elles selon le vent ?
# Effet du vent à température constante (ajustement)
anova(lm(O3 ~ T12, data = ozone), m_par)
cols <- c("red", "blue", "darkgreen", "orange")
niv <- levels(ozone$vent)
par(mfrow = c(1, 2))
plot(O3 ~ T12, data = ozone, col = cols[ozone$vent], pch = 19, main = "pentes parallèles")
for (k in seq_along(niv)) {
gT <- seq(min(ozone$T12), max(ozone$T12), length.out = 2)
lines(gT, predict(m_par, data.frame(T12 = gT, vent = niv[k])), col = cols[k], lwd = 2)
}
legend("topleft", niv, col = cols, pch = 19, bty = "n", cex = 0.8)
plot(O3 ~ T12, data = ozone, col = cols[ozone$vent], pch = 19, main = "avec interaction")
for (k in seq_along(niv)) {
gT <- seq(min(ozone$T12), max(ozone$T12), length.out = 2)
lines(gT, predict(m_inter, data.frame(T12 = gT, vent = niv[k])), col = cols[k], lwd = 2)
}
par(mfrow = c(1, 1))
table(ozone$vent) # effectifs : prudence avec 4 pentes sur 50 joursExplication. Dans predict(), on fournit
vent sous forme de chaîne de caractères : R la convertit
dans les modalités du facteur utilisé lors de l’ajustement.
cols[ozone$vent] associe une couleur à chaque modalité (un
facteur est codé en interne par des entiers 1, …, I).
Lecture clinique. C’est le modèle de l’ajustement : comparer des groupes (\(\alpha_i\)) à covariable constante (\(\beta\)). Tester l’interaction, c’est chercher une modification d’effet selon le sous-groupe — avec la prudence qu’impose un faible effectif par groupe.
On regroupe tous les outils de la partie 3 dans une fonction.
diagnostic <- function(m, temps = seq_along(resid(m))) {
ts <- rstudent(m); hh <- hatvalues(m); cd <- cooks.distance(m)
nn <- length(ts); p <- length(coef(m))
op <- par(mfrow = c(2, 3), mar = c(4, 4, 2, 1)); on.exit(par(op))
plot(fitted(m), ts, pch = 20, main = "t* vs ajustés", xlab = "ajustés", ylab = "t*")
abline(h = c(-1, 1) * qt(0.975, nn - p - 1), lty = 2, col = "red")
plot(fitted(m), abs(ts), pch = 20, main = "|t*| + lowess", xlab = "ajustés", ylab = "|t*|")
lines(lowess(fitted(m), abs(ts)), col = "red", lwd = 2)
qqnorm(ts, pch = 20, main = "Q-Q plot"); qqline(ts, col = "red")
plot(temps, ts, type = "b", pch = 20, main = "t* vs temps", xlab = "temps", ylab = "t*")
abline(h = 0, lty = 2)
plot(hh, type = "h", main = "leviers", ylab = "h_ii"); abline(h = c(2, 3) * p / nn, lty = 2:3, col = "red")
plot(cd, type = "h", main = "Cook", ylab = "C_i"); abline(h = 4 / nn, lty = 2, col = "red")
Xm <- model.matrix(m)[, -1, drop = FALSE]
vif <- if (ncol(Xm) > 1) sapply(seq_len(ncol(Xm)), function(j)
1 / (1 - summary(lm(Xm[, j] ~ Xm[, -j]))$r.squared)) else NA
if (ncol(Xm) > 1) names(vif) <- colnames(Xm)
e <- resid(m)
list(DW = round(sum(diff(e)^2) / sum(e^2), 3), VIF = round(vif, 2),
aberrants = which(abs(ts) > qt(0.975, nn - p - 1)),
leviers = which(hh > 2 * p / nn), influents = which(cd > 4 / nn))
}Explication. on.exit(par(op)) restaure
les paramètres graphiques à la sortie de la fonction. Pour un facteur,
le VIF calculé ici porte sur chaque indicatrice ;
car::vif() donne le VIF généralisé (GVIF), plus adapté.
ozone.txtQuestion : quels facteurs météorologiques expliquent le pic d’ozone, et quelle prévision pour une journée donnée ?
# 1. Modèle maximal, choisi a priori
m_max <- lm(O3 ~ T12 + Ne12 + Vx + O3v + vent + nebulosite, data = ozone)
summary(m_max)
car::vif(m_max)
tapply(ozone$Ne12, ozone$nebulosite, range) # nebulosite n'est qu'un recodage de Ne12 !
# 2. Diagnostics
diagnostic(m_max)
# 3. Sélection (BIC) et comparaison emboîtée
m_sel <- step(m_max, k = log(nrow(ozone)), trace = 0)
formula(m_sel)
anova(m_sel, m_max)
# 4. Modèle retenu : interprétation
summary(m_sel)
confint(m_sel)
diagnostic(m_sel)
# 5. Prévision pour une nouvelle journée
jour <- data.frame(T12 = 22, Ne12 = 3, Vx = 5, O3v = 90, vent = "OUEST", nebulosite = "SOLEIL")
predict(m_sel, jour, interval = "confidence")
predict(m_sel, jour, interval = "prediction")À vous. (a) Pourquoi nebulosite et
Ne12 ne doivent-elles pas figurer ensemble ? (b)
Interpréter chaque coefficient de m_sel « toutes choses
égales par ailleurs », avec son IC. (c) Les p-valeurs de
m_sel sont-elles fiables après sélection (3.9) ? (d)
Commenter l’écart entre IC et IP.
ozone_complet en apprentissage
/ testLe \(R^2\) mesure l’ajustement sur les données qui ont servi à estimer. Pour juger une prévision, on estime sur une partie des données et on évalue sur l’autre (ISLR, chap. 5).
set.seed(2026)
idx <- sample(nrow(oc), round(0.7 * nrow(oc)))
app <- oc[idx, ]
test <- oc[-idx, ]
modeles <- list(
simple = maxO3 ~ T12,
meteo = maxO3 ~ T12 + Ne12 + Vx,
meteo_veille = maxO3 ~ T12 + Ne12 + Vx + maxO3v,
selection_bic = formula(m_bic),
complet = reformulate(vars, response = "maxO3"),
interactions = reformulate(paste0("(", paste(vars, collapse = " + "), ")^2"), response = "maxO3")
)
res <- t(sapply(modeles, function(f) {
m <- lm(f, data = app)
c(p = length(coef(m)),
R2 = summary(m)$r.squared,
R2aj = summary(m)$adj.r.squared,
BIC = BIC(m),
RMSE_app = sqrt(mean(resid(m)^2)),
RMSE_test = sqrt(mean((test$maxO3 - predict(m, test))^2)))
}))
round(res, 3)Avec environ 950 jours d’apprentissage, même le modèle à 79 coefficients se comporte correctement en test. Réduisons l’échantillon d’apprentissage à 150 jours (moyenne sur 50 tirages) :
petit_app <- replicate(50, {
idx_p <- sample(nrow(oc), 150)
sapply(modeles[c("meteo_veille", "complet", "interactions")], function(f) {
m <- lm(f, data = oc[idx_p, ])
sqrt(mean((oc$maxO3[-idx_p] - predict(m, oc[-idx_p, ]))^2))
})
})
round(rowMeans(petit_app), 2) # RMSE de test moyenÀ vous. (a) Avec 950 jours, classer les modèles
selon le \(R^2\), le BIC et le RMSE de
test : les trois critères sont-ils d’accord ? (b) Avec 150 jours, que
devient le modèle avec interactions ? Relier au compromis biais–variance
(2.6) : 79 coefficients estimés sur 150 jours ont une très grande
variance. (c) Quel est l’apport de maxO3v (3.5) ?
Explication. sample() tire 70 % des
lignes pour l’apprentissage. (a + b + c)^2 dans une formule
génère tous les effets principaux et toutes les interactions deux à
deux. Le RMSE (racine de l’erreur quadratique moyenne) de test estime
l’erreur de prévision sur de nouveaux jours.
Chaque appel de cas_mystere() fabrique un jeu de données
présentant un problème (ou aucun). Ajustez
y ~ x1 + x2, appliquez diagnostic(), et
identifiez le problème avant de regarder la
solution.
genere_cas <- function(type, n = 60) {
x1 <- runif(n, 0, 10); x2 <- runif(n, 0, 10)
groupe <- factor(sample(c("A", "B"), n, replace = TRUE))
eps <- rnorm(n, 0, 2)
if (type == "colinearite") x2 <- x1 + rnorm(n, 0, 0.3)
if (type == "autocorrelation") eps <- as.numeric(arima.sim(list(ar = 0.85), n = n, sd = 1.2))
y <- 5 + 2 * x1 + x2 + eps
if (type == "non_linearite") y <- 5 + 0.3 * x1^2 + x2 + eps
if (type == "heteroscedasticite") y <- 5 + 2 * x1 + x2 + eps * x1 / 2
if (type == "aberrant") { i <- sample(n, 1); y[i] <- y[i] + sample(c(-1, 1), 1) * 15 }
if (type == "levier_influent") { i <- sample(n, 1); x1[i] <- 25; y[i] <- 5 + x2[i] + rnorm(1, 0, 2) }
if (type == "variable_oubliee") y <- y + 6 * (groupe == "B")
data.frame(temps = 1:n, x1, x2, groupe, y)
}
types_cas <- c("aucun", "non_linearite", "heteroscedasticite", "aberrant", "levier_influent",
"colinearite", "autocorrelation", "variable_oubliee")
cas_mystere <- function() {
tp <- sample(types_cas, 1)
d <- genere_cas(tp)
attr(d, "solution") <- tp
d
}
# Utilisation
d_m <- cas_mystere()
m_m <- lm(y ~ x1 + x2, data = d_m)
summary(m_m)
diagnostic(m_m, temps = d_m$temps)
boxplot(rstudent(m_m) ~ d_m$groupe) # penser aux variables non incluses !
# attr(d_m, "solution") # à décommenter après avoir concluExplication. genere_cas() part du
modèle correct \(y=5+2x_1+x_2+\varepsilon\) puis introduit
un seul défaut selon type.
arima.sim(list(ar = 0.85), ...) génère des erreurs
autorégressives d’ordre 1 (le cas du test de Durbin-Watson).
attr(d, "solution") cache la réponse dans un attribut de
l’objet.
Signatures à rechercher.
| Problème | Ce que l’on voit | Diapos | Remède possible |
|---|---|---|---|
| non-linéarité | courbe (U) dans \(t^*\) vs ajustés ou vs \(x_1\) | 15, 81 | terme \(x_1^2\), transformation |
| hétéroscédasticité | cône, lowess croissant de \(\lvert t^*\rvert\) | 77–79 | transformation de \(Y\) (log), moindres carrés pondérés |
| aberrant | un \(\lvert t^*\rvert\) très au-delà du seuil de Bonferroni | 74–75 | vérifier la donnée, sensibilité |
| levier influent | \(h_{ii}\) élevé et Cook élevé | 86–89 | vérifier la donnée, sensibilité |
| colinéarité | VIF ≫ 10, IC larges, F global significatif | 90 | retirer / combiner des variables |
| autocorrélation | DW ≪ 2, « vagues » dans \(t^*\) vs temps | 83 | variable retardée, modèles de séries temporelles |
| variable oubliée | résidus décalés selon groupe |
81–82 | ajouter la variable (ANCOVA, partie 4) |
Pour s’entraîner autrement.
genere_cas("heteroscedasticite", n = 500) permet de voir ce
que devient chaque signature quand \(n\) augmente ;
sapply(types_cas, function(tp) diagnostic(lm(y ~ x1 + x2, genere_cas(tp)))$DW)
compare une statistique d’un cas à l’autre.
| Notion | Diapos | Calcul à la main | Fonction R |
|---|---|---|---|
| Estimateurs MC | 19–22, 52 | solve(crossprod(X), crossprod(X, Y)) |
lm(), coef() |
| Matrice du plan | 48 | cbind(1, x) |
model.matrix() |
| Valeurs ajustées, résidus | 30–31 | X %*% b, y - X %*% b |
fitted(), resid() |
| \(\hat\sigma\) | 32 | sqrt(sum(e^2) / (n - p)) |
sigma() |
| \(V(\hat\beta)\) | 28, 54 | sigma^2 * solve(crossprod(X)) |
vcov() |
| Tests \(t\), IC | 39–40 | pt(), qt() |
summary(), confint() |
| IC / IP de prévision | 41 | formules 1.10 | predict(..., interval = ) |
| \(R^2\), \(R^2_{aj}\) | 42 | SCE / SCT |
summary()$r.squared |
| Log-vraisemblance | 65–67 | optim() sur dnorm(log = TRUE) |
logLik() |
| \(h_{ii}\) | 86–88 | diag(X %*% solve(crossprod(X)) %*% t(X)) |
hatvalues() |
| \(t_i\), \(t_i^*\) | 71–73 | formules 3.1 | rstandard(), rstudent() |
| Q-Q plot | 76 | qnorm(ppoints(n)) |
qqnorm(), qqline() |
| Durbin-Watson | 83 | sum(diff(e)^2) / sum(e^2) |
lmtest::dwtest() |
| Distance de Cook | 89 | formule 3.7 | cooks.distance() |
| VIF | 90 | 1 / (1 - R2_j) |
car::vif() |
| AIC / BIC | 91 | n * log(SCR / n) + k * p |
extractAIC(), step() |
| ANOVA, modèles emboîtés | 101–104 | formule F 4.3 | anova() |
| Contraintes identifiantes | 98–100 | — | contrasts =, relevel() |
| ANCOVA | 105 | — | lm(y ~ facteur * x) |