L’exemple traité est un modèle de régression logistique avec effet
aléatoire. La démarche de l’analyse reste identique quel que soit le
modèle. Seul l’écriture du modèle diffère.
Objectif du document :
Il s’agit expliquer la démarche de construction d’un modèle
explicatif pour une variable binaire (1/0) notée Y et nommée variable à
expliquer en fonction de différentes variables explicatives et de
l’appliquer sous le logiciel R. Des effets aléatoires sont nécessaires
dans un modèle lorsque les données sont hiérarchisées (groupes
d’observations non indépendants) : plusieurs mesures des individus
que l’on mesure plusieurs fois au cours du temps…On est généralement
intéressé à comparer les différentes mesures utilisées comme répétitions
pour une expérience. L’effet aléatoire permet d’inférer sur l’effet
recherché
Plan du document :
I – Le modèle de régression logistique : quelques généralités
II – Etapes d’analyse
III – Exemple d’écriture d’un modèle glmer
I – Le modèle de régression logistique : quelques généralités
Voir Procédure régression logistique sans effet aléatoire Partie1 =>l’emoji⚠️ vous informera des changements à prendre en compte pour le modèle Glmer
Rappel des étapes à suivre
II – Etapes d’analyses
Etape
Préalable – Mise en forme de la BDD et création des variables
nécessaires pour le modèle et vérification des effectifs dans chaque
modalité.
Etape 1 - Réalisation de modèles de régression logistique univariés pour toutes les variables explicatives . (⚠️ les OR pourront être obtenus avec le package(broom.mixed) et la fonction tidy)
Etape 2 - Vérification de la linéarité des variables quantitatives.
Etape 3 - Test de l’association entre les variables explicatives : calcul des VIF
Etape 4 - Mise en œuvre de l’analyse multivariée
Etape 5 - Validité du modèle avec la vérification des résidus (⚠️besoin d’un nouveau package DHARMa)
Etape 6
– Présentations et interprétations des résultats. Plusieurs
présentations des résultats sur les facteurs de risque ayant un effet
sur la pratique du sport sont possibles : tableaux, graphiques,…
Etape 7 - Visualisation graphique des effets du modèle (⚠️avec le package ggeffects)
Etape 8 – Calcul d’un pseudo R² pour les modèles glmer (⚠️besoin d’un nouveau package (Munin) pour avoir la fonction r.squaredGLMM)
Etape 9 – Vérification de la distribution des effets aléatoires
Pour les modèles allant jusqu’à l’analyse multivariée, il est
indispensable de travailler sur une BDD complète. Les variables
d’intérêts : à expliquer et explicatives ne doivent pas posséder de
données manquantes (DM). Si tel est le cas, pensez à créer un sous jeu
de données complet avant la réalisation d’un modèle. Vérification avec
la package DataExplorer et la fonction plot_missing(BDD).
Dans cet exemple, nous allons explorer l’exemple 2 sur le cancer du
poumon à l’aide d’un ensemble de données simulées que nous avons mis en
ligne. Divers résultats ont été recueillis sur les patients, qui sont
imbriqués dans les médecins, lesquels sont à leur tour imbriqués dans
les hôpitaux. Il existe également quelques variables au niveau du
médecin, telles que l’expérience, que nous utiliserons dans notre
exemple.
Pour aller plus loin: https://stats.oarc.ucla.edu/r/dae/mixed-effects-logistic-regression/
packages
Code
library(tidyverse)#Pour les modèles glm avec effet aléatoire : library(lme4)#Automatisation de l'exploration des données (dimension, données manquantes, distribution, ...) library (DataExplorer) # pour organiser les données library(tidyverse) # pour télécharger la BDD dans le projet library(here) # pour les graphes library(ggplot2) library(rstatix) library(ggpubr) # Distribution variable continue library(FSA) # pour spécifier le package de la fonction si conflit entre 2 packages library(conflicted) conflict_prefer("filter", "dplyr") conflict_prefer("select", "dplyr") #significativité globale de la variable dans le modèle testé library(car)# collinéarité entre variables explicatives library(performance) #Au choix pour la sortie du résultat du modèle tableau et/ou graphe library(gtsummary) library(finalfit) #Permet de récupérer les OR library(broom.mixed) #Pour les comparaison 2 à 2 des OR : library(emmeans)#Pour les résidus library(DHARMa) #visualisation des effets du modèle library(ggeffects) #Pour le pseudo R2 library(MuMIn) # Si besoin d'augmenter l'itération : glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 100000))
Importation d’une BDD pour tester le Glmer
Code
# lien : https://stats.oarc.ucla.edu/r/dae/mixed-effects-logistic-regression/hdp <-read.csv("https://stats.idre.ucla.edu/stat/data/hdp.csv")hdp <-within(hdp, { Married <-factor(Married, levels =0:1, labels =c("no", "yes")) DID <-factor(DID) HID <-factor(HID) CancerStage <-factor(CancerStage)})
Sous jeu de données pour l’analyse
Code
# Réduction de cette BDD hiérarchiséeBDD<-hdp |>select(DID,remission,IL6, CRP, CancerStage , LengthofStay , Experience )
Nous utilisons ci-dessous la commande glmer pour estimer un modèle de
régression logistique à effets mixtes avec Il6, CRP, et LengthofStay
comme prédicteurs continus au niveau du patient, CancerStage comme
prédicteur catégorique au niveau du patient (I, II, III, ou IV),
Experience comme prédicteur continu au niveau du médecin, et un
intercept aléatoire par DID, l’ID du médecin.
Etape 1
– Réalisation de modèles de régression logistique univarié pour toutes
les variables explicatives (= les facteurs de risque à tester), un
modèle pour chaque variable. On retient les variables pour l’étape
suivante (modèle multivarié) si la p-value est <=0.20
Ecriture du modèle univarié : fonction glmer()
ModGlmer<-glmer(VaràExpliquer ~ VarExplicative + (1|Vareffet aléatoire), data = BDD, family = binomial)
Exemple avec la variable CancerStage (variable qualitatitive)
# Vérification des effectifs
Code
# Vérification des effectifstable(BDD$remission,BDD$CancerStage)
Ecriture du modèle : même interprétation que le modèle sans effet aléatoire
Code
ModGlmer<-glmer(remission ~ CancerStage + (1|DID), data = BDD, family = binomial)summary(ModGlmer)Anova(ModGlmer, type =3)# Récupération des Odds Ratios (OR) des variables à titre d’exempletidy(ModGlmer, conf.int=TRUE, exponentiate=TRUE, effects="fixed")
Etape 2 : verification de la linéarité
Code
########################### Test de Box-Tidwell ############################### Modèle de régression logistiqueModGlmer<-glmer(remission ~ LengthofStay + (1|DID), data = BDD, family = binomial)summary(ModGlmer)#### Récupération des logit estimésprobabilities <-predict(ModGlmer, type ="response")LogOdds<-log(probabilities/(1-probabilities))#### Test de Box-Tidwelllibrary(car)boxTidwell(LogOdds~LengthofStay, data=BDD)# Interprétations#IL6: la p-value<0.05, le lien n'est pas linéaire. On transformera la variable en classe#LengthofStay: la p-value>0.05, le lien est linéaire.# pb d'itération avec PRC
Tableau résumé des pvalues obtenues à l’analyse univariée pour chaque variable
ModGlmerMulti<-glmer(remission ~ CancerStage + LengthofStay + (1|DID), data = BDD, family = binomial)library(performance)check_collinearity(ModGlmerMulti)
Les variables peuvent être laissées ensemble dans le même modèle
Etape 4 : Mise en oeuvre de l’analyse Multivariée
Code
ModGlmerMulti<-glmer(remission ~ CancerStage + LengthofStay + (1|DID), data = BDD, family = binomial)summary(ModGlmerMulti)Anova(ModGlmerMulti, type=3)
Etape 5 - Validité du modèle avec la vérification des résidus
Le package DHARMa offre une méthode générale pour vérifier si les
résidus d’un GLMM sont distribués en accord avec le modèle spécifié et
s’il y a absence de tendance résiduelle. Le package fonctionne en
simulant des réplicats de chaque observation en fonction du modèle
ajusté, puis en déterminant un “résidu standardisé” (= position relative
de la valeur observée par rapport aux valeurs simulées). Ex.: 0 si
l’observation est plus petite que toutes les simulations, 0.5 si elle
est au milieu, etc. Si le modèle représente bien les données, chaque
valeur du résidu standardisé entre 0 et 1 devrait être également
probable, donc les résidus standardisés devraient produire une
distribution uniforme entre 0 et 1.
## Sur le graphique 1 : graphique quantile-quantile des résidus standardisés
#### - Test de Kolmogorov-Smirnonov (KS) pour vérifier si la distribution des résidus s’éloigne de la distribution théorique
#### - Test de la dispersion : présence de sur ou sous-dispersion
#### - Test des valeurs extrêmes : y a t-il un excès de valeurs extrêmes dans les résidus.
## Sur le graphique 2 :
#### Résidus standardisés en fonction des rangs des valeurs prédites.
#### Courbe = régressions quantiles non-paramétriques pour Q1, Q2 et Q3.
#### Les 3 courbes doivent être horizontales (i.e. aucune
tendance résiduelle des résidus en fonction des prédictions). Si ce
n’est pas le cas, peut-être manque d’un effet important dans le modèle.
Code
ModGlmerMulti<-glmer(remission ~ CancerStage + LengthofStay + (1|DID), data = BDD, family = binomial)library(DHARMa)resid_sim <-simulateResiduals(ModGlmerMulti)plot(resid_sim)
Etape
6 – Présentations et interprétations des résultats du modèle final :
Tableau, graphe et comparaison poshoc (idem régression logistique sans
effet aléatoire)
Code
# Ajouter le tableau de résultats avec p-value globale de chaque variablelibrary(gtsummary)ModGlmerMulti %>%tbl_regression( exponentiate =TRUE) %>% add_global_p
Etape 9– Vérification de la distribution des effets aléatoires
Il est aussi important de vérifier la distribution des effets
aléatoires (hypothèse = ditribution approximativement normale). Un
graphe de type quantile-quantile permet de faire cette vérification
graphiquement.
Attention au nombre de groupe dans le jeu de données, vérifier la
normalité sur un petit nombre de points est souvent difficile
(vérification à titre indicatif).
Code
#Récupération des effets aléatoires pour chaque groupe estimés par le modèleRandom <-ranef(ModGlmerMulti)$DID# QQplotp <-ggplot(Random, aes(sample = Random$`(Intercept)`))p +stat_qq() +stat_qq_line()+theme_bw()