3.4 ANCOVA

Avant l’apparition du modèle linéaire, une version particulière d’un mélange de régression linéaire et d’une ANOVA avec une variable indépendante quantitative et une autre variable indépendante qualitative s’appelait une ANCOVA (ANalyse de la COVariance). Un tel modèle d’ANCOVA peut naturellement également se résoudre à l’aide de la fonction lm() qui, en outre, peut faire bien plus. Nous allons maintenant ajuster un tel modèle à titre de première application concrète de tout ce que nous venons de voir sur le modèle linéaire et sur les matrices de contrastes associées.

3.4.1 Bébés à la naissance

Nous étudions la masse de nouveaux nés en fonction du poids de la mère et du fait qu’elle fume ou non. Cette analyse s’inspire de Verzani (2005). Nous avons donc ici une variable dépendante wt, la masse des bébés qui est quantitative, et deux variables indépendantes ou prédictives wt1, la masse de la mère, et smoke le fait que la mère fume ou non. Or la première de ces variables explicatives est quantitative (wt1) et l’autre (smoke) est une variable facteur à quatre niveaux (0 = la mère n’a jamais fumé, 1 = elle fume y compris pendant la grossesse, 2 = elle fumait, mais a arrêté à la grossesse, et 3 = la mère a fumé, mais a arrêté, et ce, bien avant la grossesse. Un dernier niveau 9 = inconnu encode de manière non orthodoxe les valeurs manquantes dans notre tableau de données (valeurs que nous éliminerons). De même les masses des nouveau-nés et des mères sont des unités impériales (américaines) respectivement en “onces” et en “livres”. Enfin, nous devons prendre soin de bien encoder la variable smoke comme une variable factor (ici nous ne considérons pas qu’il s’agit d’un facteur ordonné et nous voulons faire un contraste de type traitement avec comparaison à des mères qui n’ont jamais fumé). Un remaniement soigneux des données est donc nécessaire avant de pouvoir appliquer notre modèle !

SciViews::R("model", lang = "fr")
# ✓ Default language: fr
babies <- read("babies", package = "UsingR")
knitr::kable(head(babies))
id pluralty outcome date gestation sex wt parity race age ed ht wt1 drace dage ded dht dwt marital inc smoke time number
15 5 1 1411 284 1 120 1 8 27 5 62 100 8 31 5 65 110 1 1 0 0 0
20 5 1 1499 282 1 113 2 0 33 5 64 135 0 38 5 70 148 1 4 0 0 0
58 5 1 1576 279 1 128 1 0 28 2 64 115 5 32 1 99 999 1 2 1 1 1
61 5 1 1504 999 1 123 2 0 36 5 69 190 3 43 4 68 197 1 8 3 5 5
72 5 1 1425 282 1 108 1 0 23 5 67 125 0 24 5 99 999 1 1 1 1 5
100 5 1 1673 286 1 136 4 0 25 2 62 93 3 28 2 64 130 1 4 2 2 2

Ce tableau est “brut de décoffrage”. Voyez help("babies", package = "UsingR") pour de plus amples informations. Nous allons maintenant remanier tout cela correctement.

# wt = masse du bébé à la naissance en onces et 999 = valeur manquante
# wt1 = masse de la mère à la naissance en livres et 999 = valeur manquante
# smoke = 0 (non), = 1 (oui), = 2 (jusqu'à grossesse),
#       = 3 (plus depuis un certain temps) and = 9 (inconnu)
babies %>.% 
  select(., wt, wt1, smoke) %>.% # Garder seulement wt, wt1 & smoke
  filter(., wt1 < 999, wt < 999, smoke < 9) %>.% # Éliminer les valeurs manquantes
  mutate(., wt = wt * 0.02835) %>.% # Transformer le poids de livres en kg
  mutate(., wt1 = wt1 * 0.4536) %>.% # Idem de onces en kg
  mutate(., smoke = as.factor(smoke)) %->%
  Babies # Enregistrer le résultat dans Babies

knitr::kable(head(Babies))
wt wt1 smoke
3.40200 45.3600 0
3.20355 61.2360 0
3.62880 52.1640 1
3.48705 86.1840 3
3.06180 56.7000 1
3.85560 42.1848 2

Notez ici l’utilisation de l’assignation alternative %->% à la place de ->. Cette variante de l’assignation standard permet de collecter les résultats dans la nouvelle version de SciViews::R. En effet, les pipelines qui font intervenir des fonctions dites “tidy” comme select(), mutate(), filter(), group_by() |> summarise() construisent maintenant des étapes d’analyses qui sont au final combinées en un code plus performant à l’aide du package {dtplyr} en travaillant sur les data frames de classe data.table qui sont à présent utilisées.

Description des données :

skimr::skim(Babies)
Tableau 3.1: Data summary
Name Babies
Number of rows 1190
Number of columns 3
Key NULL
_______________________
Column type frequency:
factor 1
numeric 2
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
smoke 0 1 FALSE 4 0: 531, 1: 465, 3: 102, 2: 92

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
wt 0 1 3.39 0.52 1.56 3.06 3.4 3.71 4.99 ▁▃▇▅▁
wt1 0 1 58.30 9.49 39.46 51.82 56.7 62.60 113.40 ▅▇▁▁▁ 
chart(data = Babies, wt ~ wt1 %col=% smoke) +
  geom_point() +
  xlab("Masse de la mère [kg]") +
  ylab("Masse du bébé [kg]")

chart(data = Babies, wt ~ smoke) +
  geom_boxplot() +
  ylab("Masse du bébé [kg]")

chart(data = Babies, wt1 ~ smoke) +
  geom_boxplot() +
  ylab("Masse de la mère [kg]")

Visuellement, nous ne voyons pas d’effet marquant. Peut-être la condition 1 de smoke (mère qui fume pendant la grossesse) mène-t-il à des bébés moins gros, mais est-ce significatif ? Pour cela, ajustons notre modèle ANCOVA avec une matrice de traitement (choix par défaut pour une la variable factor smoke). Comme nous savons déjà utiliser lm(), c’est très simple. Cela fonctionne exactement comme avant7.

# Modèle linéaire de type (anciennement) ANCOVA
Babies_lm <- lm(data = Babies, wt ~ smoke * wt1)
summary(Babies_lm)
# 
# Call:
# lm(formula = wt ~ smoke * wt1, data = Babies)
# 
# Residuals:
#     Min      1Q  Median      3Q     Max 
# -1.9568 -0.3105  0.0133  0.3136  1.4989 
# 
# Coefficients:
#              Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  3.000663   0.128333  23.382  < 2e-16 ***
# smoke1      -0.303614   0.196930  -1.542 0.123405    
# smoke2       0.901888   0.371393   2.428 0.015314 *  
# smoke3      -0.035502   0.371379  -0.096 0.923858    
# wt1          0.008117   0.002149   3.777 0.000167 ***
# smoke1:wt1   0.001153   0.003346   0.345 0.730444    
# smoke2:wt1  -0.015340   0.006390  -2.401 0.016523 *  
# smoke3:wt1   0.001177   0.006147   0.191 0.848258    
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# Residual standard error: 0.4992 on 1182 degrees of freedom
# Multiple R-squared:  0.08248, Adjusted R-squared:  0.07705 
# F-statistic: 15.18 on 7 and 1182 DF,  p-value: < 2.2e-16
anova(Babies_lm)
# Analysis of Variance Table
# 
# Response: wt
#             Df  Sum Sq Mean Sq F value    Pr(>F)    
# smoke        3  18.659  6.2197 24.9636 1.158e-15 ***
# wt1          1   6.162  6.1621 24.7325 7.559e-07 ***
# smoke:wt1    3   1.653  0.5511  2.2117   0.08507 .  
# Residuals 1182 294.497  0.2492                      
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

L’analyse de variance (tableau du bas) montre que la masse de la mère a un effet significatif au seuil alpha de 5%, de même que la variable smoke. Par contre, il n’y a pas d’interactions entre les deux. Le fait de pouvoir mesurer des interactions entre variables qualitatives et quantitatives est ici bien évidemment un plus du modèle linéaire par rapport à ce qu’on pouvait faire avant !

Dans le tableau supérieur (coefficients), nous avons le détail du modèle linéaire complet qui tient compte de la matrice de contrastes choisie pour smoke (matrice de contrastes de type traitement qui est choisie par défaut par R pour une variable de classe factor, donc, comparaison à un modèle de base ajusté pour smoke == 0, mères n’ayant jamais fumé, premier niveau de la variable). Cette analyse est bien plus complète qu’une simple ANOVA, et nous voyons se dégager d’autres résultats lorsque les quatre conditions de smoke sont contrastées. Le résumé de l’analyse nous montre que la régression de la masse des bébés en fonction de la masse de la mère (ligne wt1 dans le tableau des coefficients), bien qu’étant significative, n’explique que 8% de la variance totale (le \(R^2\)). Suite à l’examen du nuage de points plus haut, nous n’en sommes pas surpris.

Les termes smoke1, smoke2 et smoke3 sont les contrastes appliqués par rapport au contrôle (smoke == 0). Ce sont les décalages des ordonnées à l’origine dans la régression du poids du bébé en fonction du poids de la mère, lorsque la matrice de contrastes traitement est utilisée comme ici. Par exemple, pour smoke1, l’estimation de l’ordonnée à l’origine pour la sous-population de référence smoke == 0 est de 3.00 (terme “Intercept”) -0.30 valeur de décalage estimé pour le terme smoke1, ce qui donne une ordonnée à l’origine pour la droite de la sous-population smoke == 1 de 2.70. Le même raisonnement peut être fait pour la sous-population smoke == 2(3.00 + 0.90 = 3.90) et smoke == 3 (3.00 - 0.04 = 2.96).

Les termes smoke1:wt1 à smoke3:wt1 sont les interactions du modèle, mais concrètement, cela représente la variation de la pente par rapport à la droite représentant le poids des bébés wt en fonction du poids de la mère wt1. Comme l’estimation de la pente pour la sous-population smoke == 0 de référence est de 0.0081 (terme wt1 dans le tableau), nous pouvons dire que la pente estimée pour la sous-population smoke == 1 est de 0.0081 + 0.0012 = 0.0093, pour la sous-population smoke == 2 elle est de 0.0081 - 0.0153 = -0.0072, et enfin pour la sous-population smoke == 3, elle est de 0.0081 + 0.0012 = 0.0093.

Bien. Mais maintenant, est-ce que ces décalages d’ordonnées à l’origine et de pente sont significatifs ? Nous regardons plus loin les tests t de Student associés à chaque terme (H0: terme = 0, H1: terme ≠ 0). Au seuil alpha de 5% choisi au départ avant de faire l’analyse, nous constatons que les droites pour smoke1 et smoke3 ne diffèrent pas de manière significative de celle de référence pour smoke == 0. Par contre, tant le décalage de l’ordonnée à l’origine (terme smoke2) que de la pente (terme smoke2:wt1) sont significatifs au seuil alpha de 5%. Donc, la sous-population smoke == 2 a un comportement différent en termes de variation du poids du bébé wt en fonction du poids de la mère wt1.

Pour vérifier tout cela, nous pouvons considérer le modèle de base comme une droite de régression ajustée entre wt et wt1 pour la sous-population de référence smoke == 0. Ainsi, si nous faisons :

summary(lm(data = Babies, wt ~ wt1, subset = smoke == 0))
# 
# Call:
# lm(formula = wt ~ wt1, data = Babies, subset = smoke == 0)
# 
# Residuals:
#      Min       1Q   Median       3Q      Max 
# -1.95685 -0.25825  0.01476  0.25464  1.49890 
# 
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)    
# (Intercept) 3.000663   0.123572  24.283  < 2e-16 ***
# wt1         0.008117   0.002069   3.922 9.92e-05 ***
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# Residual standard error: 0.4806 on 529 degrees of freedom
# Multiple R-squared:  0.02826, Adjusted R-squared:  0.02642 
# F-statistic: 15.38 on 1 and 529 DF,  p-value: 9.924e-05

