- Particuliers : 125 euros/j
- Etudiants et Chômeurs (sur justificatifs) : 70 euros
La régression linéaire simple avec le logiciel R
Dans cet article, tourné une nouvelle fois sur la pratique, je vous propose 10 étapes pour mener à bien une régression linéaire simple avec le logiciel R.
Pour rappel, la régression linéaire simple est une méthode statistique classique, qui est employée pour évaluer la significativité du lien linéaire entre deux variables numériques continues. Autrement dit, on utilise la régression linéaire simple lorsqu’on souhaite évaluer si deux variables numériques continues sont liées de façon significative, en faisant l’hypothèse que leur relation est de type linéaire.
Les aspects théoriques de cette méthode ont déjà été abordé dans un précédent article. Pour des compléments d’informations, je vous recommande le chapitre “Regression” du Rbook , ainsi que le chapire La régression linéaire: lien entre deux variables du livre “contes et stat R” de Lise Vaudor.
Table des matières
L'importation des données
J’ ai l’habitude d’organiser mon travail sous la forme d’un projet R. Si vous ne savais pas ce qu’est un projet R, vous pouvez consulter cet article.
Dans le working directory de ce projet R, je crée un dossier “data” dans lequel je place le jeu de données à analyser, au format csv.
Si vous ne savez pas comment organiser vos données dans un tableur, vous pouvez consulter cet article pour obtenir 12 conseils pour le faire efficacement.
Ici, à titre d’exemple, j’ai créé un projet R, nommé “36_RLS”. Dans le dossier “data”, j’ai placé le jeu de données que nous allons utilisé dans cet article. Il s’agit du jeu de données Prestige du package “car”.
Pour importer un jeu de données au format csv, placé dans le dossier “data”, j’utilise le package “here”, et sa fonction du même nom. Cette fonction permet d’écrire le chemin d’accès du jeu de données, de façon automatique..
library(here)
Prestige <- read.csv2(here("data", "Prestige.csv")) Passage du jeu de données en format tibble
Ce format “tibble” permet d’accéder à un affichage des données plus pratique que celui du data.frame, puisque le type (entier, double, facteur, chaîne de caractères, etc…) des variables est directement affiché.
library(tidyverse)
Prestige <- as.tibble(Prestige)
as.tibble(Prestige)
## # A tibble: 102 x 7
## X education income women prestige census type
##
## 1 gov.administrators 13.1 12351 11.2 68.8 1113 prof
## 2 general.managers 12.3 25879 4.02 69.1 1130 prof
## 3 accountants 12.8 9271 15.7 63.4 1171 prof
## 4 purchasing.officers 11.4 8865 9.11 56.8 1175 prof
## 5 chemists 14.6 8403 11.7 73.5 2111 prof
## 6 physicists 15.6 11030 5.13 77.6 2113 prof
## 7 biologists 15.1 8258 25.6 72.6 2133 prof
## 8 architects 15.4 14163 2.69 78.1 2141 prof
## 9 civil.engineers 14.5 11377 1.03 73.1 2143 prof
## 10 mining.engineers 14.6 11023 0.94 68.8 2153 prof
## # ... with 92 more rows Evaluation visuelle de la linéarité
La régression linéaire simple permet d’évaluer la significativité du lien linéaire entre deux variables. La forme linéaire entre le deux variables est donc pré-supposée. Autrement dit, on fait l’hypothèse que la forme de la relation entre les variables est linéaire. Néanmoins, il est préférable de vérifier si cette hypothèse est acceptable, ou non, car si ce n’est pas le cas, les résultats de l’analyse n’auront pas de sens.
Ici, nous allons tester la relation entre la variable prestige (il s’agit d’un score de prestige relatif à la profession) et la variable éducation (qui reflète le niveau d’étude).
Pour évaluer de façon visuelle, la linéarité entre deux variables, on peut utiliser la fonction “scatterplot” du package “car”.
library(car)
scatterplot(prestige~education, data=Prestige) 
La ligne en trait plein est la droite de régression linéaire (définie par la méthode des moindres carrés) entre les deux variables. La ligne centrale en pointillé est la courbe de régression locale de type lowess.
Elle indique la tendance globale entre les deux variables. Les deux lignes extérieures représentent un intervalle de confiance de la courbe lowess.
Remarque : Une visualisation très pédagogique de la méthode est disponible sur ce post du blog de Lise Vaudor.
Ici, la droite de régression est comprise dans l’intervalle de confiance de la courbe lowess, l’hypothèse de linéarité est donc acceptable.
En revanche, si on étudie la forme du lien entre la variable prestige et la variable income, on s’aperçoit que la droite de régression n’est pas comprise dans l’intervalle de confiance de la courbe lowess ; l’hypothèse de linéarité est alors plus critiquable.
scatterplot(prestige~income, data=Prestige) 
Réalisation de la régression linéaire simple
Pour déterminer la droite de régression, on ajuste un modèle linéaire simple aux données, à l’aide de la fonction “lm”, comme ceci :
prest.lm1 <- lm(prestige~education, data=Prestige) Evaluation des hypothèses de validité des résultats
Le test d’évaluation de la significativité du lien linéaire entre les deux variables est valide, si les résidus :
- sont indépendants
- sont distribués selon une loi Normale de moyenne 0
- sont distribués de façon homogènes, c’est à dire, avec une variance
constante
Evaluation de l'hypothèse d'indépendance des résidus
En général, l’hypothèse d’indépendance des résidus est validée ou rejetée en fonction du protocole expérimental. Un exemple fréquent de non indépendance se rencontre lorsque la variable prédictive (en X) est une variable indiquant le temps, comme des années ou des mois, par exemple. Dans ce cas, on observe une auto-corrélation des résidus.
On parle d’auto-corrélation des résidus lorsque, par exemple, le résidu d’un point quelconque est liée à celui du point suivant dans le tableau de données. La présence d’une auto-corrélation peut être mise en évidence par un lag plot.
acf(residuals(prest.lm1), main="prest.lm1") 
Les pointillées horizontaux, sont les intervalles de confiance du coefficient de corrélation égal à 0. Les traits verticaux représentent les coefficients de corrélation entre les résidus de chaque point et ceux des points de la ligne suivante (lag=1), ou ceux séparés de deux lignes (lag=2) etc…
Ici, le plot nous montre qu’une auto-corrélation significative est présente pour les lags 1,2,3,5,7,10 et 12.
Le test de Durbin-Watson peut être employé pour évaluer la présence d’une auto-corrélation pour un lag de valeur 1. L’hypothèse d’indépendance des résidus est rejetée lorsque la p-value du test est inférieure à 0.05.
durbinWatsonTest (prest.lm1)
## lag Autocorrelation D-W Statistic p-value
## 1 0.2752512 1.43917 0.002
## Alternative hypothesis: rho != 0 Ici, le test nous indique qu’il existe une auto-corrélation significative entre les résidus d’une ligne du tableau de données et ceux de la ligne suivante. Cela peut sembler surprenant, mais comme nous le verrons plus loin sur le “residuals vs. fitted plot”, certains résidus semblent aller deux par deux.
Evaluation de l'hypothèse de normalité des résidus
Cette hypothèse peut s’évaluer graphiquement à l’aide d’un QQplot. Si les résidus sont bien distribués le long de la droite figurant sur le plot, alors l’hypothèse de normalité est acceptée. A l’inverse, s’ils s’en écartent, alors l’hypothèse de normalité est rejetée.
plot(prest.lm1,2) 
Ici, les points sont plutôt bien alignés sur la droite.
Le test de Shapiro-Wilk peut également être employé pour évaluer la normalité des résidus. L’hypothèse de normalité est rejetée si la p-value est inférieure à 0.05.
shapiro.test(residuals(prest.lm1))
##
## Shapiro-Wilk normality test
##
## data: residuals(prest.lm1)
## W = 0.98065, p-value = 0.1406 Ici, la normalité est accepté.
En revanche, si on réalisait une régression linéaire simple entre la variable prestige et la variable income, alors l’hypothèse de normalité serait plus discutable.
prest.lm2 <- lm(prestige~income, data=Prestige)
plot(prest.lm2,2 
Alors que le QQplot est à mon sens globalement satisfaisant, le test de Shapiro-Wilk rejette l’hypothèse de normalité :
shapiro.test(residuals(prest.lm2))
##
## Shapiro-Wilk normality test
##
## data: residuals(prest.lm2)
## W = 0.97169, p-value = 0.02729 Evaluation de l'hypothèse d'homogénéité des résidus
Là encore, cette hypothèse peut se vérifier de façon visuelle, pour cela il faut réaliser un “residuals vs fitted plot”. Les “fitted” correspondent aux réponses prédites par le modèle, pour les valeurs observées de la variable prédicitive. Si on s’intéresse à la régression linéaire simple entre la variable prestige et la variable éducation, les “fitted” correspondent aux valeurs de prestige prédites par le modèle pour les valeurs d’éducation présentes dans les données.
Plus précisément, on utilise pour ce plot la racine carrée des résidus standardisés. Mais pas de panique, R fait cela tout seul 😉
plot(prest.lm1, 3) 
Ici, la courbe rouge, qui est aussi une courbe de régression locale, est globalement plate. Ceci montre que les résidus ont tendance à être répartis de façon homogène tout le long du gradient des valeurs de prestige prédites. Et donc que l’hypothèse d’homogénéité des résidus est acceptée.
Comme dit précédemment, il apparaît sur ce plot, que certains résidus semblent aller deux par deux, confirmant ainsi le défaut d’indépendance :
Quand on regarde les libellés des professions on s’aperçoit qu’effectivement certaines sont très proches comme “civil.engineers ” et “mining.engineers “, ou encore “primary.school.teachers ” et “secondary.school.teachers”..
Pour revenir à l’évaluation de l’hypothèse d’homogénéité, si la courbe rouge est globalement ascendante ou descendante, alors cela signifie que la dispersion des résidus n’est pas homogène. Et dans ce cas, l’hypothèse est rejetée. En voici un exemple.
Il est également possible d’évaluer cette hypothèse en employant le test de Breush-Pagan. L’hypothèse d’homogénéité est rejetée si la p-value est inférieure à 0.05
ncvTest(prest.lm1)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.6327545 Df = 1 p = 0.4263468 Ici, le test ne rejette pas non plus l’hypothèse d’homogénéité.
Pour information, la p-value du test de Breush-Pagan réalisées sur les données du plot précédent (celui montrant une hétéorgénéité des résidus) est de 0.0012.
Evaluation à posteriori de l'hypothèse de linéarité
Cette hypothèse peut s’évaluer sur les résidus à l’aide du plot suivant :
plot(prest.lm1,1) 
Ici, le plot nous montre que lorsque les réponses prédites par le modèle (fitted values) augmentent, les résidus restent globalement uniformément distribués de part et d’autre de 0. Cela montre, qu’en moyenne, la droite de régression, est bien adaptée aux données, et donc que l’hypothèse de linéarité est acceptable.
Dans le cas contraire, les résidus pourraient avoir tendance à être successivement tous positifs, puis tous négatifs, puis à nouveau tous positifs, comme ci dessous :
Lorsqu’on regarde la droite de régression correspondant à ces résidus, on voit bien que la linéarité n’est pas satisfaite:
Visualisation des résultats
Paramètres et tests
Si les hypothèses de linéarité, de normalité, d’homogénéité et d’indépendance des résidus sont satisfaites, alors les résultats de la régression sont valides, et on peut donc les interpréter.
Ici, en réalité, ce n’est pas le cas du modèle de régression entre le niveau de prestige et le niveau d’éducation, car nous avons vu que les résidus étaient auto-corrélés. Les lignes de code suivantes servent donc seulement à l’illustration.
La fonction “summary” permet d’accéder facilement aux résultats de la régression.
summary(prest.lm1)
##
## Call:
## lm(formula = prestige ~ education, data = Prestige)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.0397 -6.5228 0.6611 6.7430 18.1636
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10.732 3.677 -2.919 0.00434 **
## education 5.361 0.332 16.148 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.103 on 100 degrees of freedom
## Multiple R-squared: 0.7228, Adjusted R-squared: 0.72
## F-statistic: 260.8 on 1 and 100 DF, p-value: < 2.2e-16 La partie `Residuals` des résultats permet d’évaluer rapidement la normalité des résidus. Lorsque les résidus sont distribués selon une loi Normale, la médiane doit être autour de 0 (c’est le cas ici), et les valeurs absolues de Q1 (premier quartile) et Q3 (troisième quartile) doivent être proches. C’est le cas, et cela était attendu puisque l’hypothèse de normalité était validée.
La première ligne de la partie coefficients concerne l’ordonnée à l’origine, alors que la seconde ligne concerne la pente. Concernant les colonnes :
- la première colonne rapporte l’estimation des coefficients des
paramètres, - la seconde colonne l’estimation de leur erreur standard,
- la troisième colonne est la statistique T
- la dernière colonne rapporte la p-value du test évaluant l’égalité à
0 des coefficients. Si la p-value est inférieure au seuil de
significativité généralement utilisé de 0.05, alors on concluera que
le paramètre est significativement différent de 0.
Le rejet d’une ou plusieurs des hypothèses n’impacte pas l’estimation des paramètres, mais seulement l’estimation de leur erreur standard, et par conséquence la statistique et la p-value des tests T employés pour tester l’égalité à 0 de ces paramètres.
En général, seul le coefficient de la pente a vraiment un intérêt. Ici, il est égal à 5.361. Cela signifie que lorsque le niveau d’éducation augmente d’une unité, alors, le niveau du score de prestige augmente de 5.361 unités.
Le test de l’égalité à 0 de l’ordonnée à l’origine n’a, lui non plus, généralement pas d’intérêt. En revanche, celui de la pente est important. En effet, ce n’est pas parce que le coefficient de la pente est fort qu’il est significativement différent de 0. La significativité de la pente dépend de la dispersion des points autour de la droite de régression. C’est donc le niveau de la p-value qui va permettre de trancher. Une p-value inférieure à 0.05, permet de conclure à un lien linéaire significatif entre la variable réponse et la variable prédictive. Le sens de la relation est donné par le signe du coefficient : s’il est positif, la relation linéaire est croissante, s’il est négatif, alors la relation linéaire est décroissante.
Intervalles de confiance des coefficients des paramètres
confint(prest.lm1)
## 2.5 % 97.5 %
## (Intercept) -18.027220 -3.436744
## education 4.702223 6.019533 Pour rappel, l’intervalle de confiance à 95% de la pente est une étendue de valeurs qui a une probabilité de 95% de contenir la vraie pente (celle de la population).
Les données influentes
Avant de conclure définitivement sur la significativité ou pas de la relation linéaire entre la variable réponse et la variable prédictive, il peut être intéressant de rechercher s’il existe des données influentes, et quel est leur impact sur les résultats.
Le package “car” dispose d’une fonction très utile pour détecter ces données différentes, il s’agit de la fonction “influenceIndexPlot”.
Cette fonction renvoie 4 graphs :
- le premier (en partant du bas), celui des hat value, reflète l’effet de levier (ou poids) de chaque donnée sur sa propre estimation. Une donnée est considérée comme atypique lorsque cette valeur est inférieure à 0.05.
- le second plot celui des p-value de Bonferroni permet de mettre en évidence les outliers. Est considérée comme outlier une donnée ayant une p-value inférieure à 0.05.
- le troisième plot, celui des résidus studentizés permet également de mettre en évidence les outliers
- le dernier plot, celui des distance de Cook permet d’évaluer l’influence des données sur les paramètres de régression. La distance de Cook mesure le changement dans l’estimation des paramètres de régression lorsque la donnée n’est pas prise en compte par les moindres carrés. Plus la distance est élevée, plus la modification des paramètres de régression est importante. Le seuil de 1 est couramment utilisé pour considérer qu’une donnée est très influente. Vous trouverez plus de détail sur l’effet de levier et les distances de Cook là.
influenceIndexPlot(prest.lm1) 
La donnée 53 est mise en évidence par 3 des 4 plots. Les deuxième et troisième graph (en partant du bas) la désigne comme potentiel outlier. De son côté, le plot des distances de Cook montre que son influence est parmi les deux plus fortes. Néanmoins, elle est inférieure à 1, ce qui nous amène à penser que son influence sur les paramètres du modèle n’est pas vraiment problématique. La pvalue de Bonferroni (plot 2) peut être obtenue avec la fonction outlierTest.
outlierTest(prest.lm1)
## No Studentized residuals with Bonferonni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferonni p
## 53 -2.98896 0.0035306 0.36012 La p-value ajustée par la méthode de Bonferonni est très éloignée du seuil de 0.05. La donnée 53 ne peut donc pas être considérée comme outlier.
Le plot des distances de cook met également en évidence une certaine influence de la données 67 sur les paramètres du modèle de régression. Par acquis de conscience, on peut refaire tourner le modèle sans ces données 53 et 67, afin visualiser leur impact. Cela peut se faire facilement en employant la fonction “comparCeofs” toujours du package “car”.
prest.lm1bis <- lm(prestige~education, data=Prestige[-c(53,67),])
compareCoefs(prest.lm1 ,prest.lm1bis)
## Calls:
## 1: lm(formula = prestige ~ education, data = Prestige)
## 2: lm(formula = prestige ~ education, data = Prestige[-c(53, 67), ])
##
## Model 1 Model 2
## (Intercept) -10.73 -11.26
## SE 3.68 3.53
##
## education 5.361 5.417
## SE 0.332 0.318
## Comme attendu, on voit que les données 53 et 67 n’ont que peu d’influence sur les coefficients des paramètres du modèle, ainsi que sur leur erreur standard, puisque les valeurs ne varient pas beaucoup.
Predictions
Parfois, on souhaite, obtenir la réponse prédite par le modèle de régression, pour une valeur spécifique de la variable prédictive, qui n’a pas été observée. Par exemple, imaginons que je souhaite obtenir le score de prestige prédit par le modèle, pour un niveau d’éducation de 10.25, qui n’a pas été observé. Pour cela, il est nécessaire de créer un data frame contenant une variable éducation, avec les valeurs souhaitées. Ce nouveau data frame est ensuite passé dans l’argument “newdata” de la fonction “predict”.
Une seule prédiction
my_df <- data.frame(education=c(10.25))
predict(prest.lm1, newdata=my_df)
## 1
## 44.21701 Il est possible également possible d‘obtenir facilement les intervalles de confiance et de prédiction, pour cette valeur prédite.
predict(prest.lm1, newdata=my_df, interval="prediction")
## fit lwr upr
## 1 44.21701 26.06518 62.36885
predict(prest.lm1, newdata=my_df, interval="confidence")
## fit lwr upr
## 1 44.21701 42.40008 46.03395 L’intervalle de confiance à 95% représente une plage de valeurs ayant une probabilité de 0.95 de contenir la vraie valeur de prestige, pour le niveau d’éducation utilisé. C’est l’intervalle de confiance de la droite de régression pour le niveau d’éducation donné.
L’intervalle de prédiction à 95% est, de son côté, une plage de valeurs ayant une probabilité de 0.95 de contenir une nouvelle valeur de prestige, qui serait observée pour le niveau d’éducation donné. Si on refaisait une enquête, si le niveau d’éducation était à 10.25, alors la probabilité que l’intervalle de prédiction contiennent cette nouvelle observation en terme de prestige, serait de 0.95.
Plusieurs prédictions
Pour obtenir plusieurs prédictions, c’est à dire pour plusieurs niveaux d’éducation, il suffit simplement de les entrer dans le vecteur :
my_df <- data.frame(education=c(10.25, 11.25, 12.25))
predict(prest.lm1, newdata=my_df,interval="confidence")
## fit lwr upr
## 1 44.21701 42.40008 46.03395
## 2 49.57789 47.75810 51.39768
## 3 54.93877 52.89190 56.98564 Manipulation des composants de la régression
Cette étape n’est pas forcément nécessaire, mais dans certaines situations, on peut avoir besoin, de manipuler certains éléments de la régression comme les résidus, les fitted, et la matrice de variance – covariance.
Les résidus
Comme déjà vu, précédemment, les résidus de la régression peuvent être obtenus avec la fonction “‘residuals”.
# on copie le jeu de données Prestige dans un nouveau jeu de données nomme my_pres
my_pres <- Prestige
# on ajoute une colonne comportant les résidus
my_pres$res <-residuals(prest.lm1)
# visualisation du nouveau jeu de données
head(my_pres)
## # A tibble: 6 x 8
## X education income women prestige census type res
##
## 1 gov.administrators 13.1 12351 11.2 68.8 1113 prof 9.25
## 2 general.managers 12.3 25879 4.02 69.1 1130 prof 14.1
## 3 accountants 12.8 9271 15.7 63.4 1171 prof 5.67
## 4 purchasing.officers 11.4 8865 9.11 56.8 1175 prof 6.31
## 5 chemists 14.6 8403 11.7 73.5 2111 prof 5.86
## 6 physicists 15.6 11030 5.13 77.6 2113 prof 4.49 Les fitted
Les fitted, correspondent, aux prédictions du modèle de régression, mais pour les valeurs observées de la variable prédictive :
# on ajoute une colonne comportant les fitted
my_pres$fitted <-fitted(prest.lm1)
# visualisation du nouveau jeu de données
head(my_pres)
## # A tibble: 6 x 9
## X education income women prestige census type res fitted
##
## 1 gov.administr~ 13.1 12351 11.2 68.8 1113 prof 9.25 59.5
## 2 general.manag~ 12.3 25879 4.02 69.1 1130 prof 14.1 55.0
## 3 accountants 12.8 9271 15.7 63.4 1171 prof 5.67 57.7
## 4 purchasing.of~ 11.4 8865 9.11 56.8 1175 prof 6.31 50.5
## 5 chemists 14.6 8403 11.7 73.5 2111 prof 5.86 67.6
## 6 physicists 15.6 11030 5.13 77.6 2113 prof 4.49 73.1 Remarque : Les fitted ajoutés au residuals sont égaux aux observations (variables prestige). Si on ne fournit pas de data frame à l’argument newdata, la fonction predict renvoie les fitted.
head(predict(prest.lm1))
## 1 2 3 4 5 6
## 59.54913 54.99238 57.72643 50.48924 67.64405 73.11215
head(fitted(prest.lm1))
## 1 2 3 4 5 6
## 59.54913 54.99238 57.72643 50.48924 67.64405 73.11215 La matrice de variance - covariance
La matrice de variance-covariance des paramètres du modèle peut s’obtenir avec la fonction “vcov” :
vcov(prest.lm1)
## (Intercept) education
## (Intercept) 13.520978 -1.1835055
## education -1.183505 0.1102162 Représentation finale de la régression
Lorsque la pente est significative, on réalise généralement une représentation graphique pour illustrer la relation linéaire qui lie les deux variables. Cette visualisation doit être en adéquation avec la régression : si une transformation a été utilisée dans le modèle de régression (log10 de la réponse par exemple), elle doit également figurer dans la représentation graphique.
En plus des données observées, et de la droite de régression, il est également de coutume de représenter l’intervalle de confiance à 95% de la pente. La fonction “geom_smooth” du package ggplot2 permet de le faire facilement, car c’est l’option par défaut. On peut également ajouter l’équation de la droite, à l’aide de la fonction “annotate“.
ggplot(Prestige, aes(y=prestige, x=education))+
geom_point()+
geom_smooth(colour="red", method="lm", fill="red") +
ylab("Prestige")+
xlab("education") +
theme_classic()+
annotate("text", x = 9, y = 80, label = "prestige = -10.73 + 5.36 * education\n (pval<0.001)") 
Il est encore possible d’ ajouter l’intervalle de prédiction sur le plot. Pour cela, il est nécessaire, au préalable, de stocker les valeurs des bornes inférieures et supérieures de cet intervalle, calculées sur les fitted. Puis de les ajouter au tableau de données.
int_pred <- predict(prest.lm1, interval="prediction")
my_prest2 <-cbind(Prestige, int_pred)
head(my_prest2)
## X education income women prestige census type fit
## 1 gov.administrators 13.11 12351 11.16 68.8 1113 prof 59.54913
## 2 general.managers 12.26 25879 4.02 69.1 1130 prof 54.99238
## 3 accountants 12.77 9271 15.70 63.4 1171 prof 57.72643
## 4 purchasing.officers 11.42 8865 9.11 56.8 1175 prof 50.48924
## 5 chemists 14.62 8403 11.68 73.5 2111 prof 67.64405
## 6 physicists 15.64 11030 5.13 77.6 2113 prof 73.11215
## lwr upr
## 1 41.33302 77.76523
## 2 36.81573 73.16903
## 3 39.52816 75.92469
## 4 32.33470 68.64379
## 5 49.31584 85.97226
## 6 54.67820 91.54609
ggplot(my_prest2, aes(y=prestige, x=education))+
geom_point()+
geom_smooth(colour="red", method="lm", fill="red") +
geom_line(aes(y=lwr), color = "red", linetype = "dashed")+
geom_line(aes(y=upr), color = "red", linetype = "dashed")+
ylab("Prestige")+
xlab("education") +
theme_classic()+
annotate("text", x = 9, y = 80, label = "prestige = -10.73 + 5.36 * education\n (pval<0.001)") Il est encore possible d’ ajouter l’intervalle de prédiction sur le plot. Pour cela, il est nécessaire, au préalable, de stocker les valeurs des bornes inférieures et supérieures de cet intervalle, calculées sur les fitted. Puis de les ajouter au tableau de données.
int_pred <- predict(prest.lm1, interval="prediction")
my_prest2 <-cbind(Prestige, int_pred)
head(my_prest2)
## X education income women prestige census type fit
## 1 gov.administrators 13.11 12351 11.16 68.8 1113 prof 59.54913
## 2 general.managers 12.26 25879 4.02 69.1 1130 prof 54.99238
## 3 accountants 12.77 9271 15.70 63.4 1171 prof 57.72643
## 4 purchasing.officers 11.42 8865 9.11 56.8 1175 prof 50.48924
## 5 chemists 14.62 8403 11.68 73.5 2111 prof 67.64405
## 6 physicists 15.64 11030 5.13 77.6 2113 prof 73.11215
## lwr upr
## 1 41.33302 77.76523
## 2 36.81573 73.16903
## 3 39.52816 75.92469
## 4 32.33470 68.64379
## 5 49.31584 85.97226
## 6 54.67820 91.54609
ggplot(my_prest2, aes(y=prestige, x=education))+
geom_point()+
geom_smooth(colour="red", method="lm", fill="red") +
geom_line(aes(y=lwr), color = "red", linetype = "dashed")+
geom_line(aes(y=upr), color = "red", linetype = "dashed")+
ylab("Prestige")+
xlab("education") +
theme_classic()+
annotate("text", x = 9, y = 80, label = "prestige = -10.73 + 5.36 * education\n (pval<0.001)") 
Pour aller plus loin
De façon générale, si vous souhaitez plus d’informations sur la régression linéaire, et sur les fonctions qui sont utilisées dans cet article, je vous recommande le livre “An R Companion to Applied Regression” de John Fox et Sanford Weisberg
Ce livre décrit avec beaucoup de détails mais de façon très pédagogique, l’utilisation du package “car”.
Voilà, j’espère que cet article vous permettra d’être plus à l’aise pour réaliser des régressions linéaires simples avec le logiciel R. Si vous avez des points bloquants, indiquez les moi en commentaire, et j’essaierai d’y répondre. Dans un prochain article, je vous parlerai des alternatives possibles lorsque les hypothèses ne sont pas satisfaites.
Si cet article vous a plu, ou vous a été utile, et si vous le souhaitez, vous pouvez soutenir ce blog en faisant un don sur sa page Tipeee 
Poursuivez votre lecture
- Comment évaluer la relation entre deux variables numériques continues
- Importer facilement vos données dans le logiciel R
- 12 conseils pour organiser efficacement vos données dans un tableur
- ANOVA à un facteur : quand les hypothèses ne sont pas satisfaites
- ANOVA à un facteur : partie 2 – la pratique
- ANOVA à un facteur : partie 1
Note : Je touche une petite commission (entre 3 et 6%) si vous passez par les liens Amazon de cet article pour acheter le livre mentionné. cela m’aide à entretenir ce blog, merci si vous le faites ! 
Merci pour cet enseignement de qualité. Je comprends mieux les choses.
Bonjour,
merci pour tout ce travail pédagogique. J’étais surpris de trouver dans les contes….. de Lise Vaudor que R² est un indicateur de la qualité d’ajustement du modèle linéaire ce qui semble largement contesté par un de vos articles. J’avoue que cela prolonge un peu ma confusion d’ignorant.
bonjour,
je suis étudiante et je n’arrive pas a bien comprendre les cours de stats.
j’ai du mal avec les modeles et les équations, par exemple:
Trouvez pour chacun des systèmes d’équations suivants un diagramme conceptuel qui puisse y correspondre.
a) M = iM a1X a2W a3XW eM
Y = iY bM c’1X c’2W c’3XW eY
b) M = iM aX eM
Y = iY b1M b2W b3Z b4MW b5MZ eY
Bonjour,
Est-ce que les conditions de validité à vérifier sont les mêmes pour les modèles de régression linéaires multiples ?
Par avance merci,
Cordialement,
Bonjour,
oui toujours les mêmes hypothèses à vérifier : indépendance, normalité et homogénéité des résidus.
Merci beaucoup pour cette réponse (ainsi que pour le blog). En revanche la relation de forme linéaire entre les variables devient impossible à vérifier dans les régressions multiples ?
Cordialement,
Antoine
Bonjour
Je construit mes graphiques avec ggplot2
je cherche à afficher sur le graphique l’équation de la regression linéaire via un objet :
regression = lm(EFR_FVC~temps, data = SCS002)$coefficients
cat(“Equation de la droitenCVF = “,round(regression[1],3),” “,round(regression[2],3),”* Temps”)
Que je colle dans le code du graphique :
plotregfinal annotate(“text”, x = 20, y = 2.8, label = cat(“Equation de la droitenCVF = “,round(regression[1],3),” “,round(regression[2],3),”* Temps”))
et cela ne marche pas
avez vous un moyen de le faire automatique ?
merci
Bonjour Benjamin,
vous trouverez une solution à la fin du chapitre 5.9.3 du livre R graphics cookbook : https://r-graphics.org/chapter-scatter#RECIPE-SCATTER-FITLINES-TEXT
Bonne continuation.
Bonjour,
Je suis débutant en stat, pouvez vous me parler de l’importance de la matrice de variance – covariance pour un modèle de regression?
Merci bcp
Bonjour,
vous trouverez des infos en bas de cette page : http://www.jybaudot.fr/Stats/covariance.html.
Bonne continuation
Bonjour
Dans le paragraphe 7 vous comparez l’influence des outliers sur les paramètres du modèle avec outlierTest et compareCoefs. J’ai plusieurs questions :
1- Si j’utilise outlierTest et que la p-value de Bonferroni est significative, est-ce assez? Dans ce cas, je fais quoi de mes outliers? Je les teste un à un ou avec compareCoefs? Je les élimine?
2- A partir de quelle degré de différence, estimez-vous que les coefficients du modèle ont été influencés par les outliers? Vous dites dans le paragraphe que c’est ok pour les valeurs 53 et 67 car elles ont fait peu varier les coefficients, mais quel seuil de variation est acceptable (10%? 15?…) ?
Merci d’avance pour vos réponses et merci pour votre travail.
Bonjour,
je rencontre une petite difficulté un peu bête. Le package car ne semble pas etre dispo pour R 3.3.3 Une autre facon d’acceder à sa fonction scatterlot ?
merci
bertrand
Bonjour,
merci beaucoup pour ce pas à pas détaillé.
Pour l’estimation graphique de linéarité j’ai tenté de retrouver avec ggplot2 la sortie de scatterplot par defaut :
ggplot(z, aes(y=valeur, x=ann))
geom_point()
geom_smooth(method=loess, method.args=list(degree=1, span=0.75), level=0.999)
geom_smooth(method=lm, se=FALSE)
Note : petit détail, il est inutile de mettre ici span à sa valeur par défaut (0.75), mais je le laisse dans la formule parce que pour une raison que j’ignore si il n’est pas introduit par method.args il n’est pas pris en compte……..
Le graphique produit est évasé aux extremités contrairement à celui fourni par scatter plot (tient compte du fait que la regression local ne se fait que “d’un coté”).
Mon intérêt pour l’usage ggplot2 ici réside dans la possibilité de facet_wrap et facet_grid.
merci encore.
Un énorme merci pour ces explications claires et bien détaillées! J’en ressors moins frustrée!
Merci Mélody, ravie si cela vous aide.
Bonne continuation.
Bonsoir,
Merci pour votre excellent travail ! Sauriez-vous me conseiller un livre pour MANOVA?
Je maitrise bien one-way, un peu moins bien two-ways mais chercherais a allez plus loins !
Merci d’avance!
Bonjour,
vous me posez une colle, j’avoue que j’ai très peu utilisée cette méthode.
Peut être cette page : https://www.r-bloggers.com/multiple-analysis-of-variance-manova/
Bonne continuation.
Bonjour, Claire
Merci pour ce poste très instructif, pourriez vous faire un poste dans le même style pour la régression non linéaire et les modèle à effet mixte
Merci encore
Bonjour,
Bravo pour la qualité de l’information
Je souhaiterais faire exactement la même méthode que dans la partie 10 (représentation finale de la régression) pour mettre sur un même graphique les points, le modèle, intervalle de confiance et intervalle de prédiction
Pour des raisons expérimentales il faut que ma régression passe par 0
Habituellement sur R, il faut utiliser le script ” Y ~ -1 X” pour obtenir la pente de l’équation Y=aX sans ordonnée à l’origine
Je n’arrive pas avec le package ggplot2 à faire passer ma régression par 0
Lorsque je programme
geom_smooth(method=”lm”, formula=y~-1 x)
j’obtiens un intervalle de confiance qui me semble faux
Avez-vous une solution ?
Merci par avance
Bonjour,
là comme ça, non, je ne sais pas comment faire. Il faudrait creuser un peu le sujet…je le garde dans un coin de ma tête, pourquoi pas pour un prochain article.
Si entre temps, vosu trouvez la solution, merci d’avance de la partager.
J’ai tenté d’apporter une réponse dans cet article : Représentez un intervalle de confiance et de prédiction
Bravo et merci pour ce site clair net et précis!!!!!!!!! Les explications sont vraiment super!
Merci Renaud !
Merci pour cet excellent travail et vous souhaitons bonne continuité. cependant je fais comment pour obtenir le document pédagogique que vous avez recommandé?
Bonjour,
vous pouvez acheter le livre sur Amazon ou autre. Je ne connais pas de lien vers une version gratuite.
Bonne continuation
Bonjour,
Grand merci pour ce travail,je viens de tirer une idee dans la vérification des hypothèses pour valider le model.
mon travail est aussi semblable a la votre.
Bonjour,
Je vous remercie pour ce travail de qualité.
Je ne trouve pas votre article sur les aspects théoriques de la régression linéaire simple, est ce que vous seriez me dire ou je peux le trouver ?
Bonjour,
je pense que c’est ici : https://delladata.fr/comment-evaluer-si-deux-variables-numeriques-continues-sont-liees/
Bonne continuation
Bonjour, merci pour cet article très intéressant !
Bonjour,
Quand j’applique le outlierTest, il m’arrive parfois que me soit retourné “Bonferroni p NA”. A quoi cela est-il du ?
Et quelle est l’interprétation dans ce cas la, les valeurs sont-elles considérées comme outliers ?
Bonjour,
C’est un peu difficile de vous répondre sans voir les données …
Est ce que vous avez des données manquantes ? Quelle est la taille de l’échantillon ? Est ce que vous détectez beaucoup ou aucun outlier ?