Programmer en R/Estimer un modèle de régression linéaire

Programmer en R

• Introduction
  • Faire ses premiers pas
  • Installer R
  • Choisir son environnement de travail
  • Installer des librairies
  • Se documenter
• Manipulation et mise en forme des données
  • Manipuler les nombres
  • Manipuler les vecteurs
  • Manipuler les chaînes de caractères
  • Manipuler les facteurs
  • Manipuler les dates et les heures
  • Manipuler un tableau de données
  • Manipuler une liste
• Importation et exportation des données
  • Importer un fichier CSV
  • Résoudre les problèmes d'encodage des caractères
  • Importer un fichier TSV
  • Importer un fichier texte à largeur fixe
  • Importer des données depuis un tableur
  • Import un fichier JSON
  • Extraire des données d'une page web
• Visualisation de données
  • Comprendre la grammaire des graphiques
  • Représenter une carte
  • Visualiser une carte choroplèthe
• Modélisation statistique
  • Estimer un modèle de régression linéaire
  • Estimer un modèle de régression circulaire
• Communication
• Programmation avancée
  • Programmer une boucle
  • Programmer une fonction
• Outils mathématiques
  • Optimiser une fonction
  • Manipuler les matrices
  • Les fonctions statistiques de base
  • Les loi de probabilités, ajustement et test

Manière simple

Considérons une loi affine à deux dimensions,

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

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

u₂ = 2,024,39·u₁ - 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

Tracé obtenu.

Considérons deux jeux de données u₁ et u₂ que nous supposons corrélées par une loi affine

u₂ = a + b×u₁ + ε

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 :

u₂ = 0,030,07 + 1,939,65×u₁ + ε.

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

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

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 b_i, avec

β_i = b_i ± t_{1-α/2, n}×s_i

on peut donc calculer le coefficient de student

$t={\frac {|\beta _{i}-b_{i}|}{s_{i}}}$

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

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

où

• 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

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

u₂ = a + b⋅u₁² + ε

# ********************
# 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

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 = β₀ + β₁⋅x₁ + β₁⋅x₂ + ε

alors les vecteurs x₁, x₂ et y doivent avoir même dimension et

y[i ] = β₀ + β₁⋅x₁[i ] + β₁⋅x₂[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