Régression logistique

Exercice 3.1

Les données logistclient contiennent des données simulées pour un cas fictif de promotion pour des clients.

Estimez le modèle logistique pour la probabilité que promo=1 avec les variables explicatives nachats, sexe et tclient.

Code

library(hecmulti)
data(logistclient, package = "hecmulti")
# Modèle logistique
mod <- glm(promo ~ nachats + sexe + tclient, 
           family = binomial(link = "logit"),
           data = logistclient)
# Coefficients du modèle et statistique de Wald
summary(mod)


Call:
glm(formula = promo ~ nachats + sexe + tclient, family = binomial(link = "logit"), 
    data = logistclient)

Coefficients:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)        -0.98615    0.24634  -4.003 6.25e-05 ***
nachats             0.20775    0.03509   5.921 3.20e-09 ***
sexe               -0.26796    0.13872  -1.932   0.0534 .  
tclientoccasionnel  0.06873    0.15638   0.439   0.6603    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1384.5  on 999  degrees of freedom
Residual deviance: 1330.6  on 996  degrees of freedom
AIC: 1338.6

Number of Fisher Scoring iterations: 4

Interprétez les coefficients du modèle à l’échelle de la cote en terme de pourcentage d’augmentation ou de diminution.

Code

exp(coef(mod)) # Rapport de cote

       (Intercept)            nachats               sexe tclientoccasionnel 
         0.3730111          1.2309041          0.7649396          1.0711432

La cote pour l’offre promotionnelle (oui versus non) des hommes est 23.5% plus faible que celle des femmes, ceteris paribus
La cote des clients occasionnels est 7.11% supérieure à celle des clients fréquents, ceteris paribus. De manière équivalente, le rapport de cotes pour tclient fréquent sur occasionnel est de1/1.0711 = 0.934: les clients fréquents ont une cote 6.6% inférieure à celle des clients occasionnels, toute chose étant égale par ailleurs.
La cote de nachats augmente de 23.1% pour chaque augmentation du nombre d’achats dans le dernier mois, ceteris paribus

Testez si l’effet de nachats est statistiquement significatif à niveau $\alpha = 0.05$.

L’intervalle de confiance à 95% pour le rapport de cote de nachats, basé sur la vraisemblance profilée, est de $[1.15, 1.32]$; comme 1 est exclu, cette différence est statistiquement significative.

On obtiendrait la même conclusion avec la statistique du test de rapport de vraisemblance, ici $37.237$ pour 1 degré de liberté. La probabilité, si $\beta_{\text{nachats}}=0$, d’obtenir une telle différence d’ajustement est inférieure à $10^{-4}$, bien en deça du seuil de significativité. On rejette l’hypothèse nulle et on conclut que le nombre d’achat est important pour expliquer si une personne s’est prévalue de l’offre promotionnelle.

Code

exp(confint(mod)) #IC pour rapport de cote

                       2.5 %    97.5 %
(Intercept)        0.2292101 0.6025441
nachats            1.1500139 1.3197437
sexe               0.5825467 1.0037294
tclientoccasionnel 0.7889239 1.4568553

Code

car::Anova(mod, type = 3) # tests de rapport de vraisemblance

Analysis of Deviance Table (Type III tests)

Response: promo
        LR Chisq Df Pr(>Chisq)    
nachats   37.237  1  1.046e-09 ***
sexe       3.737  1    0.05322 .  
tclient    0.193  1    0.66011    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Choisissez un point de coupure pour la classification pour maximiser le taux de bonne classification.
1. Pour le point de coupure choisi, construisez une matrice de confusion.
2. Calculez la sensibilité, la spécificité et le taux de bonne classification manuellement. Vérifiez vos réponses avec la sortie du tableau.

Code

set.seed(60602)
# Prédictions par validation croisée 
# (moyenne de 10 réplications, K=10 plis)
pred <- hecmulti::predvc(mod)
# Extraire la variable réponse binaire 0/1
resp <- logistclient$promo

Notez que les prédictions obtenues par validation croisée sont aléatoires, donc les résultats (aire sous la courbe, valeur-$p$ du test d’adéquation, etc.) peuvent varier si vous n’utilisez pas le même germe aléatoire.

