> Statistiques > Analyse de données > Régression linéaire
Régression linéaire
C'est un modèle de type y = b0 + b1.x1 + b2.x2 + e si par exemple on a 2 variables explicatives x1 et x2 et y est la variable dépendante. e est l'erreur sur la prédiction et il s'agit de trouver les meilleurs coefficients b0, b1 et b2 pour minimiser e
La fonction lm permet de calculer la régression linéaire d'une variable dépendante numérique en fonction de variables explicatives.
Exemple de régression linéaire toute simple :
x <- c(4, 6, 3, 5, 1, 9)
y <- 2 * x + 1 + rnorm(6, 0, 0.3)
lm(y ~ x)
donne par exemple :
Call:
lm(y ~ x)
Coefficients:
(Intercept) x
0.8384 2.0289
Appel de lm avec un data frame :
- fr <- data.frame(x = c(4, 6, 3, 5, 1, 9), y = c(8, 11, 5, 10, 3, 17), z = c("a", "a", "b", "a", "b", "a")); lin <- lm(y ~ x, fr)
- On peut tracer le graphe avec la droite par : plot(y ~ x, fr); abline(lm(y ~ x, fr)).
- on peut limiter la régression à seulement certaines observations du frame : lin <- lm(y ~ x, fr subset = z == "a")
L'objet renvoyé par lm est de la classe "lm" et a les attributs suivants :
lin <- lm(y ~ x, fr) :
- lin$coefficients (ou coefficients(lin) ou encore coef(lin)) : le vecteur des coefficients de la droite.
- lin$fitted.values (ou fitted(lin)) : les valeurs calculées par la régression d'après les variables explicatives.
- lin$residuals (ou residuals(lin)) : le vecteur des résiduels (c'est à dire, les valeurs effectives moins les valeurs calculées lin$fitted.values).
- lin$model (ou model.frame(lin) : les données d'origine.
- summary(lin) permet d'avoir les coefficients, leur deviation standard, la valeur de la statistique t (coefficient divisé par la déviation standard) ainsi que la probabilité que le coefficient soit significativement différent de 0.
- Multiple R-squared = 0.92 : 92% de la variance est expliquée par le modèle.
- adjusted R-squared : prend en compte aussi le nombre de variables (fixed effects) utilisés pour expliquer la variable dépendante (peut être beaucoup diminué si beaucoup de variables).
- la p-value est une probabilité conditionnelle : probabilité des données sous l'hypothèse nulle que le facteur n'a pas d'effet.
- la p-value de la dernière ligne est la p-value globale par rapport à la valeur de F pour les degrés de libertés indiqués, tandis que les autres (Pr(>|t|) sont celles de chaque facteur, c'est à dire la p-value sur le test que le coefficient puisse être 0 (c'est la même valeur si une seule variable explicative).
- summary(lin)$coefficients : permet d'accéder à la table des coefficients et de leur probabilité (à chaque fois, y compris pour l'intercept, la p-value est la probabilité de l'hypothèse nulle que la valeur vaut 0).
- summary(lin)$r.squared : le R2.
Intervalle de confiance sur les coefficients : on peut avoir l'intervalle de confiance sur les coefficients , avec une proba alpha = 0.95 :
- avec la valeur estimée +/- t(alpha/2) * déviation standard, où t(alpha/2) est le quantile de la distribution t à n degrés de liberté pour alpha/2 (n = nombre de points - nombre de coefficients à prédire, 2 si regression à une seule variable) : qt(0.975, n).
- directement avec confint(lin, level = 0.95)
Prédiction de valeurs :
predict(lin, newdata = data.frame(x = c(1, 2, 3, 4, 5))) : renvoie un vecteur des valeurs prédites (noms des variables explicatives doivent être les mêmes dans le data frame.
- predict(lin, data.frame(x = c(1, 2, 3, 4, 5)), se.fit = TRUE) : renvoie une liste, avec les attributs :
- fit : le vecteur des valeurs prédites.
- se.fit : le vecteur des erreurs standard sur les valeurs prédites.
- predict(lin, data.frame(x = c(1, 2, 3, 4, 5)), interval = "confidence", level = 0.99) : renvoie une matrice avec les colonnes fit (valeur prédite), lwr (valeur minimale de l'intervalle de confiance), upr (valeur maximale de l'intervalle de confiance) et level le niveau de confiance souhaité (0.95 par défaut si non indiqué)
- predict(lin, data.frame(x = c(1, 2, 3, 4, 5)), interval = "prediction", level = 0.99) : même chose, mais avec l'intervalle de prédiction, plus grand que l'intervalle de confiance.
- Exemple pour tracer les intervalles :
plot(eruptions ~ waiting, faithful, pch = 3, col = "blue")
pred <- predict(lin, newdata = data.frame(faithTest[order(faithTest$waiting),]), interval = "confidence")
matlines(faithTest[order(faithTest$waiting),], pred[, c("lwr")], col = "red")
matlines(faithTest[order(faithTest$waiting),], pred[, c("upr")], col = "red")
régression multiple : si fr est un data frame de variables x, y et z :
- lin <- lm(z ~ x + y) : régression multiple de z en fonction de x et y (termes en x et en y).
- lin <- lm(z ~ x + y - 1) : régression multiple de z en fonction de x et y (termes en x et en y), mais sans le terme constant (par défaut, il y a un terme constant)
- lin <- lm(z ~ x + y + 0) : identique à lin <- lm(z ~ x + y - 1)
- lin <- lm(z ~ x + y + 1) : identique à lin <- lm(z ~ x + y) qui contient implicitement le terme constant (ordonnée à l'origine).
- lin <- lm(z ~ x : y) : régression multiple de z en fonction du terme croisé entre x et y (terme en x * y).
- lin <- lm(z ~ x * y) : régression multiple de z en fonction de x, y et du terme croisé (termes en x, en y et en x * y). x * y est équivalent à x + y + x : y
- lin <- lm(z ~ x * y - x) : régression multiple de z en fonction de y et du terme croisé (termes en y et en x * y), le '-' retirant la dépendance par rapport au terme indiqué.
- lin <- lm(z ~ I(x * x) + y, fr) : régression multiple de z en fonction de x * x et de y, le I indiquant d'interpréter le contenu avec les opérateurs mathématiques usuels (c'est le '*' indiquant la multiplication et non le '*' indiquant une dépendance en x, y et le terme croisé).
- lin <- lm(z ~ . - z, fr) : régression multiple de z en fonction de toutes les variables ('.') sauf z (identique ici à z ~ x + y).
- predict(lin, data.frame(x = c(5, 7), y = c(3, 1))) : donne le vecteur des valeurs de z prédites en fonction des valeurs de x et y.
Traçage de différents graphes à partir du modèle : si par exemple
lin <- lm(y ~ x)
- plot(lin, which = 1) : residuals versus fitted values pour voir si il y a des biais dans les residuals
- plot(lin, which = 2) : QQ plot entre les residuals normalisés (moyenne 0, écart-type 1) et loi normale N(0, 1) pour vérifier leur normalité.
- on peut aussi faire plot(lin) pour voir les différents graphes les uns après les autres.
Détermination des variables les plus importantes pour le modèle : le principe est la minimisation du critère d'information d'Akaike (AIC, Akaike Information Criteria) qui indique l'équilibre entre la précision d'un modèle et sa complexité :
- si lin <- lm(y ~ . , fr), faire lin2 <- step(lin, trace = 0) pour trouver le modèle linéaire obtenu en supprimant des variables explicatives pour minimiser le critère d'AIC.
- avec lin2 <- step(lin, trace = 1) (défaut), on voit les modèles successifs et la valeur de l'AIC.
- on peut aussi travailler en rajoutant des termes plutôt qu'en les supprimant avec l'argument direction = "forward"
Intervalle de prédiction ou de confiance :
- intervalle de prédiction est l'intervalle dans lequel un point d'un nouvel échantillon va tomber avec une probabilité donnée.
- intervalle de prédiction :
x <- runif(50, 0, 10)
y <- 2 * x + rnorm(50, 0, 2)
plot(y ~ x)
lin <- lm(y ~ x)
abline(lin, col = "red")
newx <- seq(0, 10, by = 0.05)
prd <- predict(lin, newdata = data.frame(x = newx), interval = "prediction", level = 0.95)
lines(newx, prd[, 2], lty = 5, col = "orange");
lines(newx, prd[, 3], lty = 5, col = "orange");
- intervalle de confiance est l'intervalle dans lequel la droite de régression peut se trouver. Sa largeur diminue avec le nombre de points.
-
x <- runif(50, 0, 10)
y <- 2 * x + rnorm(50, 0, 2)
plot(y ~ x)
lin <- lm(y ~ x)
abline(lin, col = "red")
newx <- seq(0, 10, by = 0.05)
prd <- predict(lin, newdata = data.frame(x = newx), interval = "prediction", level = 0.95)
lines(newx, prd[, 2], lty = 5, col = "orange");
lines(newx, prd[, 3], lty = 5, col = "orange");
Régression linéaire passant par un point autre que l'origine, par exemple de coordonnées (x0, y0) :
-
set.seed(1); x <- runif(20, 0, 1); y <- 2 * x - 1 + rnorm(20, 0, 0.1)
plot(y ~ x)
x0 <- 1; y0 <- 1; xCent <- x - x0; yCent <- y - y0
lin <- lm(yCent ~ xCent + 0)
abline(- coef(lin)[1] * x0 + y0, coef(lin)[1], col = "red")
Copyright Aymeric Duclert
programmer en R, tutoriel R, graphes en R