SDD II module 3 : Modèle linéaire

Document complémentaire au module 3 du cours SDD II de 2025-2026. Distribué sous licence CC BY-NC-SA 4.0.

Veuillez vous référer au cours en ligne pour les explications et les interprétations de cette analyse.

Installer un environnement R adéquat pour reproduire cette analyse.

# Matrices de contraste par défaut
getOption("contrasts")

##         unordered           ordered 
## "contr.treatment"      "contr.poly"

# Dialecte SciViews::R avec la section dédiée à la modélisation
SciViews::R("model", lang = "FR")
# Importation des données et élimination des données manquantes pour skeleton
read("urchin_bio", package = "data.io") %>.%
  sdrop_na(., skeleton) ->
  urchin
# Graphique de base
chart(data = urchin, skeleton ~ weight %col=% origin) +
  geom_point() +
  geom_vline(xintercept = 30, col = "darkgray", lty = 2)

# Modèle linéaire avec interactions, oursins de masse supérieure ou égale à 30g
urchin_lm <- lm(data = urchin, skeleton ~ weight * origin, subset = weight >= 30)
# Graphique de ce modèle linéaire
chart(data = sfilter(urchin, weight >= 30), skeleton ~ weight %col=% origin) +
  geom_point() +
  stat_smooth(method = "lm", formula = y ~ x)

# Résumé du modèle
summary_(urchin_lm)

## 
## Call:
## lm(formula = skeleton ~ weight * origin, data = urchin, subset = weight >= 
##     30)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.15894 -0.50374 -0.05948  0.62930  2.67912 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            2.635101   0.554260   4.754 9.75e-06 ***
## weight                 0.226661   0.008803  25.747  < 2e-16 ***
## originPêcherie        -0.852076   0.899589  -0.947  0.34667    
## weight:originPêcherie  0.073522   0.018094   4.063  0.00012 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.014 on 73 degrees of freedom
## Multiple R-squared:  0.9353, Adjusted R-squared:  0.9326 
## F-statistic: 351.6 on 3 and 73 DF,  p-value: < 2.2e-16

# ANOVA du modèle linéaire
anova(urchin_lm)

## Analysis of Variance Table
## 
## Response: skeleton
##               Df Sum Sq Mean Sq F value    Pr(>F)    
## weight         1 964.34  964.34 937.606 < 2.2e-16 ***
## origin         1 103.67  103.67 100.797 2.151e-15 ***
## weight:origin  1  16.98   16.98  16.511 0.0001205 ***
## Residuals     73  75.08    1.03                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Second modèle linéaire sans décalage d'origine
urchin_lm2 <- lm(data = urchin, skeleton ~ weight + weight:origin, subset = weight >= 30)
summary_(urchin_lm2)

## 
## Call:
## lm(formula = skeleton ~ weight + weight:origin, data = urchin, 
##     subset = weight >= 30)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.25359 -0.67406 -0.05519  0.61250  2.58774 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           2.311644   0.436258   5.299 1.16e-06 ***
## weight                0.231547   0.007129  32.480  < 2e-16 ***
## weight:originPêcherie 0.057135   0.005292  10.797  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.013 on 74 degrees of freedom
## Multiple R-squared:  0.9345, Adjusted R-squared:  0.9327 
## F-statistic: 527.7 on 2 and 74 DF,  p-value: < 2.2e-16

# Anova du second modèle
anova(urchin_lm2)

## Analysis of Variance Table
## 
## Response: skeleton
##               Df Sum Sq Mean Sq F value    Pr(>F)    
## weight         1 964.34  964.34  938.91 < 2.2e-16 ***
## weight:origin  1 119.73  119.73  116.57 < 2.2e-16 ***
## Residuals     74  76.00    1.03                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Tableau formaté du résumé du modèle simplifié
summary_(urchin_lm2) |> tabularise()