Nous voyons en effet que les pentes et ordonnées à l’origine sont ici parfaitement identiques au modèle linéaire complet (mais pas les tests associés).

Faisons de même pour la régression entre wt et wt1 pour la sous-population smoke == 1 ?

summary(lm(data = Babies, wt ~ wt1, subset = smoke == 1))
# 
# Call:
# lm(formula = wt ~ wt1, data = Babies, subset = smoke == 1)
# 
# Residuals:
#      Min       1Q   Median       3Q      Max 
# -1.70870 -0.35089  0.01034  0.33576  1.39420 
# 
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)    
# (Intercept) 2.697048   0.153270  17.597  < 2e-16 ***
# wt1         0.009270   0.002632   3.522 0.000471 ***
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# Residual standard error: 0.5122 on 463 degrees of freedom
# Multiple R-squared:  0.02609, Adjusted R-squared:  0.02399 
# F-statistic:  12.4 on 1 and 463 DF,  p-value: 0.0004711

Nous avons une ordonnée à l’origine qui vaut 2,70 et une pente qui vaut 0.0093 comme calculé ci-dessus. L’estimation de la pente correspond donc bien à (Intercept) + smoke1 = 3,00 - 0,30 = 2,70. Cela se vérifie aussi pour les deux autres droites pour smoke2 et smoke3. De même, la pente de notre droite pour la sous-population smoke == 1 est égale à wt1 + smoke1:wt1 = 0.0081 + 0.0012 = 0.0093. Cela se vérifie aussi pour les deux autres sous-populations.

