C'est de la data et mon expeRtise afin d'en tirer le maximum

Calculer et rapporter les odds ratio sans se fatiguer

Très récemment, j’ai travaillé sur la construction d’un score prédictif dans le domaine médical. Il s’agissait de construire un outil d’aide au pronostic de mauvaise réponse à un traitement donné, à l’aide d’une dizaine de variables prédictives.

La méthodologie classiquement employée consiste, à :

  1. Catégoriser chaque variable prédictive numérique en une variable catégorielle de quelques modalités,
  2.  Utiliser des régressions logistiques univariées pour étudier l’association entre chaque variable prédictive et la variable à expliquer (ici la réponse au traitement : mauvaise ou bonne).
  3. Sélectionner les variables prédictives qui présentent une pvalue < 0.2 (d’autres seuils sont possibles, comme 0.15, ou 0.10).
  4. Construire un modèle de régression logistique multivarié contenant l’ensemble des variables sélectionnées à partir des régressions univariées.
  5. Utiliser une méthode de sélection des variables (par exemple en pas à pas descendant) pour aboutir à un modèle final, plus simple, ne contenant que des variables prédictives dont l’association avec la variable réponse est significative, ou proche de la significativité.
  6.  Construire le score en arrondissant les odds ratio ajustés du modèle logistique multivarié final.

Au bout du compte, cette méthodologie nécessite de faire tourner de nombreux modèles de régression. Et comme on est dans le cadre de régressions logistiques (simples ou multivariées)les associations (simples ou ajustées) entre les variables prédictives et la variable réponse, se mesurent à l’aide d’un odds ratio. Cet odds ratio est obtenu par l’exponentielle du coefficient, fourni en sorti de la régression logistique.

Rapporter tous les odds ratio, c'est fastidieux !

Sur le principe calculer les odds ratio n’est pas particulièrement difficile. Ce qui est très chronophage, et profondément ennuyant, c’est de synthétiser tous les résultats obtenus dans des tables. En effet, pour un rapport, ou une publication, il est nécessaire de présenter les odds ratio, leur intervalle de confiance, et la p-value, obtenus à chaque étape,
c’est à dire :

  • avec les régressions logistiques univariée (nommés aussi odds ratio bruts)
  • avec la régression logistique multivarié du modèle complet (nommés
    aussi appelés odds ratio ajustés),
  • avec la régression logistique multivarié du modèle final.

C’est vraiment pénible. J’avais un demi pied dans la dépression 😉 quand j’ai découvert le package finalfit !

Ce package est absolument GENIAL ! Il fait quasiment tout le travail à votre place ! Et en plus, il est tidyverse compatible !

Utiliser le package finalfit pour obtenir et présenter les odds ratio

Pour illustrer l’utilisation du package finalfit, je vais employer les données heart_disease du package funModelling:

library(funModeling)
    data("heart_disease")
    HD <- heart_disease
    head(HD)

    ##   age gender chest_pain resting_blood_pressure serum_cholestoral
    ## 1  63   male          1                    145               233
    ## 2  67   male          4                    160               286
    ## 3  67   male          4                    120               229
    ## 4  37   male          3                    130               250
    ## 5  41 female          2                    130               204
    ## 6  56   male          2                    120               236
    ##   fasting_blood_sugar resting_electro max_heart_rate exer_angina oldpeak
    ## 1                   1               2            150           0     2.3
    ## 2                   0               2            108           1     1.5
    ## 3                   0               2            129           1     2.6
    ## 4                   0               0            187           0     3.5
    ## 5                   0               2            172           0     1.4
    ## 6                   0               0            178           0     0.8
    ##   slope num_vessels_flour thal heart_disease_severity exter_angina
    ## 1     3                 0    6                      0            0
    ## 2     2                 3    3                      2            1
    ## 3     2                 2    7                      1            1
    ## 4     3                 0    3                      0            0
    ## 5     1                 0    3                      0            0
    ## 6     1                 0    3                      0            0
    ##   has_heart_disease
    ## 1                no
    ## 2               yes
    ## 3               yes
    ## 4                no
    ## 5                no
    ## 6                no 