On prend le modèle ajusté avec glm et on calcule la prédiction à l’aide de la validation croisée à 10 groupes, répétée 10 fois. La fonction predvc retourne la moyenne des prédictions (ici, des probabilités) pour chacune des 1000 observations.

Code

library(ggplot2)
tableau <- perfo_logistique(prob = pred, 
                             resp = resp)
# Graphique du taux de bonne classification
# selon le point de coupure
ggplot(data = tableau, 
       aes(x = coupe, y = pcorrect)) + 
  geom_line() + 
  theme_classic() + 
  scale_y_continuous(limits = c(0, 100),
                     expand = c(0,0)) + 
  labs(x = "point de coupure",
       y = "",
       subtitle = "Taux de bonne classification")

Code

opt <- which.max(tableau$pcorrect)
knitr::kable(tableau[opt,], digits = 2)

	coupe	VP	VN	FP	FN	pcorrect	sensi	speci	fpos	fneg
45	0.45	316	281	240	163	59.7	65.97	53.93	43.17	36.71

Ensuite, il suffit de passer les valeurs de la variable réponse et nos probabilités de succès prédites aux différentes fonctions.

Si on considère des points de coupure de 0.01 à 0.99 en incréments de 0.01, on obtient un point de coupure optimal à 0.45. On note que le taux de bonne classification change assez peu au final.

Tableau 1: Matrice de confusion avec point de coupure optimal

	$Y=1$	$Y=0$
$\widehat{Y}=1$	316	240
$\widehat{Y}=0$	163	281

Ainsi, si on fait les calculs à la main, on estime

la sensibilité, $\mathsf{VP}/(\mathsf{VP} + \mathsf{FN})$, soit 316 / (316 + 163) ou 0.66.
la spécificité, $\mathsf{VN}/(\mathsf{VN} + \mathsf{FP})$, soit 281 / (281 + 240) ou 0.539.
le taux de bonne classification $\mathsf{VN} + \mathsf{VP}/(\mathsf{VN} + \mathsf{VP} + \mathsf{FN} + \mathsf{FP})$, soit (316 + 281) / 1000 ou 0.597.

Ces valeurs coincident, à arrondi près, avec ce qui est reporté dans le tableau.

Produisez un graphique de la fonction d’efficacité du récepteur (courbe ROC) et rapportez l’aire sous la courbe estimée à l’aide de la validation croisée.

Code

roc <- courbe_roc(prob = pred, resp = resp)

On obtient une estimation de l’aire sous la courbe de 0.616.

Calculez la statistique de Spiegelhalter (1986) pour la calibration du modèle. Y a-t-il des preuves de surajustement?

Code

hecmulti::calibration(prob = pred, resp = resp)

Test de calibration de Spiegelhalter (1986)
Statistique de test: 0.55 
valeur-p: 0.581

L’hypothèse nulle est que le modèle est calibré; ici, la valeur-$p$ est près de 0.5, donc on ne rejette pas l’hypothèse nulle et on conclut qu’il n’y a pas de preuve de surajustement.

Exercice 3.2

Interprétez le coefficient pour l’ordonnée à l’origine $\alpha$ en terme de pourcentage d’augmentation ou de diminution de la cote par rapport à la référence jouer à l’extérieur.

Code

data(lnh_BT, package = "hecmulti")
# Ajuster le modèle de régression logistique
mod <- glm(vainqueur ~ ., data = lnh_BT, family = binomial)
# Extraire les coefficients et les IC
expcoef <- exp(coef(mod))
ic <- confint(mod)

L’ordonnée à l’origine (premier coefficient) représente l’avantage de jouer à domicile (si $\alpha>0$), peu importe les équipes qui jouent. Le coefficient $\widehat{\alpha} = 0$ la cote pour l’équipe à domicile est 18.3 % plus élevée que celle de l’équipe en déplacement.

Calculez un intervalle de confiance de niveau 95% pour l’ordonnée à l’origine et déterminez si jouer à domicile impacte significativement le score.

L’intervalle de confiance de vraisemblance profilée est [0.1,0.3] et n’inclut pas zéro: l’effet est donc significatif. Alternativement, l’intervalle pour l’effet multiplicatif de la cote est [1.1,1.3] différent de 1.

