-
Notifications
You must be signed in to change notification settings - Fork 0
/
08-GLM.qmd
5107 lines (4300 loc) · 302 KB
/
08-GLM.qmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# Régressions linéaires généralisées (GLM) {#sec-chap08}
Dans ce chapitre, nous présentons les modèles linéaires généralisés plus communément appelés GLM (*generalized linear models* en anglais). Il s’agit d’une extension directe du modèle de régression linéaire multiple (LM) basé sur la méthode des moindres carrés ordinaires, décrite dans le chapitre précédent. Pour aborder cette section sereinement, il est important d’avoir bien compris le concept de distribution présenté dans la [section @sec-024]. À la fin de cette section, vous serez en mesure de :
* Comprendre la distinction entre un modèle LM classique et un GLM.
* Identifier les composantes d’un GLM.
* Interpréter les résultats d’un GLM.
* Effectuer les diagnostics d’un GLM.
::: bloc_package
::: bloc_package-header
::: bloc_package-icon
:::
**Liste des *packages* utilisés dans ce chapitre**
:::
::: bloc_package-body
Dans ce chapitre, nous utilisons principalement les *packages* suivants :
* Pour créer des graphiques :
- `ggplot2`, le seul, l'unique!
- `ggpubr` pour combiner des graphiques et réaliser des diagrammes.
* Pour ajuster des modèles GLM :
- `VGAM` et `gamlss` offrent tous les deux un très large choix de distributions et de fonctions de diagnostic, mais nécessitent souvent un peu plus de code.
- `mgcv` offre moins de distributions que les deux précédents, mais est plus simple d'utilisation.
* Pour analyser des modèles GLM :
- `car` essentiellement pour la fonction `vif`.
- `DHARMa` pour le diagnostic des résidus simulés.
- `ROCR` et `caret` pour l'analyse de la qualité d'ajustement de modèles pour des variables qualitatives.
- `AER` pour des tests de sur-dispersion.
- `fitdistrplus` pour ajuster des distributions à des données.
- `LaplacesDemon` pour manipuler certaines distributions.
- `sandwich` pour générer des erreurs standards robustes pour le modèle GLM logistique binomial.
:::
:::
## Qu'est qu'un modèle GLM? {#sec-081}
Nous avons vu qu’une régression linéaire multiple (LM) ne peut être appliquée que si la variable dépendante analysée est continue et si elle est normalement distribuée, une fois les variables indépendantes contrôlées. Il s’agit d’une limite très importante puisqu’elle ne peut être utilisée pour modéliser et prédire des variables binaires, multinomiales, de comptage, ordinales ou plus simplement des données anormalement distribuées. Une seconde limite importante des LM est que l’influence des variables indépendantes sur la variable dépendante ne peut être que linéaire. L’augmentation d’une unité de *X* conduit à une augmentation (ou diminution) de $\beta$ (coefficient de régression) unités de *Y*, ce qui n’est pas toujours représentatif des phénomènes étudiés. Afin de dépasser ces contraintes, @GLMnelder ont proposé une extension des modèles LM, soit les modèles linéaires généralisés (GLM).
### Formulation d'un GLM {#sec-0811}
Puisqu’un modèle GLM est une extension des modèles LM, il est possible de traduire un modèle LM sous forme d’un GLM. Nous utilisons ce point de départ pour détailler la morphologie d’un GLM. Nous avons vu dans la section précédente qu’un modèle LM est formulé de la façon suivante (notation matricielle) :
$$
Y = \beta_0 + X\beta + \epsilon
$$ {#eq-regmultiple5B}
Avec $\beta_0$ la constante (*intercept* en anglais) et $\beta$ un vecteur de coefficients de régression pour les *k* variables indépendantes (*X*).
D'après cette formule, nous modélisons la variable *Y* avec une équation de régression linéaire et un terme d’erreur que nous estimons être normalement distribué. Nous pouvons reformuler ce simple LM sous forme d’un GLM avec l’écriture suivante :
$$
\begin{aligned}
&Y \sim Normal(\mu,\sigma)\\
&g(\mu) = \beta_0 + \beta X\\
&g(x) = x
\end{aligned}
$$ {#eq-glm1}
Pas de panique! Cette écriture se lit comme suit : la variable *Y* est issue d’une distribution normale $(Y \sim Normal)$ avec deux paramètres : $\mu$ (sa moyenne) et $\sigma$ (son écart-type). $\mu$ varie en fonction d’une équation de régression linéaire ($\beta_0 + \beta X$) transformée par une fonction de lien *g* (détaillée plus loin). Dans ce cas précis, la fonction de lien est appelée fonction identitaire puisqu’elle n’applique aucune transformation ($g(x) = x$). Notez ici que le second paramètre de la distribution normale $\sigma$ (paramètre de dispersion) a une valeur fixe et ne dépend donc pas des variables indépendantes à la différence de $\mu$. Dans ce modèle spécifiquement, les paramètres à estimer sont $\sigma$, $\beta_0$ et $\beta$.
Notez que dans la notation traditionnelle, la fonction de lien est appliquée au paramètre modélisé. Il est possible de renverser cette notation en utilisant la réciproque ($g'$) de la fonction de lien ($g$) :
$$
g(\mu) = \beta_0 + \beta X \Longleftrightarrow \mu = g'(\beta_0 + \beta X)
\text{ si : }g'(g(x)) = x
$$ {#eq-glm2}
Dans un modèle GLM, la distribution attendue de la variable *Y* est déclarée de façon explicite ainsi que la façon dont nos variables indépendantes conditionnent cette distribution. Ici, c’est la moyenne ($\mu$) de la distribution qui est modélisée, nous nous intéressons ainsi au changement moyen de *Y* provoqué par les variables *X*.
Avec cet exemple, nous voyons les deux composantes supplémentaires d’un modèle GLM :
* La distribution supposée de la variable *Y* conditionnée par les variables *X* (ici, la distribution normale).
* Une fonction de lien associant l’équation de régression formée par les variables indépendantes et un paramètre de la distribution retenue (ici, la fonction identitaire et le paramètre $\mu$).
Notez également que l’estimation des paramètres d’un modèle GLM (ici, $\beta_0$, $\beta X$ et $\sigma$) ne se fait plus avec la méthode des moindres carrés ordinaires utilisée pour les modèles LM. À la place, la méthode par maximum de vraisemblance (*maximum likelihood*) est la plus souvent utilisée, mais certains *packages* utilisent également la méthode des moments (*method of moments*). Ces deux méthodes nécessitent des échantillons plus grands que la méthode des moindres carrés. Dans le cas spécifique d'un modèle GLM utilisant une distribution normale, la méthode des moindres carrés et la méthode par maximum de vraisemblance produisent les mêmes résultats.
### Autres distributions et rôle de la fonction de lien {#sec-0812}
À première vue, il est possible de se demander pourquoi ajouter ces deux éléments puisqu’ils ne font que complexifier le modèle. Pour mieux saisir la pertinence des GLM, prenons un exemple appliqué au cas d’une variable binaire. Admettons que nous souhaitons modéliser / prédire la probabilité qu’une personne à vélo décède lors d’une collision avec un véhicule motorisé. Notre variable dépendante est donc binaire (0 = survie, 1 = décès) et nous souhaitons la prédire avec trois variables continues que sont : la vitesse de déplacement du ou de la cycliste ($x_1$), la vitesse de déplacement du véhicule ($x_2$) et la masse du véhicule ($x_3$). Puisque la variable *Y* n’est pas continue, il serait absurde de supposer qu’elle est issue d’une distribution normale. Il serait plus logique de penser qu’elle provient d’une distribution de Bernoulli (pour rappel, une distribution de Bernoulli permet de modéliser un phénomène ayant deux issues possibles comme un lancer de pièce de monnaie ([section @sec-024]). Plus spécifiquement, nous pourrions formuler l’hypothèse que nos trois variables $x_1$, $x_2$ et $x_3$ influencent le paramètre *p* (la probabilité d’occurrence de l’évènement) d’une distribution de Bernoulli. À partir de ces premières hypothèses, nous pouvons écrire le modèle suivant :
$$
\begin{aligned}
&Y \sim Bernoulli(p)\\
&g(p) = \beta_0 + \beta X\\
&g(x) = x
\end{aligned}
$$ {#eq-glm3}
Toutefois, le résultat n’est pas entièrement satisfaisant. En effet, *p* est une probabilité et, par nature, ce paramètre doit être compris entre 0 et 1 (entre 0 et 100 % de « chances de décès », ni plus ni moins). L’équation de régression que nous utilisons actuellement peut produire des résultats compris entre $-\infty$ et $+\infty$ pour *p* puisque rien ne contraint la somme $\beta_0+ \beta_1x_1+\beta_2x_2+ \beta_3x_3$ à être comprise entre 0 et 1. Il est possible de visualiser le problème soulevé par cette situation avec les figures suivantes. Admettons que nous avons observé une variable *Y* binaire et que nous savons qu’elle est influencée par une variable *X* qui, plus elle augmente, plus les chances que *Y* soit 1 augmentent (@fig-linearbinom).
```{r}
#| label: fig-linearbinom
#| fig-align: center
#| fig-cap: Exemple de données issues d'une distribution de Bernoulli
#| echo: false
#| message: false
#| warning: false
#| out-width: "70%"
source("code_complementaire/JG_helper.R")
library(ggplot2)
x1 <- rnorm(150)
z <- 1 + rnorm(150, mean = 2, sd = 0.5)*x1
pr <- 1/(1+exp(-z))
y <- ifelse(pr > 0.5,1,0)
df <- data.frame(x1 = x1,
y = y)
ggplot(data = df)+
geom_point(aes(x = x1, y = y), size = 0.5) +
labs(x = "X", y = "Y") +
ylim(c(-0.25,1.25))
```
Si nous utilisons l’équation de régression actuelle, cela revient à trouver la droite la mieux ajustée passant dans ce nuage de points (@fig-linearbinom2).
```{r}
#| label: fig-linearbinom2
#| fig-align: center
#| echo: false
#| fig-cap: Ajustement d'une droite de régression aux données issues d'une distribution de Bernoulli
#| message: false
#| warning: false
#| out-width: "70%"
tofr <- function(float){
return(gsub("." , ",", as.character(float), fixed = TRUE, useBytes = TRUE))
}
ggplot(data = df)+
geom_point(aes(x = x1, y = y), size = 0.5) +
geom_smooth(aes(x = x1, y = y), method='lm', formula= y~x, se = FALSE) +
labs(x = "X", y = "Y") +
ylim(c(-0.25,1.25))
```
Ce modèle semble bien cerner l’influence positive de *X* sur *Y*, mais la droite est au final très éloignée de chaque point, indiquant un faible ajustement du modèle. De plus, la droite prédit des probabilités négatives lorsque *X* est inférieur à −2,5 et des probabilités supérieures à 1 quand *X* est supérieur à 1. Elle est donc loin de bien représenter les données.
C’est ici qu’intervient la fonction de lien. La fonction identitaire que nous avons utilisée jusqu'ici n’est pas satisfaisante, nous devons la remplacer par une fonction qui conditionnera la somme $\beta_0+ \beta_1x_1+\beta_2x_2+ \beta_3x_3$ pour donner un résultat entre 0 et 1. Une candidate toute désignée est la fonction *sigmoidale*, plus souvent appelée la fonction *logistique*!
$$
\begin{aligned}
&Y \sim Bernoulli(p)\\
&S(p) = \beta_0 + \beta X\\
&S(x) = \frac{e^{x}}{e^x+1}
\end{aligned}
$$ {#eq-glm4}
La fonction logistique prend la forme d’un *S*. Plus la valeur entrée dans la fonction est grande, plus le résultat produit par la fonction est proche de 1 et inversement. Si nous reprenons l’exemple précédent, nous obtenons le modèle illustré à la @fig-linearbinom3.
```{r}
#| label: fig-linearbinom3
#| fig-align: center
#| fig-cap: Utilisation de la fonction de lien logistique
#| echo: false
#| message: false
#| warning: false
#| out-width: "70%"
ggplot(data = df)+
geom_point(aes(x = x1, y = y), size = 0.5) +
geom_smooth(aes(x = x1, y = y), method='lm', formula= y~x, se = FALSE) +
geom_smooth(aes(x = x1, y = y), method='glm',
method.args = list(family = "binomial"),
se = FALSE, color = "red") +
labs(x = "X", y = "Y") +
ylim(c(-0.25,1.25))
```
Une fois cette fonction insérée dans le modèle, nous constatons qu’une augmentation de la somme $\beta_0+ \beta_1x_1+\beta_2x_2+ \beta_3x_3$ conduit à une augmentation de la probabilité *p* et inversement, et que cet effet est non linéaire. Nous avons donc maintenant un GLM permettant de prédire la probabilité d’un décès lors d’un accident en combinant une distribution et une fonction de lien adéquates.
### Conditions d'application {#sec-0813}
La famille des GLM englobe de (très) nombreux modèles du fait de la diversité de distributions existantes et des fonctions de liens utilisables. Cependant, certaines combinaisons sont plus souvent utilisées que d’autres. Nous présentons donc dans les prochaines sections les modèles GLM les plus communs. Les conditions d’application varient d’un modèle à l’autre, il existe cependant quelques conditions d’application communes à tous ces modèles :
* l’indépendance des observations (et donc des erreurs);
* l’absence de valeurs aberrantes / fortement influentes;
* l’absence de multicolinéarité excessive entre les variables indépendantes.
Ces trois conditions sont également valables pour les modèles LM tel qu'abordé dans le [chapitre @sec-chap07]. La distance de *Cook* peut ainsi être utilisée pour détecter les potentielles valeurs aberrantes et le facteur d’inflation de la variance (*VIF*) pour détecter la multicolinéarité.
Les conditions d’application particulières sont détaillées dans les sections dédiées à chaque modèle.
### Résidus et déviance {#sec-0814}
Dans la section sur la régression linéaire simple, nous avons présenté la notion de résidu, soit l’écart entre la valeur observée (réelle) de *Y* et la valeur prédite par le modèle. Pour un modèle GLM, ces résidus traditionnels (aussi appelés résidus naturels) ne sont pas très informatifs si la variable à modéliser est binaire, multinomiale ou même de comptage. Lorsque l’on travaille avec des GLM, nous préférons utiliser trois autres formes de résidus, soit les résidus de Pearson, les résidus de déviance et les résidus simulés.
**Les résidus de Pearson** sont une forme ajustée des résidus classiques, obtenus par la division des résidus naturels par la racine carrée de la variance modélisée. Leur formule varie donc d’un modèle à l’autre puisque l’expression de la variance change en fonction de la distribution du modèle. Pour un modèle GLM gaussien, elle s'écrit :
$$
r_i = \frac{y_i - \mu_i}{\sigma}
$$ {#eq-glm5}
Pour un modèle GLM de Bernoulli, elle s'écrit :
$$
r_i = \frac{y_i - p_i}{\sqrt{p_i(1-p_i)}}
$$ {#eq-glm6}
avec $\mu_i$ et $p_i$ les prédictions du modèle pour l'observation *i*.
**Les résidus de déviance** sont basés sur le concept de *likelihood* présenté dans la [section @sec-02adjdistrib]. Pour rappel, le *likelihood*, ou la vraisemblance d’un modèle, correspond à la probabilité conjointe d’avoir observé les données *Y* selon le modèle étudié. Pour des raisons mathématiques (voir [section @sec-02adjdistrib]), le *log likelihood* est plus souvent calculé. Plus cette valeur est forte, moins le modèle se trompe. Cette interprétation est donc inverse à celle des résidus classiques, c’est pourquoi le *log likelihood* est généralement multiplié par −2 pour retrouver une interprétation intuitive. Ainsi, pour chaque observation *i*, nous pouvons calculer :
$$
d_i = \mbox{-2} \times log(P(y_i|M_e))
$$ {#eq-glm7}
avec $d_i$ le résidu de déviance et $P(y_i|M_e)$ la probabilité d’avoir observé la valeur $y_i$ selon le modèle étudié ($M_e$).
La somme de tous ces résidus est appelée la déviance totale du modèle.
$$
D(M_e) = \sum_{i=1}^n \mbox{-2} \times log(P(y_i|M_e))
$$ {#eq-glm8}
Il s’agit donc d’une quantité représentant à quel point le modèle est erroné vis-à-vis des données. Notez qu’en tant que telle, la déviance n’a pas d’interprétation directe en revanche, elle est utilisée pour calculer des mesures d’ajustement des modèles GLM.
**Les résidus simulés** sont une avancée récente dans le monde des GLM, ils fournissent une définition et une interprétation harmonisée des résidus pour l’ensemble des modèles GLM. Dans la section sur les LM ([section @sec-0722]), nous avons vu comment interpréter les graphiques des résidus pour détecter d’éventuels problèmes dans le modèle. Cependant, cette technique est bien plus compliquée à mettre en œuvre pour les GLM puisque la forme attendue des résidus varie en fonction de la distribution choisie pour modéliser *Y*. La façon la plus efficace de procéder est d’interpréter les graphiques des résidus simulés qui ont la particularité d’être **identiquement distribués, quel que soit le modèle GLM construit**. Ces résidus simulés sont compris entre 0 et 1 et sont calculés de la manière suivante :
* À partir du modèle GLM construit, simuler *S* fois (généralement 1 000) une variable *Y’* avec autant d’observation (*n*) que *Y*. Cette variable simulée est une combinaison de la prédiction du modèle (coefficient et variables indépendantes) et de sa dispersion (variance). Ces simulations représentent des variations vraisemblables de la variable *Y* si le modèle est correctement spécifié. En d’autres termes, si le modèle représente bien le phénomène à l’origine de la variable *Y*, alors les simulations *Y’* issues du modèle devraient être proches de la variable *Y* originale. Pour une explication plus détaillée de ce que signifie simuler des données à partir d’un modèle, référez-vous au *bloc attention* intitulé *Distinction entre simulation et prédiction* dans la [section @sec-08152].
* Pour chaque observation, nous obtenons ainsi *S* valeurs formant une distribution $Ds_i$, soit les valeurs simulées par le modèle pour cette observation.
* Pour chacune de ces distributions, nous calculons la probabilité cumulative d’observer la vraie valeur $Y_i$ d'après la distribution $Ds_i$. Cette valeur est comprise entre 0 (toutes les valeurs simulées sont plus grandes que $Y_i$) et 1 (toutes les valeurs simulées sont inférieures à $Y_i$).
Si le modèle est correctement spécifié, le résultat attendu est que la distribution de ces résidus est uniforme. En effet, il y a autant de chances que les simulations produisent des résultats supérieurs ou inférieurs à $Y_i$ si le modèle représente bien le phénomène [@RandomizedResid;@gelman2006data]. Si la distribution des résidus ne suit pas une loi uniforme, cela signifie que le modèle échoue à reproduire le phénomène à l’origine de *Y*, ce qui doit nous alerter sur sa pertinence.
### Vérification l’ajustement {#sec-0815}
Il existe trois façons de vérifier l’ajustement d’un modèle GLM :
* utiliser des mesures d’ajustement (AIC, pseudo-R^2^, déviance expliquée, etc.);
* comparer les distributions de la variable originale et celle des prédictions;
* comparer les prédictions du modèle avec les valeurs originales.
Notez d’emblée que vérifier la qualité d’ajustement d’un modèle (ajustement aux données originales) ne revient pas à vérifier la validité d’un modèle (respect des conditions d’application). Cependant, ces deux éléments sont généralement liés, car un modèle mal ajusté a peu de chances d’être valide et inversement.
#### Mesures d'ajustement {#sec-08151}
Les mesures d’ajustement sont des indicateurs plus ou moins arbitraires dont le principal intérêt est de faciliter la comparaison entre plusieurs modèles similaires. Il est nécessaire de les reporter, car dans certains cas, ils peuvent indiquer que des modèles sont très mal ajustés.
##### Déviance expliquée {#sec-081511}
Rappelons que la déviance d’un modèle est une quantité représentant à quel point le modèle est erroné. L’objectif de l’indicateur de la déviance expliquée est d’estimer le pourcentage de la déviance maximale observable dans les données que le modèle est parvenu à expliquer. La déviance maximale observable dans les données est obtenue en utilisant la déviance totale du modèle nul (notée $M_n$, soit un modèle dans lequel aucune variable indépendante n’est ajoutée et ne comportant qu’une constante). Cette déviance est maximale puisqu’aucune variable indépendante n’est présente dans le modèle. Nous calculons ensuite le pourcentage de cette déviance totale qui a été contrôlée par le modèle étudié ($M_e$).
$$
\mbox{déviance expliquée} = \frac{D(M_n) - D(M_e)}{D(M_n)} = 1- \frac{D(M_e)}{D(M_n)}
$$ {#eq-glm9}
Il s’agit donc d’un simple calcul de pourcentage entre la déviance maximale ($D(M_n)$) et la déviance expliquée par le modèle étudié ($D(M_n )-D(M_e)$). Cet indicateur est compris entre 0 et 1 : plus il est petit, plus la capacité de prédiction du modèle est faible. Attention, cet indicateur ne tient pas compte de la complexité du modèle. Ajouter une variable indépendante supplémentaire ne fait qu’augmenter la déviance expliquée, ce qui ne signifie pas que la complexification du modèle soit justifiée (**voir l'encadré sur le principe de parcimonie**, [section @sec-0732]).
##### Pseudo-R^2^ {#sec-081512}
Le R^2^ est une mesure d’ajustement représentant la part de la variance expliquée dans un modèle linéaire classique. Cette mesure n’est pas directement transposable au cas des GLM puisqu’ils peuvent être appliqués à des variables non continues et anormalement distribuées. Toutefois, il existe des mesures semblables appelées pseudo-R^2^, remplissant un rôle similaire. Notez cependant qu’ils ne peuvent pas être interprétés comme le R^2^ classique (d'une régression linéaire multiple) : **ils ne représentent pas la part de la variance expliquée**. Ils sont compris dans l’intervalle 0 et 1; plus leurs valeurs s’approchent de 1, plus le modèle est ajusté.
```{r}
#| label: tbl-pseudor2
#| tbl-cap: Principaux pseudos-$R^2$
#| echo: false
#| message: false
#| warning: false
df <- data.frame(
Nom = c("McFadden", "McFadden ajusté", "Efron", "Cox & Snell", "Nagelkerke"),
formule = c("$1-\\frac{loglike(M_e)}{loglike(M_n)}$",
"$1-\\frac{loglike(M_e)-K}{loglike(M_n)}$",
"$1-\\frac{\\sum_{i=1}^n(y_i-\\hat{y}_i)^2}{\\sum_{i=1}^n(y_i-\\bar{y}_i)^2}$",
"$1-e^{-\\frac{2}{n}({loglike(M_e) - loglike(M_n))}}$",
"$\\frac{1-e^{-\\frac{2}{n}({loglike(M_e) - loglike(M_n))}}}{1-e^{\\frac{2*loglike(M_n)}{n}}}$"),
Commentaire = c("Le rapport des *loglikelihood*, très proche de la déviance expliquée.", "Version ajustée du R^2^ de McFadden tenant compte du nombre de paramètres (*k*) dans le modèle.", "Rapport entre la somme des résidus classiques au carré (numérateur) et de la somme des écarts au carré à la moyenne (dénominateur). Notez que pour un GLM gaussien, ce pseudo-R^2^ est identique au R^2^ classique.", "Transformation de la déviance afin de la mettre sur une échelle de 0 à 1 (mais ne pouvant atteindre exactement 1).", "Ajustement du R^2^ de Cox et Snell pour que l’échelle de valeurs possibles puisse comporter 1 (attention, car les valeurs de ce R^2^ tendent à être toujours plus fortes que les autres).")
)
knitr::kable(df,
format.args = list(decimal.mark = ',', big.mark = " "),
col.names = c("Nom" , "Formule", "Commentaire"),
col.to.resize = 3,
col.width = "8cm")
```
En dehors du pseudo-R^2^ de McFadden ajusté, aucune de ces mesures ne tient compte de la complexité du modèle. Il est cependant important de les reporter, car des valeurs très faibles indiquent vraisemblablement un modèle avec une moindre capacité informative. À l’inverse, des valeurs trop fortes pourraient indiquer un problème de surajustement (**voir encadré sur le principe de parcimonie**, [section @sec-0732]).
##### Critère d'information d'Akaike (AIC) {#sec-081513}
Probablement l’indicateur le plus répandu, sa formule est relativement simple, car il s’agit seulement d’un ajustement de la déviance :
$$
AIC = D(M_e) + 2K
$$ {#eq-glm10}
avec *K* le nombre de paramètres à estimer dans le modèle (coefficients, paramètres de distribution, etc.).
L'AIC n’a pas d’interprétation directe, mais permet de comparer deux modèles imbriqués ([section @sec-0732]). Plus l'AIC est petit, mieux le modèle est ajusté. L’idée derrière cet indicateur est relativement simple. Si la déviance *D* est grande, alors le modèle est mal ajusté. Ajouter des paramètres (des coefficients pour de nouvelles variables *X*, par exemple) ne peut que réduire *D*, mais cette réduction n’est pas forcément suffisamment grande pour justifier la complexification du modèle. L'AIC pondère donc *D* en lui ajoutant 2 fois le nombre de paramètres du modèle. Un modèle plus simple (avec moins de paramètres) parvenant à une même déviance est préférable à un modèle complexe (principe de parcimonie ou du rasoir d’Ockham), ce que permet de « quantifier » l'AIC. Attention, l'AIC **ne peut pas être utilisé pour comparer des modèles non imbriqués**. Notez que d’autres indicateurs similaires comme le WAIC, le BIC et le DIC sont utilisés dans un contexte d’inférence bayésienne. Retenez simplement que ces indicateurs sont conceptuellement proches du AIC et s’interprètent (à peu de choses près) de la même façon.
#### Comparaison des distributions originales et prédites {#sec-08152}
Une façon rapide de vérifier si un modèle est mal ajusté est de comparer la forme de la distribution originale et celle capturée par le modèle. L’idée est la suivante : si le modèle est bien ajusté aux données, il est possible de se servir de celui-ci pour générer de nouvelles données dont la distribution ressemble à celle des données originales. Si une différence importante est observable, alors les résultats du modèle ne sont pas fiables, car le modèle échoue à reproduire le phénomène étudié. Cette lecture graphique ne permet pas de s’assurer que le modèle est valide ou bien ajusté, mais simplement d’écarter rapidement les mauvais candidats. Notez que cette méthode ne s’applique pas lorsque la variable modélisée est binaire, multinomiale ou ordinale.
Le graphique à réaliser comprend donc la distribution de la variable dépendante *Y* (représentée avec un histogramme ou un graphique de densité) et plusieurs distributions simulées à partir du modèle. Cette approche est plus répandue dans la statistique bayésienne, mais elle reste pertinente dans l’approche fréquentiste. Il est rare de reporter ces figures, mais elles doivent faire partie de votre diagnostic.
::: bloc_attention
::: bloc_attention-header
::: bloc_attention-icon
:::
**Distinction entre simulation et prédiction**
:::
::: bloc_attention-body
Notez ici que **simuler des données** à partir d’un modèle et **effectuer des prédictions** à partir d’un modèle sont deux opérations différentes. Prédire une valeur à partir d’un modèle revient simplement à appliquer son équation de régression à des données. Si nous réutilisons les mêmes données, la prédiction renvoie toujours le même résultat, il s’agit de la partie systématique (ou déterministe) du modèle.
Pour illustrer cela, admettons que nous avons ajusté un modèle GLM de type gaussien (fonction de lien identitaire) avec trois variables continues $X_1$, $X_2$ et $X_3$ et des coefficients respectifs de 0,5, 1,2 et 1,8 ainsi qu’une constante de 7. Nous pouvons utiliser ces valeurs pour prédire la valeur attendue de $Y$ quand $X_1= 3$, $X_2= 5$ et $X_3 = 5$ :
$\mbox{Prédiction} = \mbox{7 + 3}\times \mbox{0,5 + 5}\times \mbox{1,2 + 5}\times\mbox{1,8 = 23,5}$
En revanche, simuler des données à partir d’un modèle revient à ajouter la dimension stochastique (aléatoire) du modèle. Puisque notre modèle GLM est gaussien, il comporte un paramètre $\sigma$ (son écart-type); admettons, pour cet exemple, qu’il est de 1,2. Ainsi, avec les données précédentes, il est possible de simuler un ensemble infini de valeurs dont la distribution est la suivante : $Normal(\mu = \mbox{23,5, } \sigma = \mbox{1,2})$. 95 % du temps, ces valeurs simulées se trouveront dans l’intervalle $\mbox{[21,1-25,9]}$ ($\mu - 2\sigma \text{; } \mu + 2\sigma$), puisque cette distribution est normale. Les valeurs simulées dépendent donc de la distribution choisie pour le modèle et de l’ensemble des paramètres du modèle, pas seulement de l’équation de régression.
Si vous aviez à ne retenir qu’une seule phrase de ce bloc, retenez que la prédiction ne se réfère qu’à la partie systématique du modèle (équation de régression), alors que la simulation incorpore la partie stochastique (aléatoire) de la distribution du modèle. Deux prédictions effectuées sur des données identiques donnent nécessairement des résultats identiques, ce qui n’est pas le cas pour la simulation.
:::
:::
#### Comparaison des prédictions du modèle avec les valeurs originales {#sec-08153}
Les prédictions d’un modèle devraient être proches des valeurs réelles observées. Si ce n’est pas le cas, alors le modèle n’est pas fiable et ses paramètres ne sont pas informatifs. Dépendamment de la nature de la variable modélisée (quantitative ou qualitative), plusieurs approches peuvent être utilisées pour quantifier l’écart entre valeurs réelles et valeurs prédites.
##### Pour une variable quantitative {#sec-081531}
La mesure la plus couramment utilisée pour une variable quantitative est l’erreur moyenne quadratique (*Root Mean Square Error* – RMSE en anglais).
$$
RMSE = \sqrt{\frac{\sum_{i=1}^n(y_i - \hat{y_i})^2}{n}}
$$ {#eq-glm11}
Il s’agit de la racine carrée de la moyenne des écarts au carré entre valeurs réelles et prédites. Le RMSE est exprimé dans la même unité que la donnée originale et nous donne une indication sur l’erreur moyenne de la prédiction du modèle. Admettons, par exemple, que nous modélisons les niveaux de bruit environnemental en ville en décibels et que notre modèle de régression ait un RMSE de 3,5. Cela signifierait qu’en moyenne notre modèle se trompe de 3,5 décibels (erreur pouvant être négative ou positive), ce qui serait énorme (3 décibels correspondent à une multiplication par deux de l’intensité sonore) et nous amènerait à reconsidérer la fiabilité du modèle. Notez que l’usage d’une moyenne quadratique plutôt qu’une moyenne arithmétique permet de donner plus d’influence aux larges erreurs et donc de pénaliser davantage des modèles faisant parfois de grosses erreurs de prédiction. Le RMSE est donc très sensible à la présence de valeurs aberrantes.
À la place de la moyenne quadratique, il est possible d’utiliser la simple moyenne arithmétique des valeurs absolues des erreurs (MAE). Cette mesure est cependant moins souvent utilisée :
$$
MAE = \frac{\sum_{i=1}^n|y_i - \hat{y_i|}}{n}
$$ {#eq-glm12}
Ces deux mesures peuvent être utilisées pour comparer la capacité de prédiction de deux modèles appliqués aux mêmes données, même s’ils ne sont pas imbriqués. Elles ne permettent cependant pas de prendre en compte la complexité du modèle. Un modèle plus complexe aura toujours des valeurs de RMSE et de MAE plus faibles.
##### Pour une variable qualitative {#sec-081532}
Lorsque l’on modélise une variable qualitative, une erreur revient à prédire la mauvaise catégorie pour une observation. Il est ainsi possible de compter, pour un modèle, le nombre de bonnes et de mauvaises prédictions et d’organiser cette information dans une **matrice de confusion**. Cette dernière prend la forme suivante pour un modèle binaire :
```{r}
#| label: tbl-confusmat1
#| tbl-cap: Exemple de matrice de confusion
#| echo: false
#| message: false
#| warning: false
df <- data.frame(
C1 = c("A" , "B", "Total (%)"),
C2 = c("15" , "5", "20 (46,6)"),
C3 = c("3" , "20", "23 (53,5)"),
C4 = c("18 (41,9)" , "25 (51,1)", "43 (81,4)")
)
knitr::kable(df,
format.args = list(decimal.mark = ',', big.mark = " "),
col.names = c("Valeur prédite / Valeur réelle" , "A", "B", "Total (%)"),
align= c("l", "r", "r", "r"))
```
En colonne du @tbl-confusmat1, nous avons les catégories observées et en ligne, les catégories prédites. La diagonale représente les prédictions correctes. Dans le cas présent, le modèle a bien catégorisé 35 (15 + 20) observations sur 43, soit une précision totale de 81,4 %; huit sont mal classifiées (18,6 %); cinq avec la modalité A ont été catégorisées comme des B, soit 20 % des A, et seules trois B ont été catégorisées comme des A (13 %).
La matrice ci-dessus (@tbl-confusmat1) ne comporte que deux catégories possibles puisque la variable *Y* modélisée est binaire. Il est facile d'étendre le concept de matrice de confusion au cas des variables avec plus de deux modalités (multinomiale). Le @tbl-confusmat2 est un exemple de matrice de confusion multinomiale.
```{r}
#| label: tbl-confusmat2
#| tbl-cap: Exemple de matrice de confusion multinomiale
#| echo: false
#| message: false
#| warning: false
#| out-width: "50%"
df <- data.frame(
C1 = c("A" , "B" , "C", "D", "Total (%)"),
C2 = c("15" , "5", "2", "1" , "23 (18,1)"),
C3 = c("3", "20", "10", "0", "33 (25,7)"),
C4 = c("1", "2", "25", "5", "33 (25,7)"),
C5 = c("5", "12" , "8" , "14" , "39 (30,5)"),
C6 = c("24 (18,7)", "39 (30,4)", "45 (35,2)", "20 (15,6)", "128")
)
knitr::kable(df,
format.args = list(decimal.mark = ',', big.mark = " "),
col.names = c("Valeur prédite / Valeur réelle" , "A", "B" , "C", "D", "Total (%)"),
align= c("l", "r", "r", "r", "r", "r"))
```
Trois mesures pour chaque catégorie peuvent être utilisées pour déterminer la capacité de prédiction du modèle :
* La précision (*precision* en anglais), soit le nombre de fois où une catégorie a été correctement prédite, divisé par le nombre de fois où la catégorie a été prédite.
* Le rappel (*recall* en anglais), soit le nombre de fois où une catégorie a été correctement prédite divisé par le nombre de fois où elle se trouve dans les données originales.
* Le score *F1* est la moyenne harmonique entre la précision et le rappel, soit :
$$
\text{F1} = 2 \times \frac{\text{précision} \times \text{rappel}}{\text{précision} + \text{rappel}}
$$ {#eq-glm13}
Il est possible de calculer les moyennes pondérées des différents indicateurs (macro-indicateurs) afin de disposer d’une valeur d’ensemble pour le modèle. La pondération est faite en fonction du nombre de cas observé de chaque catégorie; l’idée étant qu’il est moins grave d’avoir des indicateurs plus faibles pour des catégories moins fréquentes. Cependant, il est tout à fait possible que cette pondération ne soit pas souhaitable. C’est par exemple le cas dans de nombreuses études en santé portant sur des maladies rares où l’attention est concentrée sur ces catégories peu fréquentes.
Le **coefficient de Kappa** (variant de 0 à 1) peut aussi être utilisé pour quantifier la fidélité générale de la prédiction du modèle. Il est calculé avec l'@eq-kappEq :
$$
k = \frac{Pr(a)-Pr(e)}{1-Pr(e)}
$$ {#eq-kappEq}
avec $Pr(a)$ la proportion d'accords entre les catégories observées et les catégories prédites, et $Pr(e)$ la probabilité d'un accord aléatoire entre les catégories observées et les catégories prédites (@eq-kappEq2).
$$
Pr(e) = \sum^{J}_{j=1} \frac{Cnt_{prédit}(j)}{n\times2} \times \frac{Cnt_{réel}(j)}{n\times2}
$$ {#eq-kappEq2}
avec *n* le nombre d'observations, $Cnt_{prédit}(j)$ le nombre de fois où le modèle prédit la catégorie *j* et $Cnt_{réel}(j)$ le nombre de fois où la catégorie *j* a été observée.
Pour l'interprétation du coefficient de Kappa, référez-vous au @tbl-Kappvals.
```{r}
#| label: tbl-Kappvals
#| tbl-cap: Inteprétation des valeurs du coefficient de Kappa
#| echo: false
#| message: false
#| warning: false
indicators <- data.frame(
K = c("< 0", "0 - 0,20", "0,21 - 0,40", "0,41 - 0,60", "0,61 - 0,80", "0,81 - 1"),
Inter = c("Désaccord", "Accord très faible", "Accord faible", "Accord modéré", "Accord fort", "Accord presque parfait")
)
knitr::kable(indicators,
format.args = list(decimal.mark = ',', big.mark = " "),
col.names = c("K" , "Interprétation"),
align= c("l", "l"))
```
Enfin, un test statistique basé sur la distribution binomiale peut être utilisé pour vérifier que le modèle atteint un niveau de précision supérieur au seuil de non-information. Ce seuil correspond à la proportion de la modalité la plus présente dans le jeu de données. Dans la matrice de confusion utilisée dans le @tbl-Kappvals, ce seuil est de 30,5 % (catégorie D), ce qui signifie qu’un modèle prédisant tout le temps la catégorie D aurait une précision de 30,5 % pour cette catégorie. Il est donc nécessaire que notre modèle fasse mieux que ce seuil.
Dans le cas de la matrice de confusion du @tbl-confusmat2, nous obtenons donc les valeurs affichées dans le @tbl-confusIndic.
```{r}
#| label: tbl-confusIndic
#| tbl-cap: Indicateurs de qualité de prédiction
#| echo: false
#| message: false
#| warning: false
indicators <- data.frame(
C1 = c("A", "B", "C", "D", "macro", "Kappa", "Valeur de p (précision > NIR)"),
C2 = c("65,2", "60,6", "75,8", "35,9", "57,8", "0,44", "< 0,0001"),
C3 = c("31,3", "25,6", "27,8", "35,0", "30,0", "", ""),
C4 = c("42,3", "36,0", "40,7" ,"35,4" ,"38,2" , "" , "")
)
knitr::kable(indicators,
format.args = list(decimal.mark = ',', big.mark = " "),
col.names = c("", "précision", "rappel", "F1"),
align= c("l", "r", "r", "r"))
```
À la lecture du @tbl-confusIndic, nous remarquons que :
* La catégorie D est la moins bien prédite des quatre catégories (faible précision et faible rappel).
* La catégorie C a une forte précision, mais un faible rappel, ce qui signifie que de nombreuses observations étant originalement des A, B ou D ont été prédites comme des C. Ce constat est également vrai pour la catégorie B.
* Le coefficient de Kappa indique un accord modéré entre les valeurs originales et la prédiction.
* La probabilité que la précision du modèle ne dépasse pas le seuil de non-information est inférieure à 0,001, indiquant que le modèle à une précision supérieure à ce seuil.
### Comparaison de deux modèles GLM {#sec-0816}
Tel qu'abordé dans le chapitre sur les régressions linéaires classiques, il est courant de comparer plusieurs modèles imbriqués ([section @sec-0732]). Cette procédure permet de déterminer si l'ajout d'une ou de plusieurs variables contribue à significativement améliorer le modèle. Il est possible d'appliquer la même démarche aux GLM à l'aide du test de rapport de vraisemblance (*likelihood ratio test*). Le principe de base de ce test est de comparer le *likelihood* de deux modèles GLM imbriqués; la valeur de ce test se calcule avec l'équation suivante :
$$
LR = 2(loglik(M_2) - loglik(M_1))
$$ {#eq-glm14}
avec $M_2$ un modèle reprenant toutes les variables du modèle $M_1$, impliquant donc que $loglik(M_2) >= loglik(M_1)$.
Avec ce test, nous supposons que le modèle $M_2$, qui comporte plus de paramètres que le modèle $M_1$, devrait être mieux ajusté aux données. Si c'est bien le cas, la différence entre les *loglikelihood* de deux modèles devrait être supérieure à zéro. La valeur calculée *LR* suit une distribution du khi-deux avec un nombre de degrés de liberté égal au nombre de paramètres supplémentaires dans le modèle $M_2$ comparativement à $M_1$. Avec ces deux informations, il est possible de déterminer la valeur de *p* associée à ce test et de déterminer si $M_2$ est significativement mieux ajusté que $M_1$ aux données. Notez qu'il existe aussi deux autres tests (test de Wald et test de Lagrange) ayant la même fonction. Il s'agit, dans les deux cas, d'approximation du test de rapport des vraisemblances dont la puissance statistique est inférieure au test de rapport de vraisemblance [@NeymanLemma].
Dans les prochaines sections, nous décrivons les modèles GLM les plus couramment utilisés. Il en existe de nombreuses variantes que nous ne pouvons pas toutes décrire ici. L'objectif est de comprendre les rouages de ces modèles afin de pouvoir, en cas de besoin, transposer ces connaissances sur des modèles plus spécifiques. Pour faciliter la lecture de ces sections, nous vous proposons une carte d’identité de chacun des modèles présentés. Elles contiennent l’ensemble des informations pertinentes à retenir pour chaque modèle.
## Modèles GLM pour des variables qualitatives {#sec-082}
Nous abordons en premier les principaux GLM utilisés pour modéliser des variables binaires, multinomiales et ordinales. Prenez bien le temps de saisir le fonctionnement du modèle logistique binomial, car il sert de base pour les trois autres modèles présentés.
### Modèle logistique binomial {#sec-0821}
Le modèle logistique binomial est une généralisation du modèle de Bernoulli que nous avons présenté dans l’introduction de cette section. Le modèle logistique binomiale couvre donc deux cas de figure :
1. La variable observée est binaire (0 ou 1). Dans ce cas, le modèle logistique binomiale devient un simple modèle de Bernoulli.
2. La variable observée est un comptage (nombre de réussites) et nous disposons d’une autre variable avec le nombre de réplications de l’expérience. Par exemple, pour chaque intersection d’un réseau routier, nous pourrions avoir le nombre de décès à vélo (variable *Y* de comptage) et le nombre de collisions vélo / automobile (variable quantifiant le nombre d’expériences, chaque collision étant une expérience).
Spécifiquement, nous tentons de prédire le paramètre *p* de la distribution binomiale à l’aide de notre équation de régression et de la fonction logistique comme fonction de lien. Notez ici que cette fonction de lien influence directement l'interprétation des paramètres du modèle. Pour rappel, cette fonction est définie comme :
$$
g(x) = ln(\frac{x}{1-x})
$$
avec $ln$ étant le logarithme naturel.
Au-delà de sa propriété mathématique assurant que $g(x) \in \mathopen[0,1\mathclose]$, cette fonction offre une interprétation intéressante. La partie $\frac{x}{1-x}$ est une cote et s'interprète en termes de chances d'observer un évènement. Par exemple, dans le cas des accidents de cyclistes, si la probabilité d'observer un décès suite à une collision est de 0,1, alors la cote de cet évènement est $\frac{\frac{1}{10}}{\frac{9}{10}} = \frac{1}{9}$ soit un contre neuf. Dans un modèle GLM logistique, les coefficients ajustés pour les variables indépendantes représentent des **logarithmes de rapport de cote**, car ils comparent les chances d'observer l'évènement (y = 1) en fonction des valeurs des variables indépendantes.
```{r}
#| label: tbl-binomidentity
#| tbl-cap: Carte d'identité du modèle logistique binomial
#| echo: false
#| message: false
#| warning: false
model_formula <- "$Y \\sim Binomial(p)$ \\newline $g(p) = \\beta_0 + \\beta X$ \\newline $g(x) = log(\\frac{x}{1-x})$"
tableau <- data.frame(
C1 = c("Type de variable dépendante" , "Distribution utilisée", "Formulation", "Fonction de lien", "Paramètre modélisé", "Paramètres à estimer", "Conditions d'application"),
C2 = c("Variable binaire (0 ou 1) ou comptage de réussite à une expérience (ex : 3 réussites sur 5 expériences)", "Binomiale", model_formula , "Logistique", "p", "$\\beta_0$, $\\beta$", "Non-séparation complète, absence de sur-dispersion ou de sous-dispersion")
)
knitr::kable(tableau,
format.args = list(decimal.mark = ',', big.mark = " "),
col.names = NULL,
col.to.resize = c(2),
col.width = "8cm")
```
#### Interprétation des paramètres {#sec-08211}
Les seuls paramètres à estimer du modèle sont les coefficients $\beta$ et la constante $\beta_0$. La fonction de lien logistique transforme la valeur de ces coefficients, en conséquence, **ils ne peuvent plus être interprétés directement**. $\beta_0$ et $\beta$ sont exprimés dans une unité particulière: des logarithmes de rapports de cote (*log odd ratio*). Le rapport de cote est relativement facile à interpréter contrairement à son logarithme. Pour l’obtenir, il suffit d’utiliser la fonction exponentielle (l’inverse de la fonction logarithme) pour passer des log rapports de cote à de simples rapports de cote.
Donc si $exp(\beta)$ est inférieur à 1, il réduit les chances d’observer l’évènement et inversement si $exp(\beta)$ est supérieur à 1.
Par exemple, admettons que nous ayons un coefficient $\beta_1$ de 1,2 pour une variable $X_1$ dans une régression logistique. Il est nécessaire d'utiliser son exponentiel pour l'interpréter de façon intuitive. $exp\mbox{(1,2)} = \mbox{3,32}$, ce qui signifie que lorsque $X_1$ augmente d'une unité, les chances d'observer 1 plutôt que 0 comme valeur de *Y* sont multipliées par 3,32. Admettons maintenant que $\beta_1$ vaille −1,2, nous calculons donc $exp\mbox{(-1,2) = 0,30}$, ce qui signifie qu'à chaque augmentation d'une unité de $X_1$, les chances d'observer 1 plutôt que 0 comme valeur de *Y* sont multipliées par 0,30. En d'autres termes,les chances d'observer 1 plutôt que 0 sont divisées par 3,33 ($\mbox{1}/\mbox{0,30} = \mbox{3,33}$), soit une diminution de 70 % ($\mbox{1}-\mbox{0,3} = \mbox{0,7}$) des chances d'observer 1 plutôt que 0.
::: bloc_aller_loin
::: bloc_aller_loin-header
::: bloc_aller_loin-icon
:::
**Les rapports de cotes**
:::
::: bloc_aller_loin-body
Le rapport de cote ou rapport des chances est une mesure utilisée pour exprimer l’effet d’un facteur sur une probabilité. Il est très utilisé dans le domaine de la santé, mais aussi des paris. Prenons un exemple concret avec le port du casque à vélo. Si sur 100 accidents impliquant des cyclistes portant un casque, nous observons seulement 3 cas de blessures graves à la tête, contre 15 dans un second groupe de 100 cyclistes ne portant pas de casque, nous pouvons calculer le rapport de cote suivant :
$$
\frac{p(1-q)}{q(1-p)} = \frac{\mbox{0,15} \times (\mbox{1}-\mbox{0,03})}{\mbox{0,03} \times (\mbox{1}-\mbox{0,15})} = \mbox{5,71}
$$
avec *p* la probabilité d'observer le phénomène (ici la blessure grave à la tête) dans le groupe 1 (ici les cyclistes sans casque) et *q* la probabilité d'observer le phénomène dans le groupe 2 (ici les cyclistes avec un casque). Ce rapport de cote indique que les cyclistes sans casques ont 5,71 fois plus de risques de se blesser gravement à la tête lors d’un accident comparativement aux cyclistes portant un casque.
:::
:::
#### Conditions d'application {#sec-08212}
La non-séparation complète signifie qu’aucune des variables *X* n'est, à elle seule, capable de parfaitement distinguer les deux catégories 0 et 1 de la variable *Y*. Dans un tel cas de figure, les algorithmes d’ajustement utilisés pour estimer les paramètres des modèles sont incapables de converger. Notez aussi l’absurdité de créer un modèle pour prédire une variable *Y* si une variable *X* est capable à elle seule de la prédire à coup sûr. Ce problème est appelé un effet de Hauck-Donner. Il est assez facile de le repérer, car la plupart du temps les fonctions de R signalent ce problème (message d'erreur sur la convergence). Sinon, des valeurs extrêmement élevées ou faibles pour certains rapports de cote peuvent aussi indiquer un effet de Hauck-Donner.
La sur-dispersion est un problème spécifique aux distributions n’ayant pas de paramètre de dispersion (binomiale, de Poisson, exponentielle, etc.), pour lesquelles la variance dépend directement de l'espérance. La sur-dispersion désigne une situation dans laquelle les résidus (ou erreurs) d'un modèle sont plus dispersés que ce que suppose la distribution utilisée. À l’inverse, il est aussi possible (mais rare) d’observer des cas de sous-dispersion (lorsque la dispersion des résidus est plus petite que ce que suppose la distribution choisie). Ce cas de figure se produit généralement lorsque le modèle parvient à réaliser une prédiction trop précise pour être crédible. Si vous rencontrez une forte sous-dispersion, cela signifie souvent que l’une de vos variables indépendantes provoque une séparation complète. La meilleure option, dans ce cas, est de supprimer la variable en question du modèle. La variance attendue d’une distribution binomiale est $nb \times p \times(1-p)$, soit le produit entre le nombre de tirages, la probabilité de réussite et la probabilité d'échec. À titre d'exemple, si nous considérons une distribution binomiale avec un seul tirage et 50 % de chances de réussite, sa variance serait : $1 \times \mbox{0,5} \times \mbox{(1}-\mbox{0,5}) = \mbox{0,25}$.
Plusieurs raisons peuvent expliquer la présence de sur-dispersion dans un modèle :
* il manque des variables importantes dans le modèle, conduisant à un mauvais ajustement et donc une sur-dispersion des erreurs;
* les observations ne sont pas indépendantes, impliquant qu’une partie de la variance n’est pas contrôlée et augmente les erreurs;
* la probabilité de succès de chaque expérience varie d’une répétition à l’autre (différentes distributions).
La conséquence directe de la sur-dispersion est la sous-estimation de la variance des coefficients de régression. En d’autres termes, la sur-dispersion conduit à sous-estimer notre incertitude quant aux coefficients obtenus et réduit les valeurs de *p* calculées pour ces coefficients. Les risques de trouver des résultats significatifs à cause des fluctuations d'échantillonnage augmentent.
Pour détecter une sur-dispersion ou une sous-dispersion dans un modèle logistique binomial, il est possible d'observer les résidus de déviance du modèle. Ces derniers sont supposés suivre une distribution du khi-deux avec *n−k* degrés de liberté (avec *n* le nombre d'observations et *k* le nombre de coefficients dans le modèle). Par conséquent, la somme des résidus de déviance d'un modèle logistique binomiale divisée par le nombre de degrés de liberté devrait être proche de 1. Une légère déviation (jusqu'à 0,15 au-dessus ou au-dessous de 1) n'est pas alarmante; au-delà, il est nécessaire d'ajuster le modèle.
Notez que si la variable *Y* modélisée est exactement binaire (chaque expérience est indépendante et n'est composée que d'un seul tirage) et que le modèle utilise donc une distribution de Bernoulli, le test précédent pour détecter une éventuelle sur-dispersion n'est pas valide. @hilbe2009logistic parle de sur-dispersion implicite pour le modèle de Bernoulli et recommande notamment de toujours ajuster les erreurs standards des modèles utilisant des distributions de Bernoulli, binomiale et de Poisson. L'idée ici est d'éviter d'être trop optimiste face à l'incertitude du modèle sur les coefficients et de l'ajuster en conséquence. Pour cela, il est possible d'utiliser des quasi-distributions ou des estimateurs robustes [@zeileis2004econometric]. Notez que si le modèle ne souffre pas de sur ou sous-dispersion, ces ajustements produisent des résultats équivalents aux résultats non ajustés.
#### Exemple appliqué dans R {#sec-08213}
**Présentation des données**
Pour illustrer le modèle logistique binomial, nous utilisons ici un jeu de données proposé par l’Union européenne : [l’enquête de déplacement sur la demande pour des systèmes de transports innovants](https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/P82V9X). Pour cette enquête, un échantillon de 1 000 individus représentatifs de la population a été constitué dans chacun des 26 États membres de l’UE, soit un total de 26 000 observations. Pour chaque individu, plusieurs informations ont été collectées relatives à la catégorie socioprofessionnelle, le mode de transport le plus fréquent, le temps du trajet de son déplacement le plus fréquent et son niveau de sensibilité à la cause environnementale. Nous modélisons ici la probabilité qu’un individu déclare utiliser le plus fréquemment le vélo comme moyen de transport. Les variables explicatives sont résumées au @tbl-binomdata. Il existe bien évidemment un grand nombre de facteurs individuels qui influence la prise de décision sur le mode de transport. Les résultats de ce modèle ne doivent donc pas être pris avec un grand sérieux; il est uniquement construit à des fins pédagogiques, sans cadre conceptuel solide.
```{r}
#| label: tbl-binomdata
#| tbl-cap: Variables indépendantes utilisées pour prédire le mode de transport le plus utilisé
#| echo: false
#| message: false
#| warning: false
tableau <- data.frame(
C1 = c("Pays", "Sexe", "Age", "Education", "StatutEmploi", "Revenu", "Residence", "Duree", "ConsEnv"),
C2 = c("Pays de résidence", "Sexe biologique", "Âge biologique", "Niveau d’éducation maximum atteint", "Employé ou non", "Niveau de revenu autodéclaré", "Lieu de résidence", "Durée du voyage le plus fréquent autodéclarée (en minutes)", "Préoccupation environnementale"),
C3 = c("Variable multinomiale", "Variable binaire", "Variable continue", "Variable multinomiale", "Variable binaire", "Variable multinomiale", "Variable multinomiale", "Variable continue", "Variable ordinale"),
C4 = c("Le nom d’un des 26 pays membres de l’UE", "Homme ou femme", "L’âge en nombre d’années variant de 16 à 84 ans dans le jeu de données", "Premier cycle, secondaire inférieur (classes supérieures de l’école élémentaire), secondaire, troisième cycle", "Employé ou non", "Très faible revenu, faible revenu, revenu moyen, revenu élevé, revenu très élevé, sans reponse", "Zone rurale, petite ou moyenne ville (moins de 250 000 habitants), grande ville (entre 250 000 et 1 million d’habitants) , aire métropolitaine (plus d’un million d’habitants)", "Nombre de minutes", "Échelle de Likert de 1 à 10")
)
knitr::kable(tableau,
format.args = list(decimal.mark = ',', big.mark = " "),
col.names = c("Nom de la variable" , "Signification" , "Type de variable" , "Mesure"),
col.to.resize = c(2,4),
col.width = "5cm",
align= c("l", "l", "l", "l"))
```
**Vérification des conditions d'application**
La première étape de la vérification des conditions d'application est de calculer les valeurs du facteur d'inflation de variance (VIF) pour s'assurer de l'absence de multicolinéarité trop forte entre les variables indépendantes. L'ensemble des valeurs de VIF sont inférieures à 5, indiquant l'absence de multicolinéarité excessive dans le modèle.
```{r}
#| message: false
#| warning: false
#| out-width: "50%"
library(car)
# Chargement des données
dfenquete <- read.csv("data/glm/enquete_transport_UE.csv", encoding = 'UTF-8')
dfenquete$Pays <- relevel(as.factor(dfenquete$Pays), ref = "Allemagne")
# Vérification du VIF
model1 <- glm(y ~
Pays + Sexe + Age + Education + StatutEmploi + Revenu +
Residence + Duree + ConsEnv,
family = binomial(link="logit"),
data = dfenquete
)
vif(model1)
```
La seconde étape de vérification est le calcul des distances de Cook et l'identification d'éventuelles valeurs aberrantes (@fig-cookdistGLMbinom).
```{r}
#| label: fig-cookdistGLMbinom
#| fig-align: center
#| fig-cap: Distances de Cook pour le modèle binomial avec toutes les observations
#| warning: false
#| message: false
#| out-width: "60%"
# Calcul et représentation des distances de Cook
cookd <- data.frame(
dist = cooks.distance(model1),
oid = 1:nrow(dfenquete)
)
ggplot(cookd) +
geom_point(aes(x = oid, y = dist ), color = rgb(0.1,0.1,0.1,0.4), size = 1)+
geom_hline(yintercept = 0.002, color = "red")+
labs(x = "observations", y = "distance de Cook") +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank())
```
Le calcul de la distance de Cook révèle un ensemble d'observations se démarquant nettement des autres (délimitées dans la @fig-cookdistGLMbinom) par la ligne rouge). Nous les isolons dans un premier temps pour les analyser.
```{r}
#| message: false
#| warning: false
#| out-width: "80%"
# Isoler les observations avec de très fortes valeurs de Cook
# valeur seuil choisie : 0,002
cas_etranges <- subset(dfenquete, cookd$dist>=0.002)
cat(nrow(cas_etranges),'observations se démarquant dans le modèle')
print(cas_etranges)
```
À la lecture des valeurs pour ces 19 cas étranges, nous remarquons que la plupart des observations proviennent de Malte et de Chypre. Ces deux petites îles constituent des cas particuliers en Europe et devraient vraisemblablement faire l'objet d'une analyse séparée. Nous décidons donc de les retirer du jeu de données. Deux autres observations étranges sont observables en Slovaquie et au Luxembourg. Dans les deux cas, les répondants ont renseigné des temps de trajet fantaisistes de respectivement 775 et 720 minutes. Nous les retirons donc également de l'analyse.
```{r}
#| label: fig-cookdistGLMbinom2
#| fig-align: center
#| fig-cap: Distances de Cook pour le modèle binomial sans les valeurs aberrantes
#| warning: false
#| message: false
#| out-width: "60%"
# Retirer les observations aberrantes
dfenquete2 <- subset(dfenquete, (dfenquete$Pays %in% c("Malte", "Chypre")) == F &
dfenquete$Duree < 400)
# Réajuster le modèle
model2 <- glm(y ~
Pays + Sexe + Age + Education + StatutEmploi + Revenu +
Residence + Duree + ConsEnv,
family = binomial(link="logit"),
data = dfenquete2)
# Recalculer la distance de Cook
cookd <- data.frame(
dist = cooks.distance(model2),
oid = 1:nrow(dfenquete2)
)
ggplot(cookd) +
geom_point(aes(x = oid, y = dist ), color = rgb(0.1,0.1,0.1,0.4), size = 1)+
labs(x = "observations", y = "distance de Cook") +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank())
```
Après avoir retiré ces valeurs aberrantes, nous n'observons plus de nouveaux cas singuliers avec les distances de Cook (@fig-cookdistGLMbinom2).
La prochaine étape de vérification des conditions d'application est l'analyse des résidus simulés. Nous commençons donc par calculer ces résidus et afficher leur histogramme (@fig-residsimbinomHistoUnif).
```{r}
#| label: fig-residsimbinomHistoUnif
#| fig-align: center
#| fig-cap: Distribution des résidus simulés pour le modèle binomial
#| out-width: "60%"
#| message: false
#| warning: false
library(DHARMa)
# Extraire les probabilités prédites par le modèle
probs <- predict(model2, type = "response")
# Calculer 1000 simulations a partir du modele ajuste
sims <- lapply(1:length(probs), function(i){
p <- probs[[i]]
vals <- rbinom(n = 1000, size = 1, prob = p)
})
matsim <- do.call(rbind, sims)
# Utiliser le package DHARMa pour calculer les résidus simulés
sim_res <- createDHARMa(simulatedResponse = matsim,
observedResponse = dfenquete2$y,
fittedPredictedResponse = probs,
integerResponse = TRUE)
ggplot()+
geom_histogram(aes(x = residuals(sim_res)),
bins = 30, fill = "white", color = rgb(0.3,0.3,0.3))+
labs(x = "résidus simulés", y = "fréquence")
```
L'histogramme indique clairement que les résidus simulés suivent une distribution uniforme (@fig-residsimbinomHistoUnif). Il est possible d'aller plus loin dans le diagnostic en utilisant la fonction `plot` sur l'objet `sim_res`. La partie de droite de la figure ainsi obtenue (@fig-residsimbinom) est un diagramme de quantiles-quantiles (ou Q-Q plot). Les points du graphique sont supposés suivre une ligne droite matérialisée par la ligne rouge. Une déviation de cette ligne indique un éloignement des résidus de leur distribution attendue. Trois tests sont également réalisés par la fonction :
* Le premier (Test de Kolmogorov-Smirnov, *KS test*) permet de tester si les points dévient significativement de la ligne droite. Dans notre cas, la valeur de *p* n'est pas significative, indiquant que les résidus ne dévient pas de la distribution uniforme.
* Le second test permet de vérifier la présence de sur ou sous-dispersion. Dans notre cas, ce test n'est pas significatif, n'indiquant aucun problème de sur-dispersion ou de sous-dispersion.
* Le dernier test permet de vérifier si des valeurs aberrantes sont présentes dans les résidus. Une valeur non significative indique une absence de valeurs aberrantes.
Le second graphique permet de comparer les résidus et les valeurs prédites. L'idéal est donc d'observer une ligne droite horizontale au milieu du graphique qui indiquerait une absence de relation entre les valeurs prédites et les résidus (ce que nous observons bien ici).
```{r}
#| eval: false
plot(sim_res)
```
```{r}
#| label: fig-residsimbinom
#| fig-align: center
#| fig-cap: Diagnostic des résidus simulés par le package DHARMa
#| out-width: "80%"
#| echo: false
#| message: false
#| warning: false
invisible(capture.output({
out <- "images/Chap08/fig-residsimbinom.png"
png(filename = out,
res = 300, width = 20, height = 10, units = "cm")
plot(sim_res)
dev.off()
}))
knitr::include_graphics(out, dpi = 300)
```
L'analyse approfondie des résidus nous permet donc de conclure que le modèle respecte les conditions d'application et que nous pouvons passer à la vérification de la qualité d'ajustement du modèle.
**Vérification de la qualité d'ajustement**
Pour calculer les différents R^2^ d'un modèle GLM, nous proposons la fonction suivante :
```{r}
#| message: false
#| warning: false
#| out-width: "100%"
rsqs <- function(loglike.full, loglike.null, full.deviance, null.deviance, nb.params, n){
# Calcul de la déviance expliquée
explained_dev <- 1-(full.deviance / null.deviance)
K <- nb.params
# R2 de McFadden ajusté
r2_faddenadj <- 1- (loglike.full - K) / loglike.null
Lm <- loglike.full
Ln <- loglike.null
# R2 de Cox and Snell
Rcs <- 1 - exp((-2/n) * (Lm-Ln))
# R2 de Nagelkerke
Rn <- Rcs / (1-exp(2*Ln/n))
return(
list("deviance expliquee" = explained_dev,
"McFadden ajuste" = r2_faddenadj,
"Cox and Snell" = Rcs,
"Nagelkerke" = Rn
)
)
}
```
Nous l'utilisons pour l'ensemble des modèles GLM de ce chapitre. Dans le cas du modèle binomial, nous obtenons :
```{r}
#| message: false
#| warning: false
#| out-width: "100%"
# Ajuster un modele null avec seulement une constante
model2.null <- glm(y ~1,
family = binomial(link="logit"),
data = dfenquete2)
# Calculer les R2
rsqs(loglike.full = as.numeric(logLik(model2)), # loglikelihood du modèle complet
loglike.null = as.numeric(logLik(model2.null)), # loglikelihood du modèle nul
full.deviance = deviance(model2), # déviance du modèle complet
null.deviance = deviance(model2.null), # déviance du modèle nul
nb.params = model2$rank, # nombre de paramètres dans le modèle
n = nrow(dfenquete2) # nombre d'observations
)
```
La déviance expliquée par le modèle est de 8,8 %, les pseudos R^2^ de McFadden (ajusté), d'Efron et de Nagelkerke sont respectivement 0,084, 0,069 et 0,124. Toutes ces valeurs sont relativement faibles et indiquent qu'une large partie de la variabilité de *Y* reste inexpliquée.
Pour vérifier la qualité de prédiction du modèle, nous devons comparer les catégories prédites et les catégories réelles de notre variable dépendante et construire une matrice de confusion. Cependant, un modèle GLM binomial prédit **la probabilité d’appartenance au groupe 1** (ici les personnes utilisant le vélo pour effectuer leur déplacement le plus fréquent). Pour convertir ces probabilités prédites en catégories prédites, il faut choisir un seuil de probabilité au-delà duquel nous considérons que la valeur attendue est 1 (cycliste) plutôt que 0 (autre). Un exemple naïf serait de prendre le seuil 0,5, ce qui signifierait que si le modèle prédit qu’une observation a au moins 50 % de chance d’être une personne à vélo, alors nous l’attribuons à cette catégorie. Cependant, cette méthode est rarement optimale; il est donc plus judicieux de fixer le seuil de probabilité en trouvant le point d’équilibre entre la sensibilité (proportion de 1 correctement identifiés) et la spécificité (proportion de 0 correctement identifiés). Ce point d’équilibre est identifiable graphiquement en calculant la spécificité et la sensibilité de la prédiction selon toutes les valeurs possibles du seuil.
```{r}
#| label: fig-equlibresensispeci
#| fig-cap: Point d'équilibre entre sensibilité et spécificité
#| fig-align: center
#| out-width: "80%"
#| #| message: false
#| warning: false
library(ROCR)
# Obtention des prédictions du modèle
prob <- predict(model2, type = "response")
# Calcul de la sensibilité et de la spécificité (package ROCR)
predictions <- prediction(prob, dfenquete2$y)
sens <- data.frame(x = unlist(ROCR::performance(predictions, "sens")@x.values),
y = unlist(ROCR::performance(predictions, "sens")@y.values))
spec <- data.frame(x = unlist(ROCR::performance(predictions, "spec")@x.values),
y = unlist(ROCR::performance(predictions, "spec")@y.values))
# Trouver numériquement la valeur seuil (minimiser la différence absolue
# entre sensibilité et spécificité)
real <- dfenquete2$y
find_cutoff <- function(seuil){
pred <- ifelse(prob>seuil,1,0)
sensi <- sum(real==1 & pred==1) / sum(real==1)
spec <- sum(real==0 & pred==0) / sum(real==0)
return(abs(sensi-spec))
}
prob_seuil <- optimize(find_cutoff, interval = c(0,1), maximum = FALSE)$minimum
cat("Le seuil de probabilité à retenir équilibrant",
"la sensibilité et la spécificité est de", prob_seuil)
# Affichage du graphique
ggplot() +
geom_line(data = sens, mapping = aes(x = x, y = y)) +
geom_line(data = spec, mapping = aes(x = x, y = y, col="red")) +
scale_y_continuous(sec.axis = sec_axis(~., name = "Spécificité")) +
labs(x = "Seuil de probabilité", y = "Sensibilité") +
geom_vline(xintercept = prob_seuil, color = "black", linetype = "dashed") +
annotate(geom = "text", x = prob_seuil, y = 0.01, label = round(prob_seuil,3))+
theme(axis.title.y.right = element_text(colour = "red"), legend.position = "none")
```
Nous constatons à la @fig-equlibresensispeci que si la valeur du seuil est 0 %, alors la prédiction a une sensibilité parfaite (le modèle prédit toujours 1, donc tous les 1 sont détectés); à l’inverse, si le seuil choisi est 100 %, alors la prédiction à une spécificité parfaite (le modèle prédit toujours 0, donc tous les 0 sont détectés). Dans notre cas, la valeur d’équilibre est d'environ 0,148, donc si le modèle prédit une probabilité au moins égale à 14,8 % qu’un individu utilise le vélo pour son déplacement le plus fréquent, nous devons l’attribuer à la catégorie *cycliste*. Avec ce seuil, nous pouvons convertir les probabilités prédites en classes prédites et construire notre matrice de confusion.
```{r}
#| message: false
#| warning: false
#| out-width: "60%"
library(caret) # pour la matrice de confusion
# Calcul des catégories prédites
ypred <- ifelse(predict(model2, type = "response")>0.148,1,0)
info <- confusionMatrix(as.factor(dfenquete2$y), as.factor(ypred))
# Affichage des valeurs brutes de la matrice de confusion
print(info)
```
Les résultats proposés par le *package* `caret` sont exhaustifs; nous vous proposons ici une façon de les présenter dans deux tableaux : l'un présente la matrice de confusion (@tbl-confusmatbinom) et l'autre, les indicateurs de qualité de prédiction (@tbl-confusmatbinom2).
```{r}
#| label: tbl-confusmatbinom
#| tbl-cap: Matrice de confusion pour le modèle binomial
#| echo: false
#| message: false
#| warning: false
mat <- info[[2]]
rs <- rowSums(mat)
rp <- round(rowSums(mat) / sum(mat) * 100,1)
cs <- colSums(mat)
cp <- round(colSums(mat) / sum(mat) * 100,1)
mat2 <- cbind(mat, rs, rp)
mat3 <- rbind(mat2, c(cs, sum(mat), NA), c(cp, NA, NA))
mat4 <- cbind(c("0 (prédit)", "1 (prédit)", "Total", "%"), mat3)
options(knitr.kable.NA = "")
knitr::kable(mat4,
format.args = list(decimal.mark = ',', big.mark = " "),
col.names = c("" , "0 (réel)", "1 (réel)", "Total", "%"),
align= c("l", "r", "r", "r", "r"))
```
D’après ces indicateurs, nous constatons que le modèle a une capacité de prédiction relativement faible, mais tout de même significativement supérieure au seuil de non-information. La valeur de rappel pour la catégorie 1 (cycliste) est faible, indiquant que le modèle a manqué un nombre important de cyclistes lors de sa prédiction.
```{r}
#| label: tbl-confusmatbinom2
#| tbl-cap: Matrice de confusion pour le modèle binomial
#| echo: false
#| message: false
#| warning: false
precision <- diag(mat) / rowSums(mat)
rappel <- diag(mat) / colSums(mat)
F1 <- 2*((precision*rappel)/(precision + rappel))
pval <- round(info[[3]][[6]],4)
if(pval == 0){
pval <- "< 0,0001"
}
macro_scores <- c(weighted.mean(precision, colSums(mat)),
weighted.mean(rappel, colSums(mat)),
weighted.mean(F1, colSums(mat)))
final_table <- rbind(cbind(precision, rappel, F1), c(NA, NA, NA, NA), macro_scores)
final_table <- rbind(final_table, c(info[[3]][[2]], NA, NA), c(info[[3]][[6]], NA, NA))
rnames <- c(rownames(mat),"" , "macro" , "Kappa" , "Valeur de p (précision > NIR)")
final_table <- cbind(rnames, round(final_table,2))
final_table[6,2] <- pval
options(knitr.kable.NA = "")
knitr::kable(final_table,
format.args = list(decimal.mark = ',', big.mark = " "),
col.names = c("" , "Précision", "Rappel", "F1"),
align= c("l", "r", "r", "r"))
```
**Interprétation des résultats du modèle**
L'interprétation des résultats d'un modèle binomial passe par la lecture des rapports de cotes (exponentiel des coefficients) et de leurs intervalles de confiance. Nous commençons donc par calculer la version robuste des erreurs standards des coefficients.
```{r}
#| message: false
#| warning: false
library(dplyr)
library(sandwich) # pour calculer les erreurs standards robustes
covModel2 <- vcovHC(model2, type = "HC0") # méthode HC0, basée sur les résidus
stdErrRobuste <- sqrt(diag(covModel2)) # extraire la diagonale
# Extraction des coefficients
coeffs <- model2$coefficients
# Recalcul des scores Z
zvalRobuste <- coeffs / stdErrRobuste
# Recalcul des valeurs de P
pvalRobuste <- 2 * pnorm(abs(zvalRobuste), lower.tail = FALSE)
# Calcul des rapports de cote
oddRatio <- exp(coeffs)
# Calcul des intervalles de confiance à 95 % des rapports de cote
lowerBound <- exp(coeffs - 1.96 * stdErrRobuste)
upperBound <- exp(coeffs + 1.96 * stdErrRobuste)
# Étoiles pour les valeurs de p
starsp <- case_when(pvalRobuste <= 0.001 ~ "***",
pvalRobuste > 0.001 & pvalRobuste <= 0.01 ~ "**",
pvalRobuste > 0.01 & pvalRobuste <= 0.05 ~ "*",
pvalRobuste > 0.05 & pvalRobuste <= 0.1 ~ ".",
TRUE ~ ""
)
# Compilation des résultats dans un tableau
tableau_binom <- data.frame(
coefficients = coeffs,
rap.cote = oddRatio,
err.std = stdErrRobuste,
score.z = zvalRobuste,
p.val = pvalRobuste,
rap.cote.2.5 = lowerBound,
rap.cote.97.5 = upperBound,
sign = starsp
)
```
Considérant que la variable *Pays* a 24 modalités, il est plus judicieux de présenter ses 23 rapports de cotes sous forme d’un graphique. Nous avons choisi l’Allemagne comme catégorie de référence puisqu’elle fait partie des pays avec une importante part modale pour le vélo sans pour autant constituer un cas extrême comme le Danemark.
```{r}
#| message: false
#| warning: false
#| label: fig-paysbinom
#| fig-align: center
#| fig-cap: Rapports de cote pour les différents pays de l'UE
#| out-width: "80%"
# Isoler les ligne du tableau récapitualtif pour les pays
paysdf <- subset(tableau_binom, grepl("Pays", row.names(tableau_binom), fixed = TRUE))
#paysdf$Pays <- gsub("Pays" , "", row.names(paysdf), fixed = TRUE)
paysdf$Pays <- substr(row.names(paysdf), 5, nchar(row.names(paysdf)))
ggplot(data = paysdf) +
geom_vline(xintercept = 1, color = "red")+ #afficher la valeur de référence
geom_errorbarh(aes(xmin = rap.cote.2.5, xmax = rap.cote.97.5,
y = reorder(Pays, rap.cote)), height = 0)+
geom_point(aes(x = rap.cote, y = reorder(Pays, rap.cote))) +
geom_text(aes(x = rap.cote.97.5, y = reorder(Pays, rap.cote),
label = paste("RC : ", round(rap.cote,2), sep = "")),
size = 3, nudge_x = 0.25)+
labs(x = "Rapports de cote", y = "Pays (référence : Allemagne)")
```
Dans la @fig-paysbinom, la barre horizontale pour chaque pays représente l’intervalle de confiance de son rapport de cotes (le point); plus cette ligne est longue, plus grande est l’incertitude autour de ce paramètre. Lorsque les lignes de deux pays se chevauchent, cela signifie qu’il n’y a pas de différence significative au seuil 0,05 entre les rapports de cotes des deux pays. La ligne rouge tracée à x = 1, représente le rapport de cotes du pays de référence (ici l’Allemagne). Nous constatons ainsi que comparativement à un individu vivant en Allemagne, ceux vivant au Danemark et aux Pays-Bas ont 2,4 fois plus de chances d’utiliser le vélo pour leur déplacement le plus fréquent. Les Pays de l’Ouest (France, Luxembourg, Royaume-Uni, Irlande) et du Sud (Grèce, Italie, Espagne, Portugal) ont en revanche des rapports de cotes plus faibles. En France, les chances qu’un individu utilise le vélo pour son trajet le plus fréquent sont 3,22 (1/0,31) fois plus faibles que si l’individu vivait en Allemagne.
Pour le reste des coefficients et des rapports de cotes, nous les rapportons dans le @tbl-coeffbinom.
```{r}
#| label: tbl-coeffbinom
#| tbl-cap: Résultats du modèle binomial