Introduction à la régression polynomiale
Définition
La régression polynomiale est une approche statistique qui est employée pour modéliser une relation de forme non-linéaire entre la réponse (y) et la ou les variables explicatives (x).
Pour prendre en charge cette forme non-linéaire de la relation entre y et x, ces modèles de régression intègrent des polynômes dans leurs équations :
Bien que les modèles polynomiaux permettent de modéliser des relations de formes non-linéaires (courbure, sinusoïde, etc.), ils appartiennent à la famille des modèles linéaires. Car dans le terme “modèle linéaire”, l’adjectif “linéaire” fait référence aux paramètres du modèle et au fait que leurs effets sont additionnés.
En biostatistiques, les modèles polynomiaux les plus utilisés (en tout cas par moi) sont ceux de degré 2 (quadratique), et plus rarement de degré 3 (cubique), c‘est-à-dire de la forme :
\[y=\alpha + \beta_1x\; +\;\beta_2x^2 + \epsilon\]
Et :
\[y=\alpha + \beta_1x\; +\;\beta_2x^2 +\;\beta_3x^2 \epsilon\]
Remarque : la régression linéaire est une régression polynomiale de degré 1 !
Les formes modélisées par la régression polynomiale
La régression polynomiale de degré 2, permet de modéliser des relations de formes diverses :
Ces exemples sont issus du R book (Crawley, M. J. (2012). The R book. John Wiley & Sons.)
Voici encore un exemple de régression polynomiale de degré 3 :
A quoi ça sert la régression polynomiale ?
A mon sens, il y a deux grands cas d’utilisation de la régression polynomiale.
Le premier, c’est lorsqu’on souhaite réellement (pas grossièrement) évaluer la linéarité de la relation entre une réponse (y) et une variable explicative (x), ou à l’inverse évaluer une courbure.
Dans cette situation, on va ajuster un modèle de régression linéaire, puis un modèle de régression polynomiale de degré 2, et enfin, on va comparer les ajustements à l’aide d’un test F, car les modèles sont emboîtés. Si l’ajustement du modèle polynomial est meilleur, alors la linéarité est rejetée au profit de la courbure. Dans le cas inverse, c’est la courbure qui est rejetée au profit de la linéarité.
La seconde situation, c’est lorsqu’on souhaite construire un modèle de prédiction. Dans cette situation, ce que l’on recherche, c’est obtenir des prédictions précises. Et dans ce cas-là on préférera un modèle plus complexe (quadratique par exemple) qu’un modèle qui explique simplement la relation entre y et x (linéaire)
Tutoriel avec R
Nous allons ici nous placer dans un contexte d’évaluation stricte de la linéarité entre la variable `mpg` (miles par gallon) et `disp` (volume du cylindre en cubic inches).
Visualisation des données
library(tidyverse)
ggplot(mtcars, aes(y=mpg,x=disp))+
geom_point()+ scale_y_continuous(limits=c(0,35))+
theme_classic() 
Régression linéaire simple
Visualisation
La courbe peut être ajoutée à l’aide de la ligne geom_smooth(method="lm", colour="blue")
ggplot(mtcars, aes(y=mpg,x=disp))+
geom_point()+
scale_y_continuous(limits=c(0,35))+
theme_classic()+
geom_smooth(method="lm", colour="blue") 
Le modèle de régression linéaire ne semble pas aberrant.
Ajustement
mod1 <- lm(mpg~disp, data=mtcars)
summary(mod1)
##
## Call:
## lm(formula = mpg ~ disp, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.8922 -2.2022 -0.9631 1.6272 7.2305
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 29.599855 1.229720 24.070 < 2e-16 ***
## disp -0.041215 0.004712 -8.747 9.38e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.251 on 30 p-value: 9.38e-10degrees of freedom
## Multiple R-squared: 0.7183, Adjusted R-squared: 0.709
## F-statistic: 76.51 on 1 and 30 DF, Vérification des hypothèses de normalité et d'homogénéité des résidus
Pour plus de détails, consultez mon article dédié à la régression linéaire simple :
par(mfrow=c(2,2))
plot(mod1) 
Visuellement, la normalité et l’homogénéité des résidus semblent souffrir de quelques défauts. Nous allons employer les tests de Shapiro-Wilk (normalité) et de Breusch-Pagan (homogénéité) :
shapiro.test(residuals(mod1))
##
## Shapiro-Wilk normality test
##
## data: residuals(mod1)
## W = 0.9271, p-value = 0.03255
library(car)
ncvTest(mod1)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 2.233958, Df = 1, p = 0.1350 L’hypothèse de normalité est rejetée (p <0.05) et l’hypothèse d’homogénéité des résidus est acceptée (p >0.05).
Bien que l’hypothèse de normalité des résidus soit rejetée, nous allons considérer que l’écart n’est pas dramatique, et passer à la régression polynomiale.
Régression polynomiale de degrès 2
La forme de la relation entre `mpg` et `disp` montre une légère courbure, nous allons donc réaliser une régression polynomiale de degré 2
La régression polynomiale de degré 2, peut être réalisée à l’aide de 2 syntaxes équivalentes :
mod <- lm(y~x+I(x^2))
mod <- lm(poly(x,2, raw=TRUE)) Dans la première syntaxe, la lettre I veut dire “indicatrice”, elle permet de protéger l’équation d’opérations erronées, par R.
Dans la seconde, l’argument `raw=TRUE` permet d’obtenir une paramétrisation équivalente à celle de la première syntaxe, les résultats seront donc identiques.
Visualisation
La courbe peut être ajoutée à l’aide de la ligne geom_smooth(method="lm", colour="red", formula=y~x+I(x^2))
ggplot(mtcars, aes(y=mpg,x=disp))+
geom_point()+
scale_y_continuous(limits=c(0,35))+
theme_classic()+
geom_smooth(method="lm",
colour="red",
formula=y~x+I(x^2)) 
Visuellement, l’ajustement semble satisfaisant.
Ajustement
mod2 <- lm(mpg~disp+I(disp^2), data=mtcars)
summary(mod2)
##
## Call:
## lm(formula = mpg ~ disp + I(disp^2), data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.9112 -1.5269 -0.3124 1.3489 5.3946
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.583e+01 2.209e+00 16.221 4.39e-16 ***
## disp -1.053e-01 2.028e-02 -5.192 1.49e-05 ***
## I(disp^2) 1.255e-04 3.891e-05 3.226 0.0031 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.837 on 29 degrees of freedom
## Multiple R-squared: 0.7927, Adjusted R-squared: 0.7784
## F-statistic: 55.46 on 2 and 29 DF, p-value: 1.229e-10 Nous pouvons voir l’arrivée d’une ligne `I(disp^2)` dans les résultats. Cela permet de vérifier que tout s’est bien déroulé.
mod2bis <- lm(mpg~poly(disp,2, raw=TRUE), data=mtcars)
summary(mod2bis)
##
## Call:
## lm(formula = mpg ~ poly(disp, 2, raw = TRUE), data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.9112 -1.5269 -0.3124 1.3489 5.3946
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.583e+01 2.209e+00 16.221 4.39e-16 ***
## poly(disp, 2, raw = TRUE)1 -1.053e-01 2.028e-02 -5.192 1.49e-05 ***
## poly(disp, 2, raw = TRUE)2 1.255e-04 3.891e-05 3.226 0.0031 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.837 on 29 degrees of freedom
## Multiple R-squared: 0.7927, Adjusted R-squared: 0.7784
## F-statistic: 55.46 on 2 and 29 DF, p-value: 1.229e-10
Nous pouvons voir que les résultats sont les mêmes que ceux obtenus avec la première syntaxe.
Comme précédemment, nous vérifions les hypothèses visuellement :
par(mfrow=c(2,2))
plot(mod2bis) 
Là encore nous pouvons voir des défauts de normalité et d’homogénéité des résidus.
Nous employons alors les tests :
shapiro.test(residuals(mod2))
##
## Shapiro-Wilk normality test
##
## data: residuals(mod2)
## W = 0.93679, p-value = 0.06073
ncvTest(mod2)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.3508705, Df = 1, p = 0.55362 Cette fois, le test de Shapiro ne rejette pas la normalité des résidus, et le test de Breusch-Pagan ne rejette pas non plus l’hypothèse de leur homogénéité.
Comparaison des deux modèles
Pour comparer les ajustements de ces deux modèles, nous allons employer un test F, défini par :
\[F=\frac{\left( \frac{RSS_1-RSS_2}{nb_{param_2}-nb_{param_1}} \right)}{\left(\frac{RSS_2}{n-nb_{param_2}}\right)}\]
Avec :
- l’indice 1 qui fait référence au modèle le plus simple (la régression linaire) et l’indice 2 au modèle le plus complexe (la régression polynomiale de degré 2)
- RSS (Residuals sum of squares) : la somme des carrés résiduels
- nb_param : le nombre de paramètres des modèles : 2 pour la régression linéaire (intercept et pente), 3 pour la régression polynomiale de degré 2 (intercept, pente pour x, pente pour x^2)
- n le nombre de données (ici 32)
Les sommes des carrés résiduels de chaque modèle peuvent être obtenues en employant la fonction `Anova` du package car:
Anova(mod1)
## Anova Table (Type II tests)
##
## Response: mpg
## Sum Sq Df F value Pr(>F)
## disp 808.89 1 76.513 9.38e-10 ***
## Residuals 317.16 30
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(mod2)
## Anova Table (Type II tests)
##
## Response: mpg
## Sum Sq Df F value Pr(>F)
## disp 216.932 1 26.955 1.488e-05 ***
## I(disp^2) 83.766 1 10.408 0.003104 **
## Residuals 233.393 29
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(mod1, mod2)
## Analysis of Variance Table
##
## Model 1: mpg ~ disp
## Model 2: mpg ~ disp + I(disp^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 30 317.16
## 2 29 233.39 1 83.766 10.408 0.003104 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Pour réaliser le test F, nous employons la fonction anova() (sans majuscule).
anova(mod1, mod2)
## Analysis of Variance Table
##
## Model 1: mpg ~ disp
## Model 2: mpg ~ disp + I(disp^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 30 317.16
## 2 29 233.39 1 83.766 10.408 0.003104 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Ici la p-value <0.05, l’hypothèse nulle, qui spécifie que les ajustements des deux modèles sont égaux, est rejetée. Ainsi l’ajustement du modèle de régression polynomiale de degré 2 est significativement meilleur que celui du modèle de régression linéaire. La relation entre `mpg` et `disp` est donc mieux représentée par une courbure que par une ligne droite.
Si la p-value avait été >0.05, nous aurions conclu que modéliser la courbure n’apportait pas une amélioration suffisante par rapport à la droite, et nous aurions conservé le modèle de régression linéaire.
La statistique F du test peut être retrouvée comme ceci :
# numerateur
(317.16-233.39 )/(3-2)
## [1] 83.77
# denominateur
233.39/(32-3)
## [1] 8.047931
F=((317.16-233.39 )/(3-2))/(233.39/(32-3))
F
## [1] 10.40889 Calcul de l'intervalle de prédiction
L’intervalle de confiance du modèle de régression polynomiale peut être directement représenté sur le plot, par la fonction geom_smooth().
En revanche, si nous voulons ajouter un intervalle de prédiction, il est nécessaire d’employer le code suivant :
my_mtcars <- mtcars
int_pred <- predict(mod2, interval="prediction")
my_mtcars <-cbind(my_mtcars, int_pred)
head(my_mtcars)
## mpg cyl disp hp drat wt qsec vs am gear carb fit
## Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4 22.19873
## Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4 22.19873
## Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1 25.92346
## Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1 17.02447
## Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2 14.19996
## Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1 18.49754
## lwr upr
## Mazda RX4 16.25468 28.14279
## Mazda RX4 Wag 16.25468 28.14279
## Datsun 710 19.92156 31.92535
## Hornet 4 Drive 11.00078 23.04815
## Hornet Sportabout 8.19785 20.20208
## Valiant 12.49227 24.50282
ggplot(my_mtcars, aes(y=mpg, x=disp))+
geom_point()+
geom_smooth(colour="red",
method="lm",
fill="red",
formula=y~x+I(x^2)) +
geom_line(aes(y=lwr), color = "red", linetype = "dashed")+
geom_line(aes(y=upr), color = "red", linetype = "dashed")+
theme_classic()+
annotate("text", x = 350, y = 30, label = "mpg = 3.58 - 1.05e-01 x disp +1.25e-04 X disp²") 
Remarque : Une évaluation plus rapide et sans test statistique, de la linéarité de la relation entre une variable réponse et une variable prédictive x, peut-être réalisée à l’aide de la fonction scatterplot().
Si nous l’employons avec les données `mtcars`, nous pouvons voir que la droite se retrouve en dehors des bandes de confiance.
library(car)
scatterplot(mpg~disp, data=mtcars) 
Conclusion
J’espère qu’après cette petite introduction à la régression polynomiale vous saurez comment modéliser une relation présentant une courbure, et comment évaluer si celle-ci est vraiment nécessaire, ou si une droite est suffisante !
Si cet article vous a plu ou vous a été utile, n’oubliez pas de le partager ! Vous pouvez également soutenir le blog par un don libre sur la page Tipeee.
13 réponses
Très intéressant, merci pour ce tuto.
Erreur de copier-coller pour la 4ème équation de l’exemple de la 1ère figure (courbe croissante turquoise en bas à droite) ?
Bonjour Denis,
oui c’est bien vu ! Merci.
Bonjour,
Je vous remercie pour ce cours très clair.
Comment fait-on lorsque l’on veut modéliser une interaction avec la variable transformée en polynôme? Doit-on répéter l’interaction pour chaque “ordre” de polynôme sexe*age+ sexe *age^2+sexe*age^3 par exemple?
Comment peut on l’interpréter?
Merci!
Bonjour,
je ne sais pas du tout, je n’ai jamais croisé ce cas de figure…
Navrée de ne pas pouvoir vous aider.
Bonne continuation
Bonjour,
Merci de m’éclaire le plus précisément sur 2 questions
1/Le coefficient de détermination R2 permet il de juger de la pertinence du un ajustement polynomial d’ordre 2 au même titre qu’un ajustement affine?
2/ est ce le cas aussi pour un ajustement logarithmique , puissance
Sinon ( je pense que non …) quel sens donner au coefficient R2 sur tableur Excel dans le cas d’un ajustement logarithmique…
Merci pour la considération et le temps donné à la lecture du message
Très cordialement
Bonjour,
je ne dirai pas que le R2 à lui seul permet de juger de la pertinence d’un modèle (qu’il soit ploynomial ou d’ordre 2). Pour vous en convaincre, regardez l’article : Régression linéaire simple : le R2, info ou intox ? https://delladata.fr/regression-lineaire-simple-le-r%c2%b2-info-ou-intox/).
Il me semble que la définition du R2 est toujours identique : c’est la part de dispersion expliquée par le modèle (q’il s’agit d’une régression linéaire, polynomiale d’ordre 2). Je n’ai pas compris ce que vous entendez par ajustement logarithmique.
J’espère que cela vous aide.
Bonne continuation
Bonjour,
Je vous remercie.
Article très interrasant.
J’ai une question.
Comment faire la prediction?
Par exemple si on construit un modele polynomiale simple y en fonction de x et on veut predir x_new.
Faut il le prendre x_new pour la fonction predict?
Cordialement
Bonjour,
oui c’est ça. Vous pouvez vous inspirer des prédictions réalisées dans cet article : https://delladata.fr/la-regression-lineaire-simple-avec-le-logiciel-r/
Bonne continuation
Bonjour
Il manque peut-être l’argument
smooth=list(style=”lines”)
dans le dernier graphique, pour obtenir plus exactement celui présenté.
scatterplot(mpg~disp, data=mtcars,smooth=list(style=”lines”))
Merci
E.
Bonjour à vous,
je m’intéresse à la régression polynomiale ordonnée.
J’ai donc une variable dépendante ordinale (codée as.ordered, en 3 niveaux)
J’ai une variable explicative également ordinale (codée as.ordered en 4 niveaux)
J’utilise la fonction polr du package MASS
MOD1 = polr(VD~ VE, data)
summary(MOD1)
Mais la sortie me donne 3 modalités de ma variable explicative
variable explicative.L
variable explicative.Q
variable explicative.C
Et je ne sais pas à quoi ces modalités se réfèrent ni la modalité prise en référence.
Si je recode ma variable explicative comme un facteur (as.factor)
J’ai une sortie où 3 niveaux de ma VE sont affichés (2 à 4) et le niveau 1 sert de référence. Là je comprends, mais les coefficients ne sont plus les mêmes que lorsque ma VE est codée as. ordered
Que faut-il faire?
Merci pour votre aide.
J. Pierre
Bonjour, merci pour votre article très intéressant !!
J’ai une petite question sur l’interprétation des résultats et sorties des modèles :
Comment interprétez vous les estimates du modèle ?
## disp -1.053e-01 2.028e-02 -5.192 1.49e-05 ***
## I(disp^2) 1.255e-04 3.891e-05 3.226 0.0031 **
Est ce qu’on peut dire que disp a un effet négatif ? Dans la mesure ou I(disp^2) a un estimate positif dans le summary du modèle ?
Merci pour votre réponse et encore bravo pour votre travail !
bonjour,
merci pour ce post vraiment utile! J’ai fait la même chose en utilisant la fonction glm car je devais spécifier une distribution de Poisson.
Cela a fonctionné par contre je ne parviens pas à obtenir/voir la p-value globale de mon modèle polynomial pour le reporter. Avez vous une idée?
Aussi que faire lorsque I(disp^2) n’est pas significatif alors que la comparaison des modèle me dit que l’ajustement du modèle de régression polynomiale de degré 2 est significativement meilleur?
Merci pour l’aide,
Xavier