Programmer en R/Estimer un modèle de régression linéaire
Considérons une loi affine à deux dimensions,
Manière simple
modifier- y = a + bx.
Le plus simple pour faire la régression consiste à utiliser la commande line (au singulier).
Commençons par générer les données :
u1 <- seq(0, 1, 0.01) # abscisse, variable explicative x
u2 <- 2*u1 + rnorm(u1, 0, 0.1) # ordonnée, variable expliquée y
La régression se fait simplement par
modele <- line(u1, u2)
print(modele)
ce qui donne
Call:
line(u1, u2)
Coefficients:
[1] -0.01580 2.02439
Le modèle est donc
- y = - 0,015,80 + 2,024,39·x
Puis, nous traçons les données et la droite de régression :
plot(u1, u2)
abline(modele)
Manière complète
modifierLa manière « complète » permet d'étendre la régression à des dimensions supplémentaires (régression multilinéaire), et donne plus d'informations.
Pour effectuer une régression, il faut créer une structure de données (data frame) et un modèle. Le modèle est une formule symbolique :
y ~ xpour un modèle affine (droite avec ordonnée à l'origine) du type y = a + bx ;y ~ x + 0pour un modèle linéaire (ordonnée à l'origine nulle) du type y = bx.
La régression linéaire se fait avec la fonction lm(formule, structure de données) (linear model).
À partir du même jeu de données
u1 <- seq(0, 1, 0.01) # abscisse, variable explicative
u2 <- 2*u1 + rnorm(u1, 0, 0.1) # ordonnée, variable expliquée
nous construisons la structure de données et effectuons la régression
donnees <- data.frame(u1, u2)
modele <- lm(u2~u1, donnees)
print(modele)
summary(modele)
ce qui donne
Call:
lm(formula = u2 ~ u1, data = donnees)
Coefficients:
(Intercept) u1
-0.0158 2.0244
Call:
lm(formula = u2 ~ u1, data = donnees)
Residuals:
Min 1Q Median 3Q Max
-0.27014 -0.06658 -0.00200 0.07255 0.33107
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.01580 0.02008 -0.787 0.433
x 2.02439 0.03469 58.352 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1016 on 99 degrees of freedom
Multiple R-squared: 0.9717, Adjusted R-squared: 0.9715
F-statistic: 3405 on 1 and 99 DF, p-value: < 2.2e-16
Le modèle établi par régression est donc
- u2 = 2,024,39·u1 - 0,015,80
avec un écart quadratique de 0,101,6.
Puis nous traçons les points et la droite de régression (ici en rouge) :
plot(u1, u2)
abline(modele, col="red") # droite y = ax + b de couleur rouge
Nous pouvons ensuite utiliser ce modèle pour calculer de nouvelles valeurs de y pour des valeurs de x données, avec un intervalle de confiance. Par exemple, pour avoir y en 0,3, 0,5 et 0,7 avec un intervalle de confiance de 0,95 (95 %) :
valeurs <- c(0.3, 0.5, 0.7)
predict(modele, data.frame(u1 = valeurs), level = 0.95, interval = "confidence")
Procédure simple
modifier
Considérons deux jeux de données u1 et u2 que nous supposons corrélées par une loi affine
- u2 = a + b×u1 + ε
où ε est l'écart, qui suit une variable aléatoire centrée sur 0 et d'écart type uniforme.
La première étape consiste à calculer le coefficient de corrélation linéaire pour vérifier notre hypothèse. Cela se fait avec la fonction cor() :
cor(u1, u2)
L'interprétation de ce coefficient dépend du domaine étudié, mais plus la valeur absolue est proche de 1, plus les données « collent » au modèle affine. On s'attache souvent à avoir une valeur absolue supérieure à 0,8 ; elle doit être en tous cas supérieure à 0,5.
La régression proprement dite se fait ensuite avec la fonction line() :
modele <- line(u1, u2)
print(modele)
Cette fonction utilise la méthode médiane-médiane, dite aussi « droite robuste de Tukey ».
Par exemple :
# ********************
# Création des données
u1 <- seq(0, 1, 0.01) # abscisse, variable explicative
u2 <- 2*u1 + rnorm(u1, 0, 0.1) # ordonnée, variable expliquée
# **************************
# Coefficient de corrélation
coefcorr <- cor(u1, u2)
# ****************************
# Régression linéaire de Tukey
modele <- line(u1, u2)
# *********
# Affichage
print(paste("Coefficient de corrélation linéaire : ",
1e-3*round(1e3*coefcorr)))
print("Modèle de Tukey : ")
print(modele)
donne (aux variations du tirage aléatoire près) :
[1] "Coefficient de corrélation linéaire : 0.986 [1] "Modèle de Tukey : " Call: line(u1, u2) Coefficients: [1] 0.03007 1.93965
Le modèle est donc ici :
- u2 = 0,030,07 + 1,939,65×u1 + ε.
Nous pouvons extraire les paramètres avec modele$coefficients, qui est un vecteur, ou la fonction coef(modele) qui retourne également un vecteur. Cela permet tracer les données en superposant la droite aux données :
# *****
# Tracé
plot(u1, u2)
abline(modele)
On peut extraire les coefficients de la manière suivante :
# Première manière
a <- coef(modele)[1]
b <- coef(modele)[2]
# Seconde manière
a <- modele$coefficients[1]
b <- modele$coefficients[2]
Procédure complète
modifierUne fois cette régression faite, il faut poursuivre par les tests de nullité : les coefficients a et b sont-ils significativement non nuls ? Cela nécessite de calculer l'estimateur de la variance sur a et b, puis, le cas échéant, de faire une régression en forçant le passage de la droite par l'origine.
Pour cela, on fait appel à la fonction lm(), qui utilise elle la méthode des moindres carrés. Cette fonction nécessite de définir le modèle utilisé sous la forme d'une « formule », une description symbolique du modèle. Une formule est introduite par le tilde ~. Les deux formules que nous utilisons ici sont :
u2 ~ u1 # modèle affine
u2 ~ u1 + 0 # modèle linéaire
La régression se fait de la manière suivante :
# ******************************
# régression des moindres carrés
donnees <- data.frame(u1, u2)
modele <- lm(u2~u1, donnees)
print("Modèle des moindres carrés : ")
print(modele)
plot(u1, u2)
abline(modele)
Au passage, on remarque que les coefficients obtenus sont légèrement différents. Puis, on effectue le test de non nullité de b et de a (coefficient 3 pour un intervalle de confiance de 99,7 %) :
# ********************
# Tests de non nullité
alpha <- 0.05 # risque alpha (intervalle de confiance 1 - alpha)
sommaire <- summary(modele)
# test sur b
student_b <- sommaire$coefficients[2,4]
b_non_nul <- student_b <= (alpha/2)
print("Test de non nullité de b.")
print(paste("Hypothèse 1, les données sont corrélées : ", b_non_nul))
# Test sur a
ecart_type_a <- sommaire$coefficients[1,2]
student_a <- sommaire$coefficients[1,4]
a_non_nul = student_a <= (alpha/2)
# Affichage
print("Test de non nullité de a.")
print(paste("Hypothèse 1, a est non nul : ", a_non_nul))
ce qui nous donne ici :
[1] "Test de non nullité de b." [1] "Hypothèse 1, les données sont corrélées : TRUE" [1] "Test de non nullité de a." [1] "Hypothèse 1, a est non nul : FALSE"
On refait donc une régression linéaire en forçant a = 0 :
if (!a_non_nul) { # si l'ordonnée à l'origine est nulle
modele <- lm(u2~u1 + 0, donnees) # régression passant par l'origine
a <- 0
b <- coef(modele)
sommaire <- summary(modele)
ecart_type_b = sommaire$coefficients[1]
} else {
a <- coef(modele)[1]
b <- coef(modele)[2]
ecart_type_b = sommaire$coefficients[2,2]
}
abline(modele) # tracé du modèle
print(modele)
Puis, on calcule les erreur standard sur a, b, et sur l'estimation des y si l'on utilise le modèle (en considérant un grand nombre de valeurs et un intervalle de confiance de 99,7 %) :
# **********************
# Estimation des erreurs
coef_erreur = qnorm(1-alpha/2)
b <- modele$coefficients[2]
erreur_b <- coef_erreur*ecart_type_b
a <- modele$coefficients[1]
erreur_a <- coef_erreur*ecart_type_a
erreur_y <- coef_erreur*summary(modele)$sigma
odgea <- round(log10(erreur_a)) # ordre de grandeur pour arrondi
odgeb <- round(log10(erreur_b))
odgey <- round(log10(erreur_y))
err_a <- 10^(odgea - 2)*round(erreur_a*10^(-odgea + 2))
err_b <- 10^(odgeb - 2)*round(erreur_b*10^(-odgeb + 2))
err_y <- 10^(odgey - 2)*round(erreur_y*10^(-odgey + 2))
print(paste("Erreur sur a : ", err_a))
print(paste("Erreur sur b : ", err_b))
print(paste("Erreur sur y au centre : ", err_y))
ce qui donne
[1] "Erreur sur a : 0.0058" [1] "Erreur sur b : 0.01" [1] "Erreur sur y au centre : 0.294"
On peut afficher différents graphiques de synthèse avec
layout(matrix(1:4,2,2)) # fenêtre de 2x2 graphiques
plot(modele)
qui comprend :
- les résidus en fonction de x ;
- les résidus normalisés en fonctions de x ;
- la droite de Henry (quantiles-quantiles) des résidus (permet de vérifier si le résidu suit une loi normale centrée) ;
- les résidus normalisés en fonction du levier (x - x).
Comparaison avec un modèle a priori
modifierNous avons vérifié si les paramètres obtenus étaient différents de zéro de manière significative. Mais on peut aussi vérifier s'ils sont différent d'une valeur a priori. La commande
pt(x, n)
donne la fonction de répartition de la loi de Student à n degrés de liberté (voir Les loi de probabilités, ajustement et test), c'est-à-dire la probabilité d'avoir une valeur inférieure ou égale à x, et la fonction
pt(x, n, lower.tail = FALSE)
calcule la probabilité complémentaire. Le paramètre βi est estimé par la valeur bi, avec
- βi = bi ± t1-α/2, n×si
on peut donc calculer le coefficient de student
la fonction pt() nous renvoyant la valeur de 1-α/2, ou de α/2 avec l'option lower.tail = FALSE. Ici, nous savons que la valeur théorique est 2, donc
bth <- 2
tb <- abs(bth - b)/ecart_type_b
pr_b_similaire_bth <- 2*pt(tb, length(u1)-2, lower.tail = FALSE)
Utilisation du modèle
modifierAu delà du calcul des paramètres du modèle, R permet d'exploiter ce modèle, c'est-à-dire qu'il calcule la valeur de y pour un x donné, et fournit avec l'intervalle de confiance. Cela se fait avec la commande
predict(modele, donnees_X, level = …, interval = "confidence")
où
- donnees_X est un
data.frame; - la valeur de l'attribut
levelest le niveau de confiance (0.9pour 90 %).
Par exemple, pour avoir la valeur de y en 0,3, 0,5 et 0,7 avec un inetrvalle d econfiance de 95 % :
> valeurs <- c(0.3, 0.5, 0.7)
> predict(modele, data.frame(u1 = valeurs), level = 0.95, interval = "confidence")
fit lwr upr
1 0.6041944 0.5930618 0.6153269
2 1.0069906 0.9884364 1.0255449
3 1.4097869 1.3838109 1.4357628
Autres fonctions utiles
modifierUne fois le modèle établi, on peut récupérer :
- les résidus avec
resid(modele); - les valeurs prédites par le modèle avec
fitted(modele).
Lorsque le modèle est linéaire pour les coefficients mais pas pour les fonctions des variables explicatives, on peut vérifier graphiquement l'adéquation du modèle avec fitted(modele). Par exemple, pour un modèle de type
- u2 = a + b⋅u12 + ε
# ********************
# Création des données
u1 <- seq(0, 1, 0.01) # abscisse, variable explicative
u2 <- 2*u1^2 + rnorm(u1, 0, 0.05) # ordonnée, variable expliquée
# ******************************
# régression des moindres carrés
donnees <- data.frame(u1, u2)
modele <- lm(u2~I(u1^2), donnees)
print("Modèle des moindres carrés : ")
print(modele)
plot(u1, u2)
lines(u1, fitted(modele))
Modèle à plusieurs variables
modifier
Le modèle linéaire peut être à plusieurs variables. Le data.frame utilisé doit avoir une colonne par variable. Considérons par exemple un modèle à deux variables
- y = β0 + β1⋅x1 + β1⋅x2 + ε
alors les vecteurs x1, x2 et y doivent avoir même dimension et
- y[i ] = β0 + β1⋅x1[i ] + β1⋅x2[i ] + ε[i ].
Par exemple :
# ********************
# Création des données
# nombre de valeurs de chaque variable explicative
n <- 5
# variables explicative
x1 <- seq(0, 1, length.out = n)
x2 <- seq(0, 1, length.out = n)
# répétition des valeurs pour mixer x1 et x2
x1.data <- rep(x1, each = n)
x2.data <- rep(x2, n)
# variable expliquée
y = 1 + 2*x1.data + 3*x2.data + rnorm(n^2, 0, 0.4)
# ******************************
# régression des moindres carrés
donnees <- data.frame(x1.data, x2.data, y)
modele <- lm(y ~ ., donnees)
print("Modèle des moindres carrés : ")
print(modele)
# ************************
# représentation graphique
# tracé des données
VT <- persp(x1, x2, matrix(y, n, n))
# traces du plan modèle
mx1 <- min(x1) ; mx2 <- min(x2)
Mx1 <- max(x1) ; Mx2 <- max(x2)
Mzx <- trans3d(Mx1, mx2, fitted(modele)[n], VT)
Myz <- trans3d(mx1, Mx2, fitted(modele)[n*(n-1)+1], VT)
lines(c(Minit[1], Mzx[1]), c(Minit[2], Mzx[2]), lwd = 2)
lines(c(Minit[1], Myz[1]), c(Minit[2], Myz[2]), lwd = 2)
donne
[1] "Modèle des moindres carrés : "
Call:
lm(formula = y ~ ., data = donnees)
Coefficients:
(Intercept) x1.data x2.data
1.269 1.717 2.966