Remarque : Par commodité, je ne vais pas stratifier les variables numériques, mais les utiliser comme telles (c’est à dire comme des variables numériques).

Les odds ratio des analyses univariées

Pour obtenir les odds ratios des analyses logistiques univariées, il est nécessaire de définir :

  • la variable dépendante (ici `has_heart_disease`)
  • les variables prédictives (appelées ici “explanatory”)

Toutes les régressions logistiques univariées sont ensuite réalisées par la fonction `glmuni()`, et les résultats collectés et mis en forme par la fonction `fit2df()`

library(tidyverse)
library(finalfit)

dependent = "has_heart_disease"
explanatory = c("age","gender","chest_pain","resting_electro","fasting_blood_sugar", "max_heart_rate")


res_glm_uni <- HD%>%
    glmuni(dependent, explanatory) %>% 
    fit2df(estimate_suffix=" (univarié)")

res_glm_uni

    ##            explanatory               OR (univarié)
    ## 1                  age   1.05 (1.03-1.08, p<0.001)
    ## 2           gendermale   3.57 (2.12-6.16, p<0.001)
    ## 3          chest_pain2   0.50 (0.16-1.61, p=0.237)
    ## 4          chest_pain3   0.61 (0.22-1.77, p=0.339)
    ## 5          chest_pain4  6.15 (2.43-17.08, p<0.001)
    ## 6     resting_electro1 5.09 (0.63-104.24, p=0.163)
    ## 7     resting_electro2   2.00 (1.26-3.18, p=0.003)
    ## 8 fasting_blood_sugar1   1.15 (0.61-2.18, p=0.660)
    ## 9       max_heart_rate   0.96 (0.94-0.97, p<0.001) 

Pour obtenir un rendu élégant dans un document word ou pdf, on peut utiliser les lignes suivantes :

kable(res_glm_uni,row.names=FALSE, align=c("l", "l", "r", "r", "r", "r"))
 

Voici le rendu sous Word :

Remarque : pour plus d’informations sur la fonction kable(), vous pouvez consulter l’article Guide de démarrage en R markdown.

Les odds ratio des modèles univariés et multivarié complet

Pour cela, on utilise cette fois la fonction finalfit() :


res_glm_uni_multi <- HD %>%
    finalfit(dependent, explanatory)
 
res_glm_uni_multi
 
    ##    Dependent: has_heart_disease                     no          yes
    ## 1                           age Mean (SD)   52.6 (9.5)   56.6 (7.9)
    ## 8                        gender    female    72 (43.9)    25 (18.0)
    ## 9                                    male    92 (56.1)   114 (82.0)
    ## 2                    chest_pain         1     16 (9.8)      7 (5.0)
    ## 3                                       2    41 (25.0)      9 (6.5)
    ## 4                                       3    68 (41.5)    18 (12.9)
    ## 5                                       4    39 (23.8)   105 (75.5)
    ## 11              resting_electro         0    95 (57.9)    56 (40.3)
    ## 12                                      1      1 (0.6)      3 (2.2)
    ## 13                                      2    68 (41.5)    80 (57.6)
    ## 6           fasting_blood_sugar         0   141 (86.0)   117 (84.2)
    ## 7                                       1    23 (14.0)    22 (15.8)
    ## 10               max_heart_rate Mean (SD) 158.4 (19.2) 139.3 (22.6)
    ##               OR (univariable)         OR (multivariable)
    ## 1    1.05 (1.03-1.08, p<0.001)  1.04 (1.00-1.08, p=0.047)
    ## 8                            -                          -
    ## 9    3.57 (2.12-6.16, p<0.001) 5.77 (2.93-11.90, p<0.001)
    ## 2                            -                          -
    ## 3    0.50 (0.16-1.61, p=0.237)  0.98 (0.27-3.64, p=0.975)
    ## 4    0.61 (0.22-1.77, p=0.339)  0.93 (0.30-3.04, p=0.901)
    ## 5   6.15 (2.43-17.08, p<0.001) 8.00 (2.75-25.56, p<0.001)
    ## 11                           -                          -
    ## 12 5.09 (0.63-104.24, p=0.163) 2.94 (0.29-68.83, p=0.400)
    ## 13   2.00 (1.26-3.18, p=0.003)  1.76 (0.96-3.23, p=0.067)
    ## 6                            -                          -
    ## 7    1.15 (0.61-2.18, p=0.660)  1.22 (0.54-2.76, p=0.625)
    ## 10   0.96 (0.94-0.97, p<0.001)  0.97 (0.95-0.99, p<0.001) 

