C'est de la data et mon expeRtise afin d'en tirer le maximum

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 :

$$y=\alpha + \beta_1x\; +\;\beta_2x^2 + …+\; \beta_kx^k + \epsilon$$

Bien que les modèles polynomiaux permettent de modéliser des relations de formes non-linéaires (courbure, sinusoide, 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 :

Différentes formes de relation du modèle polynomial de degré 2![] 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 : régression polynomiale cubique

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éponse en fonction d'une variable prédictive

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")
  relation linéaire

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 degrees of freedom
## Multiple R-squared: 0.7183, Adjusted R-squared: 0.709 
## F-statistic: 76.51 on 1 and 30 DF, p-value: 9.38e-10
 

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 :

👉 La régression linéaire simple avec le logiciel R  
par(mfrow=c(2,2))
plot(mod1)
diagnostic de régression du modèle linéaire

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 équivalente :

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))
  visualisation de la régression polynomiale

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)
  diagnostic de la régression polynomiale

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ètre 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
  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:

```r
# 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²")
  régression polynomiale avec R

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)
  évaluation graphique de la linéarité  

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 !

 

Poursuivez votre lecture

3 réponses

Laisser un commentaire

Votre adresse de messagerie ne sera pas publiée. Les champs obligatoires sont indiqués avec *

Bonjour !

vous venez souvent ?

Identifiez-vous pour avoir accès à toutes les fontionnalités !

Aide mémoire off'R ;)

Enregistrez vous pour recevoir gratuitement mes fiches “aide mémoire” (ou cheat sheets) qui vous permettront de réaliser facilement les principales analyses biostatistiques avec le logiciel R et pour être informés des mises à jour du site.