Modèle linéaire
$\operatorname{skeleton} = \alpha + \beta_{1}(\operatorname{weight}) + \beta_{2}(\operatorname{weight} \times \operatorname{weight}_{\operatorname{originPêcherie}}) + \epsilon$
Terme	Valeur estimée	Ecart type	Valeur de t	Valeur de p
$\alpha$	2.3116	0.43626	5.3	1.16·10-06	***
$\beta_{1}$	0.2315	0.00713	32.5	< 2·10-16	***
$\beta_{2}$	0.0571	0.00529	10.8	< 2·10-16	***
0 <= '***' < 0.001 < '**' < 0.01 < '*' < 0.05
Etendue des résidus : [-2.254, 2.588] Ecart type des résidus : 1.013 pour 74 degrés de liberté R2 multiple : 0.9345 - R2 ajusté : 0.9327 Statistique F : 527.7 sur 2 et 74 ddl - valeur de p : < 2.22e-16

# Pour obtenir l'équation du modèle dans un document R Markdown ou Quarto:

\[\operatorname{\widehat{skeleton}} = 2.31 + 0.23(\operatorname{weight}) + 0.057(\operatorname{weight} \times \operatorname{weight}_{\operatorname{originPêcherie}})\]

# Analyse des résidus du modèle simplifié
chart$residuals(urchin_lm2)

Second exemple plus compliqué : bébés à la naissance

# Configuration de l'environnement SciViews::R en français
SciViews::R("model", lang = "fr")
# Importation des données
babies <- read("babies", package = "UsingR")
tabularise$headtail(babies[, c("wt", "wt1", "smoke")])

wt	wt1	smoke
120	100	0
113	135	0
128	115	1
123	190	3
108	125	1
...	...	...
113	100	0
128	120	0
130	150	1
125	110	0
117	129	0
Premières et dernières 5 lignes d'un total de 1236

# Remaniement des données :
# 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)
#       transformé en 0 -> "never", 3 -> "before", 2 -> "until", 1 -> "during"
#       et NA = 9 (éliminé)
babies %>.% 
  sselect(., wt, wt1, smoke) %>.% # Garder seulement wt, wt1 & smoke
  sfilter(., wt1 < 999, wt < 999, smoke < 9) %>.% # Éliminer les valeurs manquantes
  smutate(., wt = wt * 0.02835) %>.% # Transformer le poids de livres en kg
  smutate(., wt1 = wt1 * 0.4536) %>.% # Idem de onces en kg
  smutate(., smoke = recode(smoke, "0" = "never", "3" = "before",
                                   "2" = "until", "1" = "during")) %>.%
  smutate(., smoke = factor(smoke,
    levels = c("never", "before", "until", "during"))) ->
  babies2 # Enregistrer le résultat dans babies2

tabularise$headtail(babies2)

wt	wt1	smoke
3.40200	45.3600	never
3.20355	61.2360	never
3.62880	52.1640	during
3.48705	86.1840	before
3.06180	56.7000	during
...	...	...
3.20355	45.3600	never
3.62880	54.4320	never
3.68550	68.0400	during
3.54375	49.8960	never
3.31695	58.5144	never
Premières et dernières 5 lignes d'un total de 1190

# Tableaux de description des données
skimr::skim(babies2)

Data summary

Name	babies2
Number of rows	1190
Number of columns	3
_______________________
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	nev: 531, dur: 465, bef: 102, unt: 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	▅▇▁▁▁

# Graphique du poids du bébé en fonction du poids de la mère et de smoke
chart(data = babies2, wt ~ wt1 %col=% smoke) +
  geom_point() +
  xlab("wt1 : masse de la mère [kg]") +
  ylab("wt : masse du bébé [kg]")

# Boite à moustaches du poids du bébé en fonction de smoke
chart(data = babies2, wt ~ smoke) +
  geom_boxplot() +
  ylab("wt : masse du bébé [kg]")

# Boite à moustaches du poids de la mère en fonction de smoke
chart(data = babies2, wt1 ~ smoke) +
  geom_boxplot() +
  ylab("wt1 : masse de la mère [kg]")

# Modèle linéaire de type (anciennement) ANCOVA
babies2_lm <- lm(data = babies2, wt ~ smoke * wt1)
summary_(babies2_lm) |> tabularise()