Comme vous pouvez le voir la fonction finalfit() renvoie aussi les moyennes et écart type des variables numériques dans les deux groupes de la variable dépendante, ainsi que les pourcentages lorsque les variables sont catégorielles.

Et pour un rendu élégant au format word,et/ou pdf :

kable(res_glm_uni_multi,row.names=FALSE, align=c("l", "l", "r", "r", "r", "r"))
 

Si vous souhaitez ne conserver que les résultats multivariés, vous pouvez simplement supprimer la colonne des résultats univariés, en utilisant la fonction select() du package dplyr :

res_glm_multi1 <- res_glm_uni_multi %>% 
    dplyr::select(-"OR (univariable)") 

Pour plus d’information sur l’utilisation de la fonction select() et du package dplyr, vous pouvez consulter l’article : Initiation à la manipulation de données avec le package dplyr.

 

Si vous n’êtes pas intéressé par la partie descriptive, vous pouvez utiliser la fonction glmulti() à la place de la fonction finalfit() :

res_glm_multi2  <- HD%>%
    glmmulti(dependent, explanatory) %>% 
    fit2df(estimate_suffix="(multivarié)")
res_glm_multi2  

    ##            explanatory             OR(multivarié)
    ## 1                  age  1.04 (1.00-1.08, p=0.047)
    ## 2           gendermale 5.77 (2.93-11.90, p<0.001)
    ## 3          chest_pain2  0.98 (0.27-3.64, p=0.975)
    ## 4          chest_pain3  0.93 (0.30-3.04, p=0.901)
    ## 5          chest_pain4 8.00 (2.75-25.56, p<0.001)
    ## 6     resting_electro1 2.94 (0.29-68.83, p=0.400)
    ## 7     resting_electro2  1.76 (0.96-3.23, p=0.067)
    ## 8 fasting_blood_sugar1  1.22 (0.54-2.76, p=0.625)
    ## 9       max_heart_rate  0.97 (0.95-0.99, p<0.001) 

Les odds ratio de tous les modèles

Il s’agit ici d’obtenir et de présenter les odds ratio :

  • des modèles de régression univariés,
  • du modèle de régression multivarié complet,
  • d’un modèle de régression multicarié restreint (obtenu par sélection des variables).

Attention, la sélection des variables du modèles multivarié doit être réalisé par ailleurs. Pour donner un exemple, je le fais ici par une procédure descendante basée sur le score d’Akaïké (AIC), avec la fonction `stepAIC()` du package `MASS`.

glm_full <- glm(has_heart_disease ~age + gender + chest_pain + resting_electro + fasting_blood_sugar + max_heart_rate, 
    data=HD, family=binomial )

library(MASS)
glm_multi_restreint <- stepAIC(glm_full, trace=FALSE) 
library(car) 
Anova(glm_multi_restreint) 
## Analysis of Deviance Table (Type II tests) 
## 
## Response: has_heart_disease 
##              LR Chisq    Df Pr(>Chisq)    
## age               5.263  1    0.02179 *  
## gender           27.142  1  1.891e-07 ***
## chest_pain       54.337  3  9.510e-12 ***
## max_heart_rate   15.361  1  8.882e-05 *** 

Les variables incluses dans le modèles multivarié final sont :

  • age
  • gendre
  • chest_pain
  • max_heart_rate

