La régression logistique

add

\bullet Présentation:

La régression logistique, élaborée en 1944 par Joseph Berkson, permet de discriminer une variable réponse Y binaire ou polychotomique (K \geq 2 classes) à partir d’une matrice de P variables explicatives \mathbf{X} = (X ^1, \cdots, X ^P). Les données peuvent être le mixte de deux formats divers: continu et qualitatif.

La force de la régression logistique réside en la forme de la fonction lien utilisée (le logit ou le probit) et qui permet une modélisation de forme sigmoïdale incluant ainsi la notion de pente influencée par la fréquence des observations, sous forme de pondérations par secteur, lorsque nous passons d’un secteur à l’autre en fonction de la classe décrite par la réponse Y.

A ce titre, elle demeure très prisée du fait des nombreux indicateurs permettant de qualifier le modèle construit. Parmi ces indicateurs, les odds ratios ou rapports de côte en sont les plus populaires et restent très souvent utilisés en recherche clinique par exemple.

\bullet La régression logistique:

Hypothèses préliminaires: Variables explicatives \mathbf{X} continues et qualitatives et variable réponse Y qualitative. Dans le cadre de variables continues l’hypothèse de log-linéarité devrait être vérifiée.

Le modèle logit dans le cas classique Y binaire:

Le modèle le plus souvent utilisé est le modèle logit et de fonction lien,

P(Y = Cl_2) = \frac{e^{\beta_0 + \beta_1 \cdot X_i ^1 + \ldots + \beta_p \cdot X_i ^p}}{1 + e^{\beta_0 + \beta_1 \cdot X_i ^1 + \ldots + \beta_p \cdot X_i ^p}}

Qui, en fonction d’une valeur seuil fixée (0.5 en général), permet la prédiction d’un statut Y = Cl_1 ou Y = Cl_2 pour cet individu.

Il est à remarquer que le seuil de classification peut être augmenté ou diminué malgré le fait soit initialement fixé selon des critères probabilistes. Cette manipulation reste souvent obligatoire. En effet, en fonction du déséquilibre observé au sein de la répartition entre le nombre d’observations prenant pour valeur Y = Cl_1 et Y = Cl_2, le seuil de 0.5 peut s’avérer obsolète car basé sur la maximisation de la log-vraisemblance et donc sur un critère de maximisation de la classification globale.

A noter que les coefficients \beta_j élevés à l’exponentielle, e ^{\beta_j}, donnent l’odds ratio pour la j-ème variable.

De plus, lorsque nous sommes en présence de variables qualitatives, la matrice \mathbf{X} est transformée de telle manière que les modalités de chacune des variables qualitatives deviennent des variables à part entière et que celles correspondant aux modalités de référence sont supprimées de la matrice de données. Il y a alors autant de coefficients que de variables présentes dans la matrice de données.

Pour conclure le paragraphe ci-dessus, il faut garder à l’esprit que lorsque nous utilisons la régression logistique, les variables continues peuvent être insérées telles quelles dans le modèle (voit partie sur l’hypothèse de log-linéarité) contrairement aux variables qualitatives qui nécessitent de définir, pour chacune d’elles, la modalité de référence relativement à laquelle leurs odds ratios seront déterminés.

Estimation des \beta:

L’estimation des coefficients \beta = (\beta_1, \cdots, \beta_P) se fait soit par le scoring de Fisher ou l’algorithme de Newton-Raphson. Le second reste le plus populaire des deux même s’ils conduisent aux mêmes résultats.

L’algorithme de Newton-Raphson est basé sur la méthode itérative du gradient. Nous fixons l’itération initiale avec \beta^0 = (0, \ldots, 0) et posons \frac{\partial L}{\partial \beta ^k} le vecteur des dérivées partielles premières de la vraisemblance et \frac{\partial ^2 L}{\partial \beta ^k \partial (\beta ^k) ^T} la matrice des dérivées partielles secondes de la vraisemblance.

Pour \beta^k, solution courante à l’itération k, nous avons,

\beta ^{k+1} = \beta ^k - (\frac{\partial ^2 L}{\partial \beta ^k \partial (\beta ^k) ^T}) ^{-1} \times \frac{\partial L}{\partial \beta ^k}

L’algorithme converge afin d’obtenir notre vecteur d’estimateurs via les formules:

\frac{\partial L}{\partial \beta ^k} = \mathbf{X} ^T \cdot \sum_i (y_i - p(X_i; \beta ^k))

