data(evaluation, package = "hecmulti")
mod0 <- glm(choix ~ force,
data = evaluation,
family = binomial)Séance 5
Ajustement d’un modèle de régression logistique
On considère l’étude 3 de Liu et al. (2026), dont les données sont disponibles via la base de données evaluation dans le paquet hecmulti.
Les participants ont pris part à un jeu en ligne dans lequel ils devaient trouver des fruits parmi un ensemble de trois boîtes, et obtenaient un score lorsqu’ils identifiaient correctement la boîte contenant le fruit. À leur insu, le jeu était truqué et leurs scores déterminés à l’avance. On leur communiquait, en plus de leur score final, une mesure des performances des autres participants, afin de leur indiquer s’ils avaient obtenu de meilleurs ou de moins bons résultats que les autres (force). La moitié des participants était censée recevoir un bonus en fonction de leurs bonnes performances, évaluées par des juges algorithmiques et humains: ils devaient choisir entre les deux (choix). Le médiateur est une mesure de la préférence relative pour leur choix de juge (mediateur).
Nous considérons un modèle logistique additif avec deux variables explicatives pour la variable réponse choix: mediateur est une variable entière sur 0-10 et force qui est catégorielle binaire.
On peut considérer les probabilités tirées du modèle simple mod0 avec une seule variable explicative (force): il prédit simplement la proportion empirique de choix pour humain (ici la catégorie de référence) pour chaque groupe expérimental (score fort ou faible).
predict(mod0,
newdata = data.frame(force = c("fort","faible")),
type = "response") 1 2
0.4188377 0.5030181
# Comparaison avec la proportion empirique
evaluation |>
dplyr::group_by(force) |>
dplyr::summarize(proportion = mean(choix == "humain"))# A tibble: 2 × 2
force proportion
<fct> <dbl>
1 fort 0.419
2 faible 0.503
On peut également entraîner un modèle plus complexe qui ajuste pour la force de conviction envers le choix de l’agent: une valeur de 10 indique une forte préférence, 0 une absence de préférence marquée.
mod1 <- glm(choix ~ mediateur + force,
data = evaluation,
family = binomial)
summary(mod1)
Call:
glm(formula = choix ~ mediateur + force, family = binomial, data = evaluation)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.39418 0.15323 2.572 0.0101 *
mediateur -0.12207 0.02097 -5.821 5.84e-09 ***
forcefaible 0.16386 0.13308 1.231 0.2182
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1374.6 on 995 degrees of freedom
Residual deviance: 1332.5 on 993 degrees of freedom
AIC: 1338.5
Number of Fisher Scoring iterations: 4
On peut obtenir les estimations des coefficients, des erreurs-types et les valeurs-p des tests de Wald dans le tableau avec le générique S3 summary. Puisque le modèle est multiplicatif, le rapport de cote donne une augmentation ou diminution de exp(beta).
summary(mod1) # tableau résumé du modèle
Call:
glm(formula = choix ~ mediateur + force, family = binomial, data = evaluation)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.39418 0.15323 2.572 0.0101 *
mediateur -0.12207 0.02097 -5.821 5.84e-09 ***
forcefaible 0.16386 0.13308 1.231 0.2182
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1374.6 on 995 degrees of freedom
Residual deviance: 1332.5 on 993 degrees of freedom
AIC: 1338.5
Number of Fisher Scoring iterations: 4
exp(coef(mod1)) # estimations de la cote(Intercept) mediateur forcefaible
1.4831727 0.8850821 1.1780534
exp(confint(mod1)) # intervalles de confiance pour la coteWaiting for profiling to be done...
2.5 % 97.5 %
(Intercept) 1.0994222 2.0055723
mediateur 0.8491512 0.9219561
forcefaible 0.9072943 1.5289846
On peut vérifier directement si les coefficients sont significatifs, mais il est parfois utile pour l’inférence de faire un test de rapport de vraisemblance pour comparer des modèles entre eux (dans le cas de variables explicatives avec plusieurs coefficients, si on a une variable catégorielle à \(K>2\) modalités. Le générique anova fait ce test avec une séquence de modèles allant du plus simple au plus complexe.
Par défaut, si on ne met qu’un seul modèle, la décomposition de la déviance (égal à moins deux fois la log-vraisemblance, \(-2\ell\)) se fait en respectant l’ordre dans lequel les variables sont entrés dans la formule: vous pouvez toujours en cas de doute plutôt ajuster explicitement les modèles à comparer. La plupart des logiciels utilisent une décomposition de type III, qui compare le modèle si on exclut une variable ou une interaction, versus le modèle complet. Ça ne respecte pas la hiérarchie effets principaux et interactions.
mod2 <- glm(choix ~ mediateur * force,
data = evaluation,
family = binomial)
anova(mod1, mod2, test = "LR") # test de rapport de vraisemblanceAnalysis of Deviance Table
Model 1: choix ~ mediateur + force
Model 2: choix ~ mediateur * force
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 993 1332.5
2 992 1330.6 1 1.9028 0.1678
car::Anova(mod1, type = 2)Analysis of Deviance Table (Type II tests)
Response: choix
LR Chisq Df Pr(>Chisq)
mediateur 35.058 1 3.201e-09 ***
force 1.514 1 0.2185
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
On voit au résultat que le médiateur explique à lui seul la majeure partie de la réponse.
Comme pour la régression linéaire, la présence de variables explicatives très fortement corrélées peut impacter l’interprétabilité. On peut vérifier la colinéarité à l’aide des facteurs d’inflation de la variance généralisés. De fortes valeurs (le point de référence est 1 pour l’absence de collinéarité) sont indicatrices de colinéarité, ce qui n’est pas le cas ici.
car::vif(mod1)mediateur force
1.049863 1.049863
Modèle logistique pour proportions
On considère une régression logistique pour données binomiales, qui vise à modéliser le taux de succès aux examens pratique de conduite en Grande-Bretagne en fonction de la région et de la taille du centre (ici coupée en catégories pour simplifier). On commence par créer un facteur qui représente la taille du centre, et on met comme référence la région urbaine de la capitale, Londres.
data(gbconduite, package = "hecmodstat")
gbconduite <- gbconduite |>
dplyr::mutate(taille = factor(
cut(gbconduite$total, c(0, 500, 1000, Inf)),
labels = c("petit", "moyen", "grand")),
region = relevel(factor(region), ref = "London")
)Évidemment, il y a une grande corrélation entre région et nombre d’examens. On peut regarder le décompte de la taille des centres par régions. On note que les régions urbaines n’ont pas de petit centres: cela nous empêchera d’ajuster une interaction. La plupart des petits centres (moins de 500 examens par année) sont situés en Écosse.
with(gbconduite, table(taille))taille
petit moyen grand
126 54 512
with(gbconduite, table(region, taille)) taille
region petit moyen grand
London 6 4 48
East Midlands 3 3 40
East of England 0 0 54
North East England 8 5 29
North West England 2 3 63
Scotland 94 17 41
South East England 0 0 78
South West England 0 6 44
Wales 9 10 29
West Midlands 2 6 54
Yorkshire and the Hu 2 0 32
On considère le taux de réussite. Dans la fonction glm, il faudra mettre comme variable réponse à gauche du ~ une matrice à deux colonnes, avec chacune le nombre de succès (colonne 1) et d’échecs (colonne 2).
mod_gbconduite <- glm(
cbind(reussite, total - reussite) ~ sexe + region + taille,
data = gbconduite,
family = binomial(link = "logit")
)
summary(mod_gbconduite)
Call:
glm(formula = cbind(reussite, total - reussite) ~ sexe + region +
taille, family = binomial(link = "logit"), data = gbconduite)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.017615 0.014757 -1.194 0.233
sexehomme 0.288708 0.003133 92.138 < 2e-16 ***
regionEast Midlands 0.246407 0.007087 34.770 < 2e-16 ***
regionEast of England 0.215940 0.006544 33.000 < 2e-16 ***
regionNorth East England 0.404858 0.008382 48.303 < 2e-16 ***
regionNorth West England 0.207555 0.006208 33.435 < 2e-16 ***
regionScotland 0.231355 0.007335 31.540 < 2e-16 ***
regionSouth East England 0.228954 0.005720 40.024 < 2e-16 ***
regionSouth West England 0.339929 0.007176 47.368 < 2e-16 ***
regionWales 0.365998 0.008580 42.657 < 2e-16 ***
regionWest Midlands 0.045230 0.006463 6.998 2.6e-12 ***
regionYorkshire and the Hu 0.089605 0.007389 12.127 < 2e-16 ***
taillemoyen -0.259958 0.017000 -15.291 < 2e-16 ***
taillegrand -0.488828 0.014094 -34.683 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 37885 on 691 degrees of freedom
Residual deviance: 21241 on 678 degrees of freedom
AIC: 26577
Number of Fisher Scoring iterations: 3
exp(cbind(coef(mod_gbconduite), confint(mod_gbconduite)))Waiting for profiling to be done...
2.5 % 97.5 %
(Intercept) 0.9825397 0.9545468 1.0113909
sexehomme 1.3347023 1.3265307 1.3429248
regionEast Midlands 1.2794200 1.2617720 1.2973147
regionEast of England 1.2410285 1.2252139 1.2570475
regionNorth East England 1.4990897 1.4746657 1.5239217
regionNorth West England 1.2306656 1.2157833 1.2457310
regionScotland 1.2603066 1.2423166 1.2785562
regionSouth East England 1.2572840 1.2432673 1.2714607
regionSouth West England 1.4048482 1.3852274 1.4247485
regionWales 1.4419519 1.4179072 1.4664067
regionWest Midlands 1.0462682 1.0330972 1.0596059
regionYorkshire and the Hu 1.0937420 1.0780153 1.1096948
taillemoyen 0.7710837 0.7458044 0.7971993
taillegrand 0.6133447 0.5966207 0.6305111
Pour les prédictions, on travaille comme d’ordinaire avec les probabilités de succès (défaut du générique fitted): il faut cependant considérer que, pour obtenir les effectifs prédits, il faut multiplier cette dernière par le nombre d’essais au centre de conduite.
On peut procéder à l’interprétation des paramètres pour les examens de conduite britanniques: toute autre chose étant égale par ailleurs,
- la cote des hommes pour la réussite de l’examen est \(33\%\) plus élevée que celle des femmes;
- Londres est la région avec le plus faible taux de succès, même en prenant en compte le volume des centres de test;
- la cote pour la réussite est \(50\%\) plus élevée en Angleterre du Nord-Est (versus la référence Londres) et \(44.7\%\) plus élevée au Pays de Galles qu’à Londres, etc.
- la cote du taux de réussite est \(63\%\) plus grande dans les petits centres que dans les grands (\(1/0.614\)).
Tous les paramètres sont statistiquement significatifs. Ce n’est pas une surprise vu la taille énorme de l’échantillon.
predlogistique <- fitted(mod_gbconduite) * gbconduite$total