Donc, notre modèle complet ne fait rien d’autre que d’ajuster les quatre droites correspondant aux relations linéaires entre wt et wt1, mais en décompose les effets, niveau par niveau de la variable qualitative smoke en fonction de la matrice de contraste que l’on a choisie. En bonus, nous avons la possibilité de tester si chacune des composantes (tableau coefficient de summary()) ou si globalement chacune des variables (tableau obtenu avec anova()) a un effet significatif ou non dans le modèle.

Le graphique correspondant est le même que si nous avions ajusté les quatre régressions linéaires indépendamment l’une de l’autre (mais les tests et les enveloppes de confiance diffèrent).

chart(data = Babies, wt ~ wt1 %col=% smoke) +
  geom_point() +
  stat_smooth(method = "lm", formula = y ~ x) +
  xlab("Masse de la mère [kg]") +
  ylab("Masse du bébé [kg]")

chart(data = Babies, wt ~ wt1 | smoke) +
  geom_point() +
  stat_smooth(method = "lm", formula = y ~ x) +
  xlab("Masse de la mère [kg]") +
  ylab("Masse du bébé [kg]")

Sur ces deux graphiques, nous voyons très bien que les droites pour smoke == 1 et smoke == 3 sont parallèles à la droite de référence smoke == 0, mais avec un léger décalage qui n’apparaît pas significatif à 5% dans notre analyse, toutefois. Par contre, la droite relative à la sous-population smoke == 2 a une pente très différente ainsi que son ordonnée à l’origine. Et ces différences sont significatives au seuil alpha de 5% (mais ne l’auraient pas été, si nous avions choisi 1% à la place). Il faudrait investiguer plus à fond pour déterminer si cette différence est reproductible (vérifier sur un autre échantillon), et si elle n’est pas due à un biais d’échantillonnage. Enfin, il faudrait, si cela se confirme, expliquer pourquoi le poids des bébés varie différemment à partir du moment au une mère arrête de fumer au début de la grossesse…

