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

modifier

La 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 ~ x pour un modèle affine (droite avec ordonnée à l'origine) du type y = a + bx ;
y ~ x + 0 pour 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
 
Tracé obtenu.

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

modifier

Une 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

modifier

Nous 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

modifier

Au 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")

  • donnees_X est un data.frame ;
  • la valeur de l'attribut level est le niveau de confiance (0.9 pour 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

modifier

Une 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 + bu12 + ε
# ********************
# 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
 
Régression linéaire pour un modèle à deux variables.

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 + β1x1 + β1x2 + ε

alors les vecteurs x1, x2 et y doivent avoir même dimension et

y[i ] = β0 + β1x1[i ] + β1x2[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