Modèle linéaire
$\operatorname{wt} = \alpha + \beta_{1}(\operatorname{smoke}_{\operatorname{before}}) + \beta_{2}(\operatorname{smoke}_{\operatorname{until}}) + \beta_{3}(\operatorname{smoke}_{\operatorname{during}}) + \beta_{4}(\operatorname{wt1}) + \beta_{5}(\operatorname{smoke}_{\operatorname{before}} \times \operatorname{wt1}) + \beta_{6}(\operatorname{smoke}_{\operatorname{until}} \times \operatorname{wt1}) + \beta_{7}(\operatorname{smoke}_{\operatorname{during}} \times \operatorname{wt1}) + \epsilon$
Terme	Valeur estimée	Ecart type	Valeur de t	Valeur de p
$\alpha$	3.00066	0.12833	23.3818	< 2·10-16	***
$\beta_{1}$	-0.03550	0.37138	-0.0956	9.24·10-01
$\beta_{2}$	0.90189	0.37139	2.4284	1.53·10-02	*
$\beta_{3}$	-0.30361	0.19693	-1.5417	1.23·10-01
$\beta_{4}$	0.00812	0.00215	3.7768	1.67·10-04	***
$\beta_{5}$	0.00118	0.00615	0.1914	8.48·10-01
$\beta_{6}$	-0.01534	0.00639	-2.4006	1.65·10-02	*
$\beta_{7}$	0.00115	0.00335	0.3446	7.30·10-01
0 <= '***' < 0.001 < '**' < 0.01 < '*' < 0.05
Etendue des résidus : [-1.957, 1.499] Ecart type des résidus : 0.4992 pour 1182 degrés de liberté R2 multiple : 0.08248 - R2 ajusté : 0.07705 Statistique F : 15.18 sur 7 et 1182 ddl - valeur de p : < 2.22e-16

anova(babies2_lm) |> tabularise()

Terme	Ddl	Somme des carrés	Carrés moyens	Valeur de Fobs.	Valeur de p
smoke	3	18.66	6.220	24.96	1.16·10-15	***
wt1	1	6.16	6.162	24.73	7.56·10-07	***
smoke:wt1	3	1.65	0.551	2.21	8.51·10-02	.
Résidus	1 182	294.50	0.249
0 <= '***' < 0.001 < '**' < 0.01 < '*' < 0.05

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

# Modèle linéaire sans interactions (+ au lieu de * dans la formule)
babies2_lm2 <- lm(data = babies2, wt ~ smoke + wt1)
summary_(babies2_lm2) |> tabularise()

$\operatorname{wt} = \alpha + \beta_{1}(\operatorname{smoke}_{\operatorname{before}}) + \beta_{2}(\operatorname{smoke}_{\operatorname{until}}) + \beta_{3}(\operatorname{smoke}_{\operatorname{during}}) + \beta_{4}(\operatorname{wt1}) + \epsilon$
Terme	Valeur estimée	Ecart type	Valeur de t	Valeur de p
$\alpha$	3.03005	0.09286	32.630	< 2·10-16	***
$\beta_{1}$	0.03549	0.05407	0.656	5.12·10-01
$\beta_{2}$	0.02267	0.05651	0.401	6.88·10-01
$\beta_{3}$	-0.23794	0.03182	-7.478	1.46·10-13	***
$\beta_{4}$	0.00762	0.00153	4.966	7.85·10-07	***
0 <= '***' < 0.001 < '**' < 0.01 < '*' < 0.05
Etendue des résidus : [-1.955, 1.494] Ecart type des résidus : 0.4999 pour 1185 degrés de liberté R2 multiple : 0.07733 - R2 ajusté : 0.07422 Statistique F : 24.83 sur 4 et 1185 ddl - valeur de p : < 2.22e-16

anova(babies2_lm2) |> tabularise()