Comme toujours, lorsqu’un effet n’est pas significatif, nous pouvons décider de simplifier le modèle. Mais attention ! Toujours considérer que les composantes sont interdépendantes. Donc, éliminer une composante du modèle peut avoir des effets parfois surprenants sur les autres.

Voyons ce que cela donne si nous éliminons les interactions, par exemple, parce que nous avons pu démontrer que le comportement différent pour smoke == 2 était dû au hasard dans une étude complémentaire qui n’aurait pas confirmé cela. Ici, nous le faisons juste dans un but pédagogique, pour monter que changer de modèle peut amener à des observations très différentes sur les mêmes données. Dans ce cas, nous ajustons des droites toutes parallèles avec uniquement un décalage de leur ordonnée à l’origine matérialisé par smoke1, smoke2 et smoke3 par rapport au modèle de référence ajusté pour la population smoke == 0 (notez l’utilisation, du signe + dans la formule, là où nous utilisions le signe * dans la modèle précédent).

# Modèle linéaire sans interactions
Babies_lm2 <- lm(data = Babies, wt ~ smoke + wt1)
summary(Babies_lm2)
# 
# Call:
# lm(formula = wt ~ smoke + wt1, data = Babies)
# 
# Residuals:
#      Min       1Q   Median       3Q      Max 
# -1.95453 -0.30780  0.01289  0.31108  1.49443 
# 
# Coefficients:
#              Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  3.030052   0.092861  32.630  < 2e-16 ***
# smoke1      -0.237938   0.031816  -7.478 1.46e-13 ***
# smoke2       0.022666   0.056508   0.401    0.688    
# smoke3       0.035486   0.054068   0.656    0.512    
# wt1          0.007617   0.001534   4.966 7.85e-07 ***
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# Residual standard error: 0.4999 on 1185 degrees of freedom
# Multiple R-squared:  0.07733, Adjusted R-squared:  0.07422 
# F-statistic: 24.83 on 4 and 1185 DF,  p-value: < 2.2e-16
anova(Babies_lm2)
# Analysis of Variance Table
# 
# Response: wt
#             Df  Sum Sq Mean Sq F value    Pr(>F)    
# smoke        3  18.659  6.2197  24.887 1.285e-15 ***
# wt1          1   6.162  6.1621  24.657 7.853e-07 ***
# Residuals 1185 296.150  0.2499                      
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Maintenant que nous avons éliminé les interactions, nous avons toujours une régression significative entre wt et wt1 (mais avec un \(R^2\) toujours très faible de 7,7%), mais maintenant, nous faisons apparaître un effet significatif du contraste avec smoke1 au seuil alpha de 5%. Et du coup, les effets des deux variables deviennent plus clairs dans notre tableau de l’ANOVA. D’une manière générale, un modèle avec interactions est toujours plus difficile à analyser du point de vue du tableau de l’ANOVA.