Fournissez un tableau avec le classement des cinq premières équipes qui ont la plus grande chance de succès selon le modèle.

Il suffit de prendre les cinq équipes qui ont les plus grands coefficients de régression (en vérifiant que ces derniers sont supérieurs à zéro, le coefficient fantôme de la catégorie de référence).

Code

# Calculer les coefficients (moins ordonnée à l'origine)
# trier et garder les cinq plus grandes valeurs
classements <- sort(x = c("Anaheim_Ducks" = 0, coef(mod)[-1]),
                    decreasing = TRUE)
top5 <- head(names(classements), 5)
# Créer un tableau avec les noms
# remplacer les barres de soulignement par des espaces
knitr::kable(x = gsub("_"," ", top5), 
             col.names = "équipe",)

Tableau 2: Classement des cinq meilleurs équipes de la LNH.

équipe
Florida Panthers
Colorado Avalanche
Carolina Hurricanes
Toronto Maple Leafs
Minnesota Wild

Pour chaque match, utilisez le modèle logistique pour prédire l’équipe gagnante.
- Construisez une matrice de confusion (1 pour une victoire de l’équipe à domicile, 0 sinon) avec un point de coupure de 0.5 (assignation à l’événement ou à la classe la plus probable) et rapportez cette dernière.
- Calculez le taux de bonne classification, la sensibilité et la spécificité à partir de votre matrice de confusion.

Code

set.seed(60602)
# Prédire par validation croisée avec point de coupure 0.5
pred <- hecmulti::predvc(modele = mod)
classif <- pred > 0.5
# Créer une matrice de confusion TRUE=1, FALSE=0
table(prediction = as.integer(classif), 
      vainqueur = as.integer(lnh_BT$vainqueur))

          vainqueur
prediction   0   1
         0 340 218
         1 268 486

Le taux de bonne classification est de 63%, la sensibilité de 69% et la spécificité de 55.9%

Produisez un graphique de la fonction d’efficacité du récepteur et rapportez l’aire sous la courbe. Commentez sur la qualité prédictive globale du modèle.

Code

roc <- hecmulti::courbe_roc(prob = pred, resp = lnh_BT$vainqueur)

L’aire sous la courbe est de 0.66, soit plus qu’une assignation aléatoire de 0.5. Le modèle a un faible pouvoir prédictif, mais ça peut suffire pour des paris sportifs si le modèle bat la concurrence.

Exercice 3.3

On s’intéresse à la satisfaction de clients par rapport à un produit. Cette dernière est mesurée à l’aide d’une échelle de Likert, allant de très insatisfait (1) à très satisfait (5). Les 1000 observations se trouvent dans la base de données multinom du paquet hecmulti.

Modélisez la satisfaction des clients en fonction de l’âge, du niveau d’éducation, du sexe et du niveau de revenu.

Est-ce que le modèle de régression multinomiale ordinale à cote proportionnelles est une simplification adéquate du modèle de régression multinomiale logistique? Si oui, utilisez ce modèle pour la suite. Sinon, ajustez le modèle de régression multinomiale logistique avec 1 comme catégorie de référence, 1 pour revenu et sec pour éducation¹ et utilisez ce dernier pour répondre aux autres questions.

Les niveaux des facteurs non-ordonnés (catégories) dans R sont classés en ordre alphanumérique. Il faut donc modifier la catégorie de référence uniquement pour le niveau d’éducation avant d’ajuster le modèle.

Code

data(multinom, package = "hecmulti")
db <- multinom |>
  dplyr::mutate(educ = relevel(educ, ref = "sec"),
                y = ordered(y))
mod1 <- nnet::multinom(
  y ~ sexe + educ + revenu + age,
  data = db,
  trace = FALSE)

Le modèle a 28 coefficients, dont deux pour éducation et revenu et un pour sexe et age, par niveau.

Pour ajuster le modèle à cotes proportionnelles, il faut d’abord convertir la variable réponse en variable ordinale à l’aide de ordered si ce n’est pas déjà la classe de la variable. Cela permettra de spécifier l’ordre des modalités.

Code

# Ajuster modèle à cote proportionnelle
mod0 <- MASS::polr(
  ordered(y) ~ sexe + educ + revenu + age,
  data = db)
