Manière simple.
Considérons une loi affine à deux dimensions,
Le plus simple pour faire la régression consiste à utiliser la commande codice_1 (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
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" :
La régression linéaire se fait avec la fonction codice_4 "(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
avec un écart quadratique de . 
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.
Considérons deux jeux de données "u"1 et "u"2 que nous supposons corrélées par une loi affine
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 codice_5 :
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 codice_6 :
modele <- line(u1, u2)
print(modele)
Cette fonction utilise la méthode médiane-médiane, dite aussi « droite robuste de Tukey ».
Par exemple :
u1 <- seq(0, 1, 0.01) # abscisse, variable explicative
u2 <- 2*u1 + rnorm(u1, 0, 0.1) # ordonnée, variable expliquée
coefcorr <- cor(u1, u2)
modele <- line(u1, u2)
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 :
Nous pouvons extraire les paramètres avec codice_7, qui est un vecteur, ou la fonction codice_8 qui retourne également un vecteur. Cela permet tracer les données en superposant la droite aux données :
plot(u1, u2)
abline(modele)
On peut extraire les coefficients de la manière suivante :
a <- coef(modele)[1]
b <- coef(modele)[2]
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 codice_9, 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 codice_10. 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 :
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 %) :
alpha <- 0.05 # risque alpha (intervalle de confiance 1 - alpha)
sommaire <- summary(modele)
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))
ecart_type_a <- sommaire$coefficients[1,2]
student_a <- sommaire$coefficients[1,4]
a_non_nul = student_a <= (alpha/2)
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 %) :
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 :
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 "bi", avec
on peut donc calculer le coefficient de student
la fonction codice_11 nous renvoyant la valeur de 1-α/2, ou de α/2 avec l'option codice_12. 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ù
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 :
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 codice_17. Par exemple, pour un modèle de type
u1 <- seq(0, 1, 0.01) # abscisse, variable explicative
u2 <- 2*u1^2 + rnorm(u1, 0, 0.05) # ordonnée, variable expliquée
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.
Le modèle linéaire peut être à plusieurs variables. Le codice_13 utilisé doit avoir une colonne par variable. Considérons par exemple un modèle à deux variables
alors les vecteurs "x"1, "x"2 et "y" doivent avoir même dimension et 
Par exemple :
n <- 5 
x1 <- seq(0, 1, length.out = n) 
x2 <- seq(0, 1, length.out = n)
x1.data <- rep(x1, each = n) 
x2.data <- rep(x2, n)
y = 1 + 2*x1.data + 3*x2.data + rnorm(n^2, 0, 0.4)
donnees <- data.frame(x1.data, x2.data, y)
modele <- lm(y ~ ., donnees)
print("Modèle des moindres carrés : ")
print(modele)
VT <- persp(x1, x2, matrix(y, n, n))
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