Le graphique correspondant est l’ajustement de droites parallèles les unes aux autres pour les quatre sous-populations en fonction de smoke. Ce graphique est difficile à réaliser. Il faut ruser, et les détails du code vont au-delà de ce cours (il n’est pas nécessaire de les comprendre à ce stade).

offsets <- c(0, -0.238, 0.0227, 0.0355)
cols <- scales::hue_pal()(4)
chart(data = Babies, wt ~ wt1) +
  geom_point(aes(col = smoke)) +
  purrr::map2(offsets, cols, function(offset, col)
    geom_smooth(method = lm, formula = y + offset ~ x, col = col)) +
  xlab("Masse de la mère [kg]") +
  ylab("Masse du bébé [kg]")

Voyons ce que donne l’analyse post hoc des comparaisons multiples (nous utilisons ici le “snippet” disponible à partir de ... -> hypothesis tests -> hypothesis tests: means -> hmanovamult : anova - multiple comparaisons [multcomp]) que nous avons déjà employée et qui reste valable ici.

summary(anovaComp. <- confint(multcomp::glht(Babies_lm2,
  linfct = multcomp::mcp(smoke = "Tukey"))))
# 
#    Simultaneous Tests for General Linear Hypotheses
# 
# Multiple Comparisons of Means: Tukey Contrasts
# 
# 
# Fit: lm(formula = wt ~ smoke + wt1, data = Babies)
# 
# Linear Hypotheses:
#            Estimate Std. Error t value Pr(>|t|)    
# 1 - 0 == 0 -0.23794    0.03182  -7.478   <1e-04 ***
# 2 - 0 == 0  0.02267    0.05651   0.401    0.977    
# 3 - 0 == 0  0.03549    0.05407   0.656    0.908    
# 2 - 1 == 0  0.26060    0.05704   4.568   <1e-04 ***
# 3 - 1 == 0  0.27342    0.05478   4.991   <1e-04 ***
# 3 - 2 == 0  0.01282    0.07199   0.178    0.998    
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# (Adjusted p values reported -- single-step method)
.oma <- par(oma = c(0, 5.1, 0, 0)); plot(anovaComp.); par(.oma); rm(.oma)

Ici, comme nous testons tous les contrastes, nous pouvons dire que la population des mères qui ont fumé pendant la grossesse smoke == 1 donne des bébés significativement moins gros au seuil alpha de 5%, et ce, en comparaison de tous les autres niveaux (mère n’ayant jamais fumé, ou ayant fumé, mais arrêté avant la grossesse, que ce soit longtemps avant ou juste avant).

Il semble évident maintenant qu’il n’est pas utile de préciser si la mère a fumé ou non avant sa grossesse. L’élément déterminant est uniquement le fait de fumer pendant la grossesse ou non. Nous pouvons le montrer également en utilisant des contrastes de Helmert, à condition de recoder smoke avec des niveaux de “gravité” croissants (“0” = n’a jamais fumé, “1” = a arrêté il y a longtemps, “2” = a arrêté juste avant la grossesse et finalement, “3” = a continué à fumer à la grossesse). Il faut donc intervertir les cas “1” et “3”. Nous pouvons utiliser recode() pour cela, mais attention, nous avons ici une variable factor, donc, ce ne sont pas des nombres, mais des chaînes de caractères (à placer entre guillemets). Une fois le recodage réalisé, il faut aussi retrier les niveaux en appelant factor(..., levels = c("0", "1", "2", "3")) sinon l’ancien ordre est conservé.

Babies %>.%
  mutate(., smoke = recode(smoke, "0" = "0", "1" = "3", "2" = "2", "3" = "1")) %>.%
  mutate(., smoke = factor(smoke, levels = c("0", "1", "2", "3"))) %->%
  Babies2

Si cela semble trop compliqué, vous pouvez aussi utiliser l’addin de réencodage dans R (QUESTIONR -> Levels Recoding or Levels Ordering). À présent que l’encodage de smoke est corrigé dans Babies2, nous pouvons modéliser à nouveau, mais cette fois-ci avec des contrastes de Helmert (notez la façon particulière de spécifier des contrastes différents de la valeur pas défaut pour une variable factor) :

Babies_lm3 <- lm(data = Babies2, wt ~ smoke + wt1,
  contrasts = list(smoke = "contr.helmert"))