# Calculer statistique de test
# (rapport de vraisemblance)
stat <- deviance(mod0) - deviance(mod1)
# À comparer à une loi khi-deux avec 
# npar1-npar0 degrés de liberté
npar0 <- length(coef(mod0)) + length(mod0$zeta)
npar1 <- length(coef(mod1))
# Calcul de la valeur-p
# Probabilité que khi-deux (df) excède 'stat'
pchisq(stat, 
       df = npar1 - npar0, 
       lower.tail = FALSE)

[1] 0.0008817018

Le modèle ordinal a 10 paramètres, contre 28 pour le modèle multinomial logistique. On peut faire un test du rapport de vraisemblance en comparant la différence des log-vraisemblance des deux modèles emboîtés: la valeur de la statistique est 42.702. La valeur-$p$ estimée est $8.8 \times 10^{-4}$$, donc on rejette l’hypothèse nulle et on conclut que le modèle à cote proportionnelle n’est pas adéquat.

Interprétez l’effet des variables éducation et sexe pour la catégorie 2 (par rapport à 1).

Il suffit de regarder les coefficients $\exp(\widehat{\beta}_{\text{sexe}}), \ldots$ associés et les interpréter en termes de rapport de cote, pour une régression logistique ordinaire.

Code

# Coefficient correspondants à sexe et éducation
exp(coef(mod1)["2", 2:4])

     sexe educcegep   educuni 
 1.424053  1.104006  1.316095

La cote pour les femmes pour insatisfait par rapport à très insatisfait est 42.4% plus élevée que pour les hommes, toute chose étant égale par ailleurs.
La cote pour les individus qui ont un diplôme collégial pour insatisfait par rapport à très insatisfait est 10.4% plus élevée que pour ceux qui on un diplôme secondaire, toute chose étant égale par ailleurs.
La cote pour les individus qui ont un diplôme universitaire pour insatisfait par rapport à très insatisfait est 31.6% plus élevée que pour ceux qui on un diplôme secondaire, toute chose étant égale par ailleurs.

Est-ce que le modèle avec une probabilité constante pour chaque item est adéquat lorsque comparé au modèle qui inclut toutes les covariables?

Pour répondre à cette question, on ajuste le modèle multinomial logistique avec uniquement une constante. Les probabilités prédites sont simplement la proportion empirique des observations de l’échantillon: on peut ainsi vérifier que le modèle a convergé en comparant les prédictions et ces proportions. Une petite différence numérique est possible puisque le modèle multinomial logistique est ajusté à l’aide d’une procédure d’optimisation numérique itérative.

Code

# Ajuster modèle avec proba constante
mod0 <- nnet::multinom(
  y ~ 1, 
  # ordonnée à l'origine seulement
  data = db,
  trace = FALSE)
# Vérifier convergence
pred0 <- predict(
  object = mod0,
  # Entrer une bd avec une ligne
  newdata = db[1,], 
  type = "prob")
# Calculer la proportion de chaque categ.
proportions_y <- table(db$y)/nrow(db)
# Calculer différences
pred0 - proportions_y


            1             2             3             4             5 
 2.278825e-05  1.157814e-05 -6.983907e-06 -5.040094e-06 -2.234239e-05

On peut ensuite calculer la statistique de rapport de vraisemblance en comparant les déviances des deux modèles.

Code

# Calculer statistique de test
# (rapport de vraisemblance)
stat <- deviance(mod0) - deviance(mod1)
# À comparer à une loi khi-deux avec 
# npar1-npar0 degrés de liberté
npar0 <- length(coef(mod0)) + length(mod0$zeta)
npar1 <- length(coef(mod1))
# Calcul de la valeur-p
# Probabilité que khi-deux (df) excède 'stat'
pchisq(stat, 
       df = npar1 - npar0, 
       lower.tail = FALSE)

[1] 0.0002755707

La statistique pour le test du rapport de vraisemblance que tous les coefficients associés aux covariables sont nuls (24 paramètres supplémentaires) est 55.41, et si le modèle sans covariable était vrai, cette statistique serait approximativement $\chi^2_{24}$. La valeur-$p$ arrondie est 0.0003, on rejette l’hypothèse nulle que tous les coefficients associés aux variables explicatives sont nuls. On conclut qu’au moins une covariable est utile pour prédire une cote par rapport au modèle avec une probabilité constante.

