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.

data(evaluation, package = "hecmulti")
mod0 <- glm(choix ~ force, 
           data = evaluation, 
           family = binomial)

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 cote
Waiting 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 vraisemblance
Analysis 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

Références

Liu, Q., Häubl, G., & Castelo, N. (2026). Consumers with weaker applications are less receptive to algorithmic evaluation. Journal of Consumer Research, ucag033. https://doi.org/10.1093/jcr/ucag033