summary(Babies_lm3)
# 
# Call:
# lm(formula = wt ~ smoke + wt1, data = Babies2, contrasts = list(smoke = "contr.helmert"))
# 
# Residuals:
#      Min       1Q   Median       3Q      Max 
# -1.95453 -0.30780  0.01289  0.31108  1.49443 
# 
# Coefficients:
#              Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  2.985106   0.091695  32.555  < 2e-16 ***
# smoke1       0.017743   0.027034   0.656    0.512    
# smoke2       0.001641   0.019599   0.084    0.933    
# smoke3      -0.064330   0.008540  -7.533 9.82e-14 ***
# wt1          0.007617   0.001534   4.966 7.85e-07 ***
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# Residual standard error: 0.4999 on 1185 degrees of freedom
# Multiple R-squared:  0.07733, Adjusted R-squared:  0.07422 
# F-statistic: 24.83 on 4 and 1185 DF,  p-value: < 2.2e-16
anova(Babies_lm3)
# Analysis of Variance Table
# 
# Response: wt
#             Df  Sum Sq Mean Sq F value    Pr(>F)    
# smoke        3  18.659  6.2197  24.887 1.285e-15 ***
# wt1          1   6.162  6.1621  24.657 7.853e-07 ***
# Residuals 1185 296.150  0.2499                      
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ici les valeurs estimées pour smoke1-3 sont à interpréter en fonction des contrastes utilisés, soit :

contr.helmert(4)
#   [,1] [,2] [,3]
# 1   -1   -1   -1
# 2    1   -1   -1
# 3    0    2   -1
# 4    0    0    3
  • smoke1 est le décalage de l’ordonnée à l’origine entre le modèle moyen établi avec les données smoke == 0 et smoke == 1 et celui pour smoke == 1 (non significatif au seuil alpha de 5%),
  • smoke2 est le décalage de l’ordonnée à l’origine pour smoke == 2 par rapport au modèle ajusté sur smoke == 0, smoke == 1etsmoke == 2 avec des pondérations respectives de -1, -1, et 2 (non significatif au seuil alpha de 5%),
  • smoke3 est le décalage de l’ordonnée à l’origine par rapport au modèle ajusté sur l’ensemble des autres observations, donc, avec smoke valant 0, 1, ou 2, et des pondérations respectives comme dans la dernière colonne de la matrice de contraste.. Donc, ce dernier contraste est celui qui nous intéresse, car il compare les cas où la mère n’a pas fumé pendant la grossesse avec le cas smoke == 3 où la mère a fumé pendant la grossesse, et il est significatif au seuil alpha de 5%. L’interprétation des valeurs estimées est plus complexe ici. Comparer ce résultat avec le modèle ajusté avec les contrastes de traitement par défaut avec smoke réencodé :
summary(lm(data = Babies2, wt ~ smoke + wt1))
# 
# Call:
# lm(formula = wt ~ smoke + wt1, data = Babies2)
# 
# Residuals:
#      Min       1Q   Median       3Q      Max 
# -1.95453 -0.30780  0.01289  0.31108  1.49443 
# 
# Coefficients:
#              Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  3.030052   0.092861  32.630  < 2e-16 ***
# smoke1       0.035486   0.054068   0.656    0.512    
# smoke2       0.022666   0.056508   0.401    0.688    
# smoke3      -0.237938   0.031816  -7.478 1.46e-13 ***
# wt1          0.007617   0.001534   4.966 7.85e-07 ***
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# Residual standard error: 0.4999 on 1185 degrees of freedom
# Multiple R-squared:  0.07733, Adjusted R-squared:  0.07422 
# F-statistic: 24.83 on 4 and 1185 DF,  p-value: < 2.2e-16

Les conclusions sont les mêmes, mais la valeur estimée pour smoke1, smoke2 et smoke3 diffèrent. Par exemple, dans ce dernier cas, smoke1 est double de la valeur avec les contrastes Helmert, ce qui est logique puisque la référence est ici la droite ajustée pour la sous-population smoke == 0 là où dans le modèle avec les contrastes de Helmert, le décalage est mesuré par rapport au modèle moyen (donc à “mi-chemin” entre les deux droites pour smoke == 0 et smoke == 1).

Naturellement, nous pouvons aussi considérer la variable smoke réencodée dans Babies2 comme une variable facteur ordonnée (ordered). Dans ce cas, c’est les contrastes polynomiaux qui sont utilisés :