Est-ce que l’effet de la variable âge est globalement significatif?

Puisqu’on modélise quatre rapport de cotes à l’aide d’un modèle logistique, d’où $\beta_{\texttt{age}_2}=\beta_{\texttt{age}_3}=\beta_{\texttt{age}_4}=\beta_{\texttt{age}_5}=0.$ On peut obtenir la valeur-$p$ avec le tableau d’analyse de déviance, qui rapporte la valeur du test de rapport de vraisemblance. On conclut que l’âge impacte la probabilité des différents items de satisfaction.

Code

car::Anova(mod1, type = 3)

Analysis of Deviance Table (Type III tests)

Response: y
       LR Chisq Df Pr(>Chisq)    
sexe    14.6448  4  0.0054975 ** 
educ     9.7083  8  0.2860960    
revenu   7.8239  8  0.4508547    
age     23.1285  4  0.0001194 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Fournissez un intervalle de confiance à niveau 95% pour l’effet multiplicatif d’une augmentation d’une unité de la variable âge pour chacune des cote par rapport à très insatisfait (1). Que concluez-vous sur l’effet de âge pour les réponses 2 à 5 par rapport à 1?

Les intervalles de confiance sont obtenus en prenant l’exponentielle des intervalles de confiance profilée pour les coefficients associés à age:

Code

exp(confint(mod1)["age",,])

               2         3         4         5
2.5 %  0.9745999 0.8723827 0.9415834 0.9879372
97.5 % 1.0108252 0.9640311 0.9893561 1.0205983

$[0.975; 1.011]$ pour $\beta_{\texttt{age}_{2|1}}$ (pas significatif),
$[0.871; 0.963]$ pour $\beta_{\texttt{age}_{3|1}}$ (significatif),
$[0.941; 0.989]$ pour $\beta_{\texttt{age}_{4|1}}$ (significatif) et
$[0.988; 1.021]$ pour $\beta_{\texttt{age}_{5|1}}$ (pas significatif).

Écrivez l’équation de la cote ajustée pour satisfait (4) par rapport à très insatisfait (1).

Pour obtenir l’équation ajustée, on utilise uniquement les coefficients pour $Y=4$ dans le tableau des coefficients.

Code

round(coef(mod1)["4",], 3)

(Intercept)        sexe   educcegep     educuni     revenu2     revenu3 
     -0.502       0.613       0.676       0.800      -0.038       0.009 
        age 
     -0.035

\[\begin{align*} \frac{\Pr(Y=4 \mid\boldsymbol{X})}{\Pr(Y=1 \mid \boldsymbol{X})} & =\exp(-0.502 - 0.035\texttt{age} + 0.676 \texttt{cegep} + 0.8\texttt{uni} \\&-0.038 \texttt{revenu}_2 +0.009 \texttt{revenu}_3 + 0.613 \texttt{sexe}) \end{align*}\]

Prédisez la probabilité qu’un homme de 30 ans qui a un diplôme collégial et qui fait partie de la classe moyenne sélectionne une catégorie donnée. Quelle modalité est la plus susceptible?

Code

# Profil du client
profil <- data.frame(
  sexe = 0, 
  age = 30, 
  educ = "cegep", 
  revenu = "2")
# Probabilité du score de satisfication
predict(mod1, 
        newdata = profil, 
        type = "prob")

         1          2          3          4          5 
0.32094804 0.22162733 0.03852415 0.12699328 0.29190719

Code

# Modalité la plus susceptible
predict(mod1, 
        newdata = profil)

[1] 1
Levels: 1 2 3 4 5

Les probabilités prédites sont (0.321; 0.222; 0.039; 0.127; 0.292). La modalité la plus susceptible est donc très insatisfait (1).

Notes de bas de page

Utilisez la fonction relevel pour changer la catégorie de référence, avec relevel(educ, ref = 'sec').↩︎

	\(Y=1\)	\(Y=0\)
\(\widehat{Y}=1\)	316	240
\(\widehat{Y}=0\)	163	281