Pour compiler les résultats des 3 niveaux d’analyse (modèles univariés, modèle multivarié complet, et modèle multivarié final), la démarche consiste à :

  1. Construire une table descriptive avec la fonction `summary_factorlist()`
  2. Construire une table avec les OR des modèles univariés avec la fonction `glmuni()`,
  3. Construire une table avec les OR du modèle multivarié complet avec la fonction `glmmulti()`,
  4. Construire une table avec les OR du modèle multivarié final avec la fonction `glmmulti()`, en utilisant un subset des variables prédictives (ici explanatory_final),
  5. Combiner toutes les tables en utilisant la fonction `finalfit_merge()`.

dependent = "has_heart_disease"
explanatory_full = c("age","gender","chest_pain","resting_electro","fasting_blood_sugar", "max_heart_rate")
explanatory_final = c("age","gender","chest_pain", "max_heart_rate")

res_summary <- HD %>% 
     summary_factorlist(dependent, explanatory_full, fit_id=TRUE) 

res_uni <- HD%>%
    glmuni(dependent, explanatory_full) %>%
    fit2df(estimate_suffix="(univarié)") 

res_multi_full <- HD%>%
    glmmulti(dependent, explanatory_full) %>%
    fit2df(estimate_suffix="(ajustés - modèle complet)") 

res_multi_final <- HD%>%
    glmmulti(dependent, explanatory_final) %>%
    fit2df(estimate_suffix="(ajustés - modèle final)") 

tab_res <- res_summary %>% 
              finalfit_merge(res_uni) %>%
              finalfit_merge(res_multi_full) %>% 
              finalfit_merge(res_multi_final) %>% 
              dplyr::select(-levels, fit_id, -index) # supression des colonnes qui ne sont pas nécessaires 
      
tab_res

tab_res
knitr::kable(tab_res, row.names=FALSE, align=c("l", "l", "r", "r", "r", "r", "r")) 

Voici le rendu sous Word :

Visualiser les odds ratio

Pour cela, il suffit d’employer la fonction or_plot()

HD %>%
    or_plot(dependent, explanatory) 

Remarques

Les modèles de régression logistique construits par les fonctions du package finalfit utilisent une distribution binomiale des erreurs. Avant d’utiliser ces fonctions, il est nécessaire de vérifier qu’aucune surdispersion n’est présente. En cas de surdispersion, les intervalles de confiance et les p-values ne sont pas valides.

Pour aller plus loin

Le package finalfit n’est pas spécifiquement dédié à la régression logistique ! Il contient des fonctions similaires à celles présentées ici et qui peuvent être employés dans le cadre de régressions linéaires simples et multiples, dans le cadre de modèles à effets mixtes ou encoredans le cadre d’analyse de survie. Je vous encourage vivement à découvrir les possibilités de ce package en consultant sa vignette.

Et si vous souhaitez une bonne introduction, pratique, à la régression logistique dans un contexte médical / épidémiologique, je vous conseille cette vidéo de Tierry Ancelle,

ainsi que les deux articles suivants :

Alors, que pensez vous des possibilités de ce package finalfit ? Dites moi, en commentaire, quelle fonction vous parait la plus bluffante ?

Si cet article vous a plu, ou vous a été utile, et si vous le souhaitez, vous pouvez soutenir ce blog en faisant un don sur sa page Tipeee 🙏

Crédits photos :Image par photosforyou de Pixabay