Babies3 <- fmutate(Babies2, smoke = as.ordered(smoke))
summary(lm(data = Babies3, wt ~ smoke + wt1))
# 
# Call:
# lm(formula = wt ~ smoke + wt1, data = Babies3)
# 
# Residuals:
#      Min       1Q   Median       3Q      Max 
# -1.95453 -0.30780  0.01289  0.31108  1.49443 
# 
# Coefficients:
#              Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  2.985106   0.091695  32.555  < 2e-16 ***
# smoke.L     -0.162480   0.026780  -6.067 1.75e-09 ***
# smoke.Q     -0.148045   0.039294  -3.768 0.000173 ***
# smoke.C     -0.044605   0.048790  -0.914 0.360787    
# wt1          0.007617   0.001534   4.966 7.85e-07 ***
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# Residual standard error: 0.4999 on 1185 degrees of freedom
# Multiple R-squared:  0.07733, Adjusted R-squared:  0.07422 
# F-statistic: 24.83 on 4 and 1185 DF,  p-value: < 2.2e-16

Notez comment R est capable d’utiliser automatiquement les contrastes adéquats (polynomiaux) lorsque la variable facteur smoke est encodée en ordered. Nous voyons ici que des contrastes tenant compte d’une variation le long des successions croissantes de niveaux de “gravité” de la variable smoke sont maintenant calculés. La ligne smoke.L du tableau Coefficients indique une variation linéaire (significative au seuil alpha de 5%), smoke.Q est une variation quadratique (également significative) et enfin smoke.C est une variation cubique. Voyez la présentation des matrices de contrastes plus haut pour bien comprendre ce qui est calculé ici.

Au final, l’élément important relatif à la variable smoke est en définitive le fait de fumer pendant la grossesse ou non, pas l’histoire de la mère avant sa grossesse en matière de tabacologie ! En modélisation, nous avons toujours intérêt à choisir le modèle le plus simple. Donc ici, cela vaut le coup de simplifier smoke à une variable à deux niveaux smoke_preg qui indique uniquement si la mère fume ou non pendant la grossesse. Ensuite, nous ajustons à nouveau un modèle plus simple avec cette nouvelle variable.

Babies %>.%
  mutate(., smoke_preg = recode(smoke, "0" = "0", "1" = "1", "2" = "0", "3" = "0")) %>.%
  mutate(., smoke_preg = factor(smoke_preg, levels = c("0", "1"))) %->%
  Babies
Babies_lm4 <- lm(data = Babies, wt ~ smoke_preg + wt1)
summary(Babies_lm4)
# 
# Call:
# lm(formula = wt ~ smoke_preg + wt1, data = Babies)
# 
# Residuals:
#      Min       1Q   Median       3Q      Max 
# -1.96243 -0.30708  0.01208  0.31051  1.48662 
# 
# Coefficients:
#              Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  3.037508   0.091896  33.054  < 2e-16 ***
# smoke_preg1 -0.245797   0.029747  -8.263 3.76e-16 ***
# wt1          0.007624   0.001531   4.981 7.25e-07 ***
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# Residual standard error: 0.4996 on 1187 degrees of freedom
# Multiple R-squared:  0.07692, Adjusted R-squared:  0.07537 
# F-statistic: 49.46 on 2 and 1187 DF,  p-value: < 2.2e-16
anova(Babies_lm4)
# Analysis of Variance Table
# 
# Response: wt
#              Df  Sum Sq Mean Sq F value    Pr(>F)    
# smoke_preg    1  18.497 18.4971  74.106 < 2.2e-16 ***
# wt1           1   6.193  6.1934  24.813 7.253e-07 ***
# Residuals  1187 296.281  0.2496                      
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

À présent, tous les termes de notre modèle sont significatifs au seuil alpha de 5%. La ligne smoke_preg1 est le décalage de l’ordonnée à l’origine du poids des bébés issus de mères fumant pendant la grossesse. Il donne donc directement la perte moyenne de poids lié aux femmes qui fument ou non pendant la grossesse. La représentation graphique de ce dernier modèle est la suivante :

offsets <- c(0, -0.246)
cols <- scales::hue_pal()(2)
chart(data = Babies, wt ~ wt1) +
  geom_point(aes(col = smoke_preg)) +
  purrr::map2(offsets, cols, function(offset, col)
    geom_smooth(method = lm, formula = y + offset ~ x, col = col)) +
  xlab("Masse de la mère [kg]") +
  ylab("Masse du bébé [kg]")

Enfin, n’oublions pas que notre modèle n’est valide que si les conditions d’application sont rencontrées, en particulier, une distribution Normale des résidus et une homoscédasticité (même variance pour les résidus). Nous vérifions cela visuellement toujours avec les graphiques d’analyse des résidus. En voici les plus importants :