Terme	Ddl	Somme des carrés	Carrés moyens	Valeur de Fobs.	Valeur de p
smoke	3	18.66	6.22	24.9	1.28·10-15	***
wt1	1	6.16	6.16	24.7	7.85·10-07	***
Résidus	1 185	296.15	0.25
0 <= '***' < 0.001 < '**' < 0.01 < '*' < 0.05

# Graphique du dernier modèle (difficile à réaliser)
offsets <- c(0, 0.03549, 0.02267, -0.23794)
cols <- scales::hue_pal()(4)
chart(data = babies2, 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("wt1 : masse de la mère [kg]") +
  ylab("wt : masse du bébé [kg]")

# Analyse post-hoc du dernier modèle
babies2_lm2_compa <- multcomp::glht(babies2_lm2,
  linfct = multcomp::mcp(smoke = "Tukey")) |> confint()
summary_(babies2_lm2_compa)

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = wt ~ smoke + wt1, data = babies2)
## 
## Linear Hypotheses:
##                      Estimate Std. Error t value Pr(>|t|)    
## before - never == 0   0.03549    0.05407   0.656    0.908    
## until - never == 0    0.02267    0.05651   0.401    0.977    
## during - never == 0  -0.23794    0.03182  -7.478   <1e-04 ***
## until - before == 0  -0.01282    0.07199  -0.178    0.998    
## during - before == 0 -0.27342    0.05478  -4.991   <1e-04 ***
## during - until == 0  -0.26060    0.05704  -4.568   <1e-04 ***
## ---
## 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(babies2_lm2_compa); par(.oma); rm(.oma)

# Nouveau modèle utilisant les contrastes de Helmert
babies2_lm3 <- lm(data = babies2, wt ~ smoke + wt1,
  contrasts = list(smoke = "contr.helmert"))
summary_(babies2_lm3) |> tabularise()

$\operatorname{wt} = \alpha + \beta_{1}(\operatorname{smoke}_{\operatorname{1}}) + \beta_{2}(\operatorname{smoke}_{\operatorname{2}}) + \beta_{3}(\operatorname{smoke}_{\operatorname{3}}) + \beta_{4}(\operatorname{wt1}) + \epsilon$
Terme	Valeur estimée	Ecart type	Valeur de t	Valeur de p
$\alpha$	2.98511	0.09170	32.5547	< 2·10-16	***
$\beta_{1}$	0.01774	0.02703	0.6563	5.12·10-01
$\beta_{2}$	0.00164	0.01960	0.0837	9.33·10-01
$\beta_{3}$	-0.06433	0.00854	-7.5328	9.82·10-14	***
$\beta_{4}$	0.00762	0.00153	4.9656	7.85·10-07	***
0 <= '***' < 0.001 < '**' < 0.01 < '*' < 0.05
Etendue des résidus : [-1.955, 1.494] Ecart type des résidus : 0.4999 pour 1185 degrés de liberté R2 multiple : 0.07733 - R2 ajusté : 0.07422 Statistique F : 24.83 sur 4 et 1185 ddl - valeur de p : < 2.22e-16

anova(babies2_lm3) |> tabularise()

Terme	Ddl	Somme des carrés	Carrés moyens	Valeur de Fobs.	Valeur de p
smoke	3	18.66	6.22	24.9	1.28·10-15	***
wt1	1	6.16	6.16	24.7	7.85·10-07	***
Résidus	1 185	296.15	0.25
0 <= '***' < 0.001 < '**' < 0.01 < '*' < 0.05

# Matrice de Helmert correspondant à notre modèle
helm <- contr.helmert(4)
rownames(helm) <- c('smoke == "never"', 'smoke == "before"',
                    'smoke == "until"', 'smoke == "during"')
colnames(helm) <- c('smoke1', 'smoke2', 'smoke3')
helm

##                   smoke1 smoke2 smoke3
## smoke == "never"      -1     -1     -1
## smoke == "before"      1     -1     -1
## smoke == "until"       0      2     -1
## smoke == "during"      0      0      3

# Encore un autre modèle, en considérant smoke comme variable ordonnée
smutate(babies2, smoke = as.ordered(smoke)) |>
  lm(data = _, wt ~ smoke + wt1) |>
  summary_() |>
  tabularise()

