Prise en compte de relations non linéaires entre les variables dans un modèle linéaire

Utiliser le package gam, c’est Jais qui en avait parlé en cours de M2.

library(package = "gam")

## Loading required package: splines
## Loaded gam 1.09.1

On va travailler sur un exemple au hasard.

airquality2 <- within(airquality, {
    Ozone <- Ozone ^ (1 / 3)
})
airquality2 <- airquality2[complete.cases(airquality2), ]

Ensuite on pose le modèle linéaire classique, en utilisant l’argument family = gaussian. On utilise la fonction s() (spline) pour indiquer qu’on souhaite transformer la variable. Attention aux données manquantes, gam n’aime pas trop.

mod1 <- gam(Ozone ~ s(Solar.R) + s(Wind) + Temp, data = airquality2, family = gaussian)

par(mfrow = c(1, 2))
termplot(mod1, rug = TRUE, se = TRUE, terms = c("s(Solar.R)", "s(Wind)"))

De base le nombre de degrès de liberté est estimé par gam mais tu peux le modifier à la main.

mod2 <- gam(Ozone ~ ns(Solar.R, df = 8) + ns(Wind, df = 2), data = airquality2, family = gaussian)

par(mfrow = c(1, 2))
termplot(mod2, rug = TRUE, se = TRUE)

La fonction summary donne deux tests : Anova for Parametric Effects et Anova for Nonparametric Effects. L’effet d’une variable transformées est découpé en deux effets : un effet paramètrique (l’effet linéaire classique) et un effet non-paramètrique (l’effet spaghetti qui s’ajoute à l’effet linéaire).

En cas de courbe en U tu peux donc avoir un effet paramètrique non significatif (en moyenne l’effet de la variable est une ligne droite) avec un effet nonparamétrique significatif (la forme de U qui s’ajoute au dessus est significative). Dans notre exemple c’est le cas pour Solar.R. Par contre pour Wind l’effet linéaire moyen est significatif (en moyenne ça descend), et en plus l’effet nonparamétrique est significatif (en plus ça descend pas droit).

Enfin les variables non transformées comme Temp n’ont du coup pas de test de l’effet non paramètrique, c’est normal.

summary(mod1)

## 
## Call: gam(formula = Ozone ~ s(Solar.R) + s(Wind) + Temp, family = gaussian, 
##     data = airquality2)
## Deviance Residuals:
##    Min     1Q Median     3Q    Max 
## -1.125 -0.306 -0.070  0.336  1.079 
## 
## (Dispersion Parameter for gaussian family taken to be 0.226)
## 
##     Null Deviance: 87.21 on 110 degrees of freedom
## Residual Deviance: 22.83 on 101 degrees of freedom
## AIC: 161.5 
## 
## Number of Local Scoring Iterations: 2 
## 
## Anova for Parametric Effects
##             Df Sum Sq Mean Sq F value  Pr(>F)    
## s(Solar.R)   1   14.9   14.87    65.8 1.2e-12 ***
## s(Wind)      1   24.3   24.27   107.4 < 2e-16 ***
## Temp         1   13.2   13.18    58.3 1.3e-11 ***
## Residuals  101   22.8    0.23                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Anova for Nonparametric Effects
##             Npar Df Npar F   Pr(F)    
## (Intercept)                           
## s(Solar.R)        3   1.52 0.21275    
## s(Wind)           3   5.95 0.00088 ***
## Temp                                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1