chart$resfitted(Babies_lm4)

chart$qqplot(Babies_lm4)

chart$scalelocation(Babies_lm4)

chart$cooksd(Babies_lm4)

Ici le comportement des résidus est sain. Des petits écarts de la normalité sur le graphique quantile-quantile s’observent peut-être, mais ce n’est pas dramatique et le modèle linéaire est robuste à ce genre de petits changements d’autant plus qu’ils apparaissent relativement symétriques en haut et et en bas de la distribution. En conclusion de cette analyse, nous pouvons dire que la masse du bébé dépend de la masse de la mère, mais assez faiblement (seulement 7,7% de la variance totale expliquée). Par contre, nous pouvons aussi dire que le fait de fumer pendant la grossesse a un effet significatif sur la réduction de la masse du bébé à la naissance (en moyenne cette réduction est de 0,246kg pour une masse moyenne à la naissance de 3,038kg, soit une réduction de 0,246 / 3,034 * 100 = 8%).

Voilà, nous venons d’analyser et d’interpréter notre premier modèle linéaire constitué d’une variable prédictive numérique et d’une variable prédictive facteur (anciennement, une ANCOVA).

3.4.2 Choix du modèle

Un modèle linéaire constitué d’une variable dépendante numérique (Y) et de deux variables explicatives dont l’une est numérique (X) et l’autre qualitative (variable “factor” F), peut s’écrire de plusieurs façons différentes. Graphiquement, il s’agit d’ajuster des droites entre X et Y qui diffèrent d’un niveau à l’autre de la variable F. Considérant que le facteur F est fixe et que les observations sont indépendantes les unes des autres, nous avons essentiellement trois cas de figure (ici, une situation fictive où F n’a que deux modalités ou niveaux nommés “A” et “B”) :

  • Les deux nuages de points sont parallèles l’un à l’autre. Les deux droites diffèrent donc seulement par un décalage de l’ordonnée à l’origine. La formule à considérer dans R est alors Y ~ X + F.

  • Les deux nuages de points s’allongent linéairement le long de droites ayant des pentes différentes, mais qui se croisent en X = 0. Les droites diffèrent donc par leurs pentes, mais pas par leurs ordonnées à l’origine. Cela s’exprime par des interactions entre X et F dans le modèle qui se notent X:F dans une formule. Nous écrirons donc notre modèle comme Y = X + X:F.

  • Les deux droites se croisent en un point quelconque différent de X = 0. Les droites diffèrent à la fois par leurs pentes et par leurs ordonnées à l’origine. Le modèle complet sera alors utilisé. Dans R nous écrirons : Y = X + F + X:F. Une forme abrégée synonyme existe aussi : Y = X * F. Ces deux formules décrivent strictement le même modèle et sont interchangeables.

En pratique, et surtout en cas de doute, vous pouvez commencer par ajuster le modèle complet. Ensuite vous analysez les différentes lignes du tableau issue de summary() appliqué à votre modèle. Vous déterminer si le test t de Student qui est réalisé sur chaque paramètre du modèle (H0: le paramètre vaut zéro, H1: le paramètre est différent de zéro) rejette l’hypothèse nulle de manière significative au seuil alpha choisi. Si pas, cela signifie que ce paramètre peut probablement être ignoré dans le modèle (c’est-à-dire, forcé à zéro dans l’équation) et vous pouvez alors passer sur une forme plus concise, soit le modèle sans interactions à droites parallèles, soit le modèle avec uniquement les interactions pour des droites qui se croisent à l’origine de l’axe des abscisses.

À vous de jouer !

Effectuez maintenant les exercices du tutoriel B03La_mod_lin (Modèle linéaire).

BioDataScience2::run("B03La_mod_lin")

Le projet suivant se fait par groupe de 2 ou de 3 (votre choix) et se poursuivra au module 4.

Réalisez en groupe le travail B03Ga_cardiovascular, partie I.

Travail en groupe de 3 pour les étudiants inscrits au cours de Science des Données Biologiques II à l’UMONS (Q1 : modélisation) à terminer avant le 2022-12-05 23:59:59.

Initiez votre projet GitHub Classroom

Voyez les explications dans le fichier README.md, partie I.

Références

Verzani, John. 2005. Using r for Introductory Statistics. Boca Raton, Florida: Chapman; Hall/CRC. https://www.crcpress.com/Using-R-for-Introductory-Statistics/Verzani/p/book/9781466590731.

  1. Pour rappel, on utilise le signe + pour indiquer un modèle sans interactions et un signe *pour spécifier un modèle complet avec interactions entre les variables.↩︎