$\operatorname{wt} = \alpha + \beta_{1}(\operatorname{smoke}_{\operatorname{.L}}) + \beta_{2}(\operatorname{smoke}_{\operatorname{.Q}}) + \beta_{3}(\operatorname{smoke}_{\operatorname{.C}}) + \beta_{4}(\operatorname{wt1}) + \epsilon$
Terme	Valeur estimée	Ecart type	Valeur de t	Valeur de p
$\alpha$	2.98511	0.09170	32.555	< 2·10-16	***
$\beta_{1}$	-0.16248	0.02678	-6.067	1.75·10-09	***
$\beta_{2}$	-0.14804	0.03929	-3.768	1.73·10-04	***
$\beta_{3}$	-0.04460	0.04879	-0.914	3.61·10-01
$\beta_{4}$	0.00762	0.00153	4.966	7.85·10-07	***
0 <= '***' < 0.001 < '**' < 0.01 < '*' < 0.05
Etendue des résidus : [-1.955, 1.494] Ecart type des résidus : 0.4999 pour 1185 degrés de liberté R2 multiple : 0.07733 - R2 ajusté : 0.07422 Statistique F : 24.83 sur 4 et 1185 ddl - valeur de p : < 2.22e-16

# Recodage de smoke en smokepreg : contraster fumer pendant la grossesse ou non
babies2 %>.%
  smutate(., smokepreg = recode(smoke, "never" = "0", "before" = "0",
                                       "until" = "0", "during" = "1")) %>.%
  smutate(., smokepreg = factor(smokepreg, levels = c("0", "1"))) ->
  babies2
# Modèle utilisant smokepreg
babies2_lm4 <- lm(data = babies2, wt ~ smokepreg + wt1)
summary_(babies2_lm4) |> tabularise()

$\operatorname{wt} = \alpha + \beta_{1}(\operatorname{smokepreg}_{\operatorname{1}}) + \beta_{2}(\operatorname{wt1}) + \epsilon$
Terme	Valeur estimée	Ecart type	Valeur de t	Valeur de p
$\alpha$	3.03751	0.09190	33.05	< 2·10-16	***
$\beta_{1}$	-0.24580	0.02975	-8.26	3.76·10-16	***
$\beta_{2}$	0.00762	0.00153	4.98	7.25·10-07	***
0 <= '***' < 0.001 < '**' < 0.01 < '*' < 0.05
Etendue des résidus : [-1.962, 1.487] Ecart type des résidus : 0.4996 pour 1187 degrés de liberté R2 multiple : 0.07692 - R2 ajusté : 0.07537 Statistique F : 49.46 sur 2 et 1187 ddl - valeur de p : < 2.22e-16

anova(babies2_lm4) |> tabularise()

Terme	Ddl	Somme des carrés	Carrés moyens	Valeur de Fobs.	Valeur de p
smokepreg	1	18.50	18.50	74.1	< 2·10-16	***
wt1	1	6.19	6.19	24.8	7.25·10-07	***
Résidus	1 187	296.28	0.25
0 <= '***' < 0.001 < '**' < 0.01 < '*' < 0.05

# Graphique du modèle linéaire avec smokepreg
offsets <- c(0, -0.2458)
cols <- scales::hue_pal()(2)
chart(data = babies2, wt ~ wt1) +
  geom_point(aes(col = smokepreg)) +
  purrr::map2(offsets, cols, function(offset, col)
    geom_smooth(method = lm, formula = y + offset ~ x, col = col)) +
  xlab("wt1 : masse de la mère [kg]") +
  ylab("wt : masse du bébé [kg]")

# Graphique principal d'analyse des résidus
chart$resfitted(babies2_lm4)

# Graphique quantile-quantile d'analyse des résidus
chart$qqplot(babies2_lm4)

# Graphique d'analyse des résidus pour vérifier l'homoscédasticité
chart$scalelocation(babies2_lm4)

# Graphique de la distance de Cook des résidus
chart$cooksd(babies2_lm4)