\frac{\partial ^2 L}{\partial \beta ^k \partial (\beta ^k) ^T} = -\mathbf{X} ^T \cdot \mathbf{W}_k \cdot \mathbf{X}

Avec \mathbf{W}_k la matrice carré de taille P \times P des poids où le j ^{eme} élément de la diagonale est:

\sum_{j = 1} ^P p(X_j; \beta ^k) \cdot \sum_{j = 1} ^P (1 - p(X_j; \beta ^k))

Cet algorithme se justifie car en réalité le système de solution des coefficients n’a pas de solution analytique et se détermine en appliquant une méthode de résolution numérique.

Divergence de l’algorithme et solution:

Nous avons divergence de l’algorithme de Newton-Raphson lorsque nous ne pouvons inverser la matrice \frac{\partial ^2 L}{\partial \beta ^k \partial (\beta ^k) ^T}.

La régression logistique pénalisée a été conçue dans le but de pallier à ce problème de divergence. Elle consiste en la détermination d’un paramètre \lambda qui sera  appliqué sur la diagonale de la matrice des dérivées partielles secondes de la vraisemblance afin de pouvoir l’inverser.

L’hypothèse de log-linéarité:

Seules les variables vérifiant l’hypothèse de log-linéarité peuvent être introduites légitimement dans le modèle.

L’hypothèse de log-linéarité consiste en une relation globalement linéaire entre la variable et l’odds ratio qui lui est associé. Cette hypothèse a pour conséquence qu’une augmentation de la variable X ^j de c multiplie le facteur risque par e ^{c \cdot \beta_j} quelque soit la tranche \left[ X_i ^j; X_{i+c} ^j \right].

De cette définition nous déduisons que dans le cas de variables qualitatives le problème ne se pose pas puisque la notion de pente entre les modalités est forcément induite par construction du modèle. La difficulté réside dans l’approche des variables continues puisqu’il n’existe pas vraiment de test statistique pour établir une telle vérification. Il est donc recommandé de transformer celles qui ne le sont pas en format qualitatif et de travailler sur la table de contingence associée.

Les tests d’apport d’une variable à un modèle:

Ils permettent de quantifier l’augmentation de la qualité d’un modèle lorsque nous lui insérons une variable supplémentaire. Ces tests sont notamment à l’origine des procédures sélectives pas-à-pas de type forward/backward/stepwise.

Nous comptons trois tests distincts ayant chacun leur propre particularité: le test de Wald construit par Abraham Wald, le test du rapport de vraisemblance publié par Jerzy Neyman et Egon Sharpe Pearson et le test du score de Rao conçut suite aux travaux de Calyampudi Radhakrishna Rao.

Nous présentons ci-dessous les formules de ces trois différentes tests.

– La statistique du test de Wald:

Z_j = sign(\beta_j) \cdot \frac{\beta_j}{\sigma_{\beta_j}}

, où \beta_j désigne le coefficient estimé de la variable X ^j et \sigma_{\beta_j} son écart-type.

La statistique du test suit un loi normale centrée-réduite et l’hypothèse H_0 est: « Le coefficient est nul ».

– La statistique du test de rapport de vraisemblance:

D = -2 ln (\frac{\mbox{Vraisemblance modele de reference}}{\mbox{Vraisemblance modele d'interet}})

La statistique de test suit une loi du \chi ^2 de paramètre (df_1 - df_2) où df_k désigne le nombre de variables du modèle k. L’hypothèse H_0 est: « Les modèles d’intérêt et d’origine sont similaires ».

– La statistique du test de Score de Rao:

R = U_{H_0} ^T (\beta) I_{H_0} ^{-1} (\beta) U_{H_0} (\beta)

U_{H_0} (.) est la matrice des dérivées partielles de la log-vraisemblance estimée sous contrainte de \beta_j = 0, (H_0) et I la matrice d’information de Fisher.

La statistique de test suit une loi du \chi ^2 de paramètre P. L’hypothèse H_0 est: « Les coefficients sont nuls ».

Le graphe ci-dessous, et provenant de l’oeuvre de Gilbert Saporta citée en partie bibliographie, présente de manière concrète ce que mesure chacun de ces trois tests:

ADD

Le cas Y polychotomique à K classes:

La régression logistique (binaire) de base n’est pas adaptée au cas multiclasse. Aussi, il a été développé une version étendue appelée régression logistique polytomique. L’idée est alors de convenir d’un pivot au sein de la variable réponse (généralement la classe la plus faible s’il y a notion d’ordre) et de faire M = K - 1 régressions logistiques binaires indépendantes. La formule de la log-vraisemblance à maximiser devient alors:

LL = \sum_i y_1(i) \times ln(\pi_1(i)) + \cdots + y_M(i) \times ln(\pi_M(i))

Avec i observation, Y le vecteur réponse et,

\pi_k = \frac{e^{C_k}}{1 + \sum_{l = 1} ^{k-1} e^{C_l}}

, pour,

C_l = ln(\frac{\pi_l}{\pi_M}) = a_{0,l} + \sum_{j = 1} a_{j,l} \cdot X ^j

\bullet Annexe théorique:

Nous présentons dans cette partie de l’article une esquisse de la démonstration de la formule des coefficients \beta du modèle logit.

Une façon naturelle d’expliciter le vecteur réponse Y à partir de nos variables explicatives \mathbf{X} est de chercher le vecteur des paramètres  \beta = (\beta_0, \ldots, \beta_P) qui ajuste l’équation:

P(Y = Cl_2 / X ^1 = x ^1, \ldots, X ^P = x ^P) = \beta_0 + \beta_1 \cdot X ^1 + \ldots + \beta_P \cdot X ^P

En introduisant le modèle logit(x) = \frac{e^x}{1 + e^x} nous avons,

P(Y = Cl_2 / X ^1 = x ^1, \ldots, X ^P = x ^P)

 = p(x ^1, \ldots, x ^P)

 = logit ^{-1} (\beta_0 + \beta_1 \cdot X ^1 + \ldots + \beta_P \cdot X ^P)

Et,

\beta_0 + \beta_1 \cdot X ^1 + \cdots + \beta_P \cdot X ^P = logit(p(X ^1, \cdots, x ^P)) = log(\frac{p(x ^1, \cdots, x ^P)}{1 + p(x ^1, \cdots, x ^P)})

En remarquant que la distribution suit une loi de Bernoulli, nous pouvons écrire l’estimation des paramètres:

P(Y = y_i / X ^1 = x_i ^1, \ldots, X ^P = x_i ^P) = p(x_i ^1, \ldots, x_i ^P) ^{y_i} (1 - p(x_i ^1, \ldots, x_i ^P))^{1 - y_i}

La vraisemblance s’écrit alors:

L(\beta_0, \ldots, \beta_P) = \Pi_{i = 1}^n p(x_i ^1, \ldots, x_i ^P) ^{y_i}(1 - p(x_i ^1, \ldots, x_i ^P))^{1 - y_i}

 = \Pi_{i = 1} ^n (\frac{p(x_i ^1, \ldots, x_i ^P)}{1 - p(x_i ^1, \ldots, x_i ^P)}) ^{y_i} \times (1 - p(x_i ^1, \ldots, x_i ^P))

, où x_i ^1, \ldots, x_i ^P sont les valeurs X ^1, \ldots, X ^P mesurées chez le i^{eme} observation \forall i \in \lbrace1, \ldots, n\rbrace.

Etant donné que,

\left\{ \begin{array}{ll} (\frac{p(x_i ^1, \ldots, x_i ^P)}{1 - p(x_i ^1, \ldots, x_i ^P)}) ^{y_i} = e^{y_i \times (\beta_0 + \beta_1 \cdot X ^1 + \ldots + \beta_P \cdot X ^P)} \\ 1 - p(x_i ^1, \ldots, x_i ^P) = \frac{1}{1 + e ^{(\beta_0 + \beta_1 \cdot X ^1 + \ldots + \beta_P \cdot X ^P)}} \end{array} \right.

, la vraisemblance s’écrit également sous la forme:

L(\beta_0, \ldots, \beta_p) = \Pi_i \frac{e^{y_i \times (\beta_0 + \beta_1 \cdot X ^1 + \ldots + \beta_P \cdot X ^P)}}{1 + e ^{(\beta_0 + \beta_1 \cdot X ^1 + \ldots + \beta_P \cdot X ^P)}}

Enfin, en rappelant que l(\beta_0, \ldots, \beta_n) = - log \lbrace L(\beta_0, \ldots, \beta_n) \rbrace, la maximisation de la vraisemblance s’apparente à la minimisation du - log. Nous procédons à l’annulation des dérivées de la log-Vraisemblance en dérivant la fonction L par rapport à chacun des estimateurs \beta_j, nous obtenons alors la formule générale:

\frac{\partial l}{\partial \beta_j} (\beta_0, \ldots, \beta_P) = - \sum_{i = 1} ^n x_i ^j (y_i - p(x_i ^1, \ldots, x_i ^P)), \forall j \in \left[ 1, \ldots, P \right]

Nos coefficients de régression sont calculés en résolvant le système d’équations suivant:

\left\{ \begin{array}{ll} \frac{\partial l}{\partial \beta_0} (\beta_0, \ldots, \beta_P) = 0 \\ \vdots \\ \frac{\partial l}{\partial \beta_j} (\beta_0, \ldots, \beta_P) = 0 \\ \vdots \\ \frac{\partial l}{\partial \beta_P} (\beta_0, \ldots, \beta_P) = 0 \end{array} \right.

, qui se résout par l’une des deux méthodes analytiques évoquées dans cet article: le scoring de Fisher ou l’algorithme de Newton-Raphson.

\bullet Exemple:

Soit le jeu de données ci-dessous,

add

Nous avons donc la variable réponse « Statut » et les deux variables explicatives \mathbf{X} = (X ^1, X ^2).

Nous construisons la table disjonctive de cette matrice et prenons soin d’éliminer la colonne de référence de chaque variable (ici les colonnes référant à la modalité « 1m » des deux variables X ^1, X ^2) et nous fixons comme seuil de tolérance \varepsilon = 0.0001 pour l’usage de l’algorithme de Newton-Raphson.

Nous débutons l’algorithme avec le vecteur \beta_0 = (0, 0, 0, 0, 0).

– A l’étape 1:

Nous avons trivialement \forall i, p(X_i; \beta ^0) = \frac{e ^{0 + 0 + 0 + 0 + 0 }}{ 1 + e ^0 } = \frac{1}{2}. Par conséquent, le terme,

\frac{\partial L}{\partial \mathbf{\beta}} = \begin{pmatrix} \sum_{i = 1} ^n (y_i - \frac{1}{2}) \\ \sum_{i = 1} ^n ((y_i - \frac{1}{2}) \times Z_i ^1) \\ \sum_{i = 1} ^n ((y_i - \frac{1}{2}) \times Z_i ^2) \\ \sum_{i = 1} ^n ((y_i - \frac{1}{2}) \times Z_i ^3) \\ \sum_{i = 1} ^n ((y_i - \frac{1}{2}) \times Z_i ^4) \end{pmatrix} = \begin{pmatrix} 0.0 \\ -0.5 \\ 0.5 \\ -0.5 \\ 2.5 \end{pmatrix}

Enfin, nous construisons la matrice ZZ = [ 1 \hspace*{5mm} Z ].

Nous avons alors le terme,

\frac{\partial ^2 L}{\partial \mathbf{\beta} \partial \mathbf{\beta} ^T} = \begin{pmatrix} \sum_i ^n (p_i \times (1-p_i)*ZZ_i ^1 \times ZZ_i ^1) & \ldots & \sum_i ^n (p_i \times (1-p_i)*ZZ_i ^1 \times ZZ_i ^p \\ \vdots & \ddots & \vdots \\ \sum_i ^n (p_i \times (1-p_i)*ZZ_i ^p \times ZZ_i ^1 & \ldots & \sum_i ^n (p_i \times (1-p_i)*ZZ_i ^p \times ZZ_i ^p \end{pmatrix}

 = \begin{pmatrix} 5.00 & 1.75 & 1.75 & 1.75 & 1.75 \\ 1.75 & 1.75 & 0.00 & 1.00 & 0.50 \\ 1.75 & 0.00 & 1.75 & 0.50 & 0.50 \\ 1.75 & 1.00 & 0.50 & 1.75 & 0.00 \\ 1.75 & 0.50 & 0.50 & 0.00 & 1.75 \\ \end{pmatrix}

Nous obtenons finalement,

\beta ^1 = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} - \begin{pmatrix} 5.00 & 1.75 & 1.75 & 1.75 & 1.75 \\ 1.75 & 1.75 & 0.00 & 1.00 & 0.50 \\ 1.75 & 0.00 & 1.75 & 0.50 & 0.50 \\ 1.75 & 1.00 & 0.50 & 1.75 & 0.00 \\ 1.75 & 0.50 & 0.50 & 0.00 & 1.75 \\ \end{pmatrix} \times \begin{pmatrix} 0.0 \\ -0.5 \\ 0.5 \\ -0.5 \\ 2.5 \end{pmatrix} = \begin{pmatrix} 1.6867470 \\ 0.1730559 \\ -0.7645126 \\ -1.2814896 \\ -2.9463308 \end{pmatrix}

Comme | \beta ^0 - \beta ^1 | = 6.852136 > \varepsilon nous continuons l’algorithme à l’étape 2.

Nous répétons l’algorithme jusqu’à convergence du vecteur de coefficient \beta ^k, ainsi nous obtenons:

– Étape 2: \beta ^2 = (-2.26578, -0.14463, 1.20629, 1.70543, 3.79744) et | \beta ^1 - \beta ^2 | = 2.3243 > \varepsilon, nous continuons à l’étape 3.

– Étape 3: \beta ^3 = (-2.40204, -0.11658, 1.32300, 1.79322, 3.97068) et | \beta ^2 - \beta ^3 | = 0.54205 > \varepsilon, nous continuons à l’étape 4.

– Étape 4: \beta ^4 = (-2.40809, -0.114772, 1.32827, 1.79688, 3.97752) et | \beta ^3 - \beta ^4 | = 0.023632 > \varepsilon, nous continuons à l’étape 5.

– Étape 5: \beta ^5 = (-2.40810, -0.11477, 1.32828, 1.79688, 3.97753) et | \beta ^4 - \beta ^5 | = 4.059e-05 < \varepsilon, l’algorithme a convergé et le vecteur \beta_5 est le vecteur final des coefficients estimés:

\begin{tabular}{|l|c|c|c|c|c|} \hline Parametres & Coefficients \\ \hline Intercept & -2.4081 \\ \hline X1 = 2m & -0.1148 \\ \hline X2 = 3m & 1.3283 \\ \hline X2 = 2m & 1.7969 \\ \hline X2 = 3m & 3.9775 \\ \hline \end{tabular}

Via cette forme des données, nous pouvons calculer les probabilités à risque d’affection à l’une des deux classes de la variable à expliquer « Statut ». Il s’agira de calculer pour chaque combinaison de modalités la probabilité d’être malade:

– Si l’individu est classe faible (Y = 1) pour X ^1 et pour X ^2:

P(Y_1 / X) = \frac{e ^{-2.4081052 + -0.1147692 \times 0 + 1.3282786 \times 0 + 1.7968847 \times 0 + 3.9775322 \times 0}}{1 + e ^{-2.4081052 + -0.1147692 \times 0 + 1.3282786 \times 0 + 1.7968847 \times 0 + 3.9775322 \times 0}} = 0.08255672

, soit 8\% de risque d’être Y = 1.

\vdots

– Si l’individu est classe forte (Y = 2) pour X ^1 et pour X ^2:

P(Y = Y_2 / X) = \frac{e ^{-2.4081052 + -0.1147692 \times 0 + 1.3282786 \times 1 + 1.7968847 \times 0 + 3.9775322 \times 1}}{1 + e ^{-2.4081052 + -0.1147692 \times 0 + 1.3282786 \times 1 + 1.7968847 \times 0 + 3.9775322 \times 1}} = 0.9477324

, soit 95\% de risque d’être Y_2.

Enfin, nous pouvons désormais construire la matrice de confusion suivante,

\begin{tabular}{|l|c|c|} \hline & Predit Y1 & Predit Y2 \\ \hline Y1 & 6 & 4 \\ \hline Y2 & 2 & 8 \\ \hline \end{tabular}

, à partir de laquelle nous pouvons calculer les taux de performance de la règle décisionnelle construite par régression logistique.

\bullet Application informatique:

Procédure SAS: https://support.sas.com/documentation/cdl/en/statug/63347/HTML/default/viewer.htm#statug_logistic_sect004.htm

Package et fonction R: https://stat.ethz.ch/R-manual/R-patched/library/stats/html/glm.html

\bullet Bibliographie:

– Application of the logistic function to bio-assay de Joseph Berkson

– Comprendre et utiliser les statistiques dans les sciences de la vie, Bruno Falissard

– The elements of statistical learning de Trevor Hastie, Robert Tibshirani et Jerome Friedman

– On the Problem of the Most Efficient Tests of Statistical Hypotheses de Jerzy Neyman et Egon Sharpe Pearson

– Probabilités, analyse de données et Statistique de Gilbert Saporta

– Applied Logistic Regression de David W. Hosmer et Stanley Lemeshow