18 réponses

  1. Merci claire pour ce riche article. Le package est vraiment génial. Merci également pour tout ce que tu fais pour nous des experts de la statistique et du logiciel R. Si tu pouvais nous aider avec un article sur l’analyse de survie.
    1 . Conditions de réalisation
    2. Définition des variables
    3. Analyse et Interprétation

  2. Merci beaucoup pour ce tutoriel très claire et très intéressant!
    J’aurai voulu savoir si avec ce package il serait possible de faire la même chose avec les risques relatifs plutôt que les OR et en prenant en compte les poids de sondage ? J’ai essayé de trouver un moyen de le combiner avec le package survey mais n’ai pas réussi :/
    Très bonne journée à vous !
    Et merci encore pour tout 🙂

  3. Merci pour cet article.
    J’avais justement un rapport d’études statistiques à présenter et jumeler les coefficients, les odds ratios, les p-value sur un même tableau. J’ai dû tout programmer tout ça moi-même. J’ignorais qu’un package existait pour le faire. Vaut mieux tard que jamais 😁😁.

  4. Bonjour,
    Merci pour ce tutoriel très clair.
    Je bloque néanmoins à l’étape
    res_glm_uni_multi %
    finalfit(dependent, explanatory)
    Ou lorsque j’applique à mon dataset j’obtiens l’erreur:

    Error: Tibble columns must have consistent lengths, only values of length one are recycled:
    * Length 82: Column `p`
    * Length 83: Columns `explanatory`, `estimate`, `confint[, 1]`, `confint[, 2]`

    Si vous avez un conseil, je suis preneur.

    Bien à vous

  5. Bonjour,
    Merci pour cette fonction qui simplifie la vie !
    Quelle test statistique utilisez vous dans l’analyse univariée pour les variables continues ? Il s’agit bien d’un test de Wilcoxon si les variables ne sont pas distribuées normalement ?

    1. Bonjour,

      je ne suis pas sûre de comprendre votre question.
      Le test de Wilcoxon permet de comparer les médianes de deux groupes. On le présente généralement comme une alternative non paramétrique au test de comparaison de deux moyennes, de Student.
      Bonne continuation

  6. Merci c’est très interessant.
    J’ai un blocage à cette étape:
    res_glm_uni %
    glmuni(dependent,explanatory) %>%
    fit2df(estimate_suffix=”univarié”)
    Error in eval(lhs, parent, parent) : objet ‘HD’ introuvable

    Quelqu’un peut-il m’aider?
    Merci beaucoup

    1. Bonjour,

      vous avez oublié de définir HD:
      library(funModeling)
      data (heart_disease)
      HD <- heart_disease

      Bonne continuation

  7. Bonjour,
    comment faire pour télécharger le tableau obtenu avec knitr::kable pour l’avoir au format word et / ou pdf pour l’intégrer dans mon outil de rédaction ?

    Bien à vous.

    1. Bonjour,

      en knittant le document (en cliquant sur le bouton avec pelote) et en choisissant un rendu word, vouas allez obtenir une document word.
      J’espère que ça vous aide…

  8. Merci pour cet article… et pour tant d’autres !

    J’ai un vague souvenir comme quoi quand une cellule d’une table de contingence contient moins de cinq cas, il faut utiliser le Fisher Exact Test pour obtenir l’OR.

    Pouvons-nous faire confiance aux OR donnés par la régression logistique s’il existe des cas de valeurs de variables catégorielles avec moins de cinq cas ? Dans votre exemple, la variable resting_electro a deux cases avec n < 5.

    Merci !

    1. Bonjour Antoine,

      je pense que ce sont les effectifs théoriques (et pas les effectifs observés) qui doivent être >5.
      Lorsque ces effectifs théoriques sont <5, effectivement on utilise le test exact de Fisher pour obtenir l’OR.

      Pour la régression logistique, à ma connaissance, il est seulement recommandé d’avoir au moins 10 fois plus d’événements que de paramètres dans le modèle. Mais je trouve votre question très pertinente…mais je n’ai pas la réponse.

Laisser un commentaire

Votre adresse de messagerie ne sera pas publiée. Les champs obligatoires sont indiqués avec *

Bonjour !

vous venez souvent ?

Identifiez-vous pour avoir accès à toutes les fontionnalités !

Aide mémoire off'R ;)

Enregistrez vous pour recevoir gratuitement mes fiches “aide mémoire” (ou cheat sheets) qui vous permettront de réaliser facilement les principales analyses biostatistiques avec le logiciel R et pour être informés des mises à jour du site.