Régression linéaire simple avec R et Python
Comparons deux scripts permettant d’obtenir les paramètres d’une régression linéaire simple stochastique, de type \(y = ax + b + \varepsilon.\). L’un avec Python, l’autre avec R.
Les données figurent sur une feuille de calcul Excel et ne sont pas reprises ici. Vous pouvez les trouver en page de RLS avec Excel ou d’erreurs d’une RLS. On compte quinze valeurs sur trois colonnes : la date, un nombre de transats loués et la recette quotidienne du loueur.
Nous devrons déterminer les paramètres de la droite de régression et leur écarts-types, l’erreur type de la régression et le coefficient de corrélation.
Python
Nous utiliserons la bibliothèque standard numpy à laquelle on ajoutera pandas. Cette dernière n’est pas la plus « statistique » de Python (ce serait plutôt statsmodels). C’est pourquoi aucune instruction ne fournit directement les coefficients et il nous faudra écrire leurs formules dans le script.
La valeur de \(y\) expliquée par le modèle, habituellement notée \(\widehat{y},\) sera nommée y_hat. À l’impression, nous retiendrons trois décimales (d’où les .3f). Nous n’avons pas demandé de nuage de points.
import pandas as pd
import numpy as np
# Lecture du fichier Excel
df = pd.read_excel("RLS.xlsx")
x = df["transats"].values
y = df["recette"].values
n = len(x)
# Régression linéaire
x_mean = np.mean(x)
y_mean = np.mean(y)
a = np.sum((x - x_mean) * (y - y_mean)) / np.sum((x - x_mean)**2)
b = y_mean - a * x_mean
# Valeurs ajustées et résidus
y_hat = b + a * x
residuals = y - y_hat
# Statistiques
sigma_hat = np.sqrt(np.sum(residuals**2) / (n - 2))
se_a = sigma_hat / np.sqrt(np.sum((x - x_mean)**2))
se_b = sigma_hat * np.sqrt(1/n + x_mean**2 / np.sum((x - x_mean)**2))
r = np.corrcoef(x, y)[0, 1]
# Résultats
print("Paramètres de la régression :")
print(f"b (intercept) = {b:.3f} (écart-type = {se_b:.3f})")
print(f"a (pente) = {a:.3f} (écart-type = {se_a:.3f})")
print(f"Erreur-type de la régression = {sigma_hat:.3f}")
print(f"Coefficient de corrélation r = {r:.3f}")
Sortie :
Paramètres de la régression :
b (intercept) = 16.921 (écart-type = 22.234)
a (pente) = 14.070 (écart-type = 1.091)
Erreur-type de la régression = 11.664
Coefficient de corrélation r = 0.963
R
Nous utiliserons presque uniquement R de base. Seul le package readxl sera chargé pour importer le jeu de données depuis Excel (on suppose qu’il a déjà été installé). La fonction lm() donne toutes les informations.
# Lecture du fichier Excel
library(readxl)
data <- read_excel("RLS.xlsx")
transats <- data$transats
recette <- data$recette
# Régression linéaire
model <- lm(recette ~ transats)
summary_model <- summary(model)
b <- coef(model)[1]
a <- coef(model)[2]
se_b <- summary_model$coefficients[1, 2]
se_a <- summary_model$coefficients[2, 2]
sigma_hat <- summary_model$sigma
r <- cor(transats, recette)
# Résultats
cat("Paramètres de la régression :\n")
cat(sprintf("b (intercept) = %.3f (écart-type = %.3f)\n", b, se_b))
cat(sprintf("a (pente) = %.3f (écart-type = %.3f)\n", a, se_a))
cat(sprintf("Erreur-type de la régression = %.3f\n", sigma_hat))
cat(sprintf("Coefficient de corrélation r = %.3f\n", r))
# Graphique
plot(transats, recette,
main = "Régression linéaire : Recette ~ Transats",
xlab = "Transats",
ylab = "Recette",
pch = 16)
abline(model, col = "red", lwd = 2)
La sortie est rigoureusement la même qu’avec Python. Quant au nuage de points…
