présentation Glmer pour Wiki

quarto lien https://www.youtube.com/watch?v=Cwg7tdSdRvY

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ée
BDD<-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.

Code
BDD_Reg<-BDD[,-1]
  BDD_Reg|> tbl_summary(by = remission)
Characteristic 0, N = 6,0041 1, N = 2,5211
IL6 3.41 (1.97, 5.48) 3.18 (1.83, 5.24)
CRP 4.40 (2.73, 6.68) 4.23 (2.59, 6.47)
CancerStage
I 1,580 (26%) 978 (39%)
II 2,357 (39%) 1,052 (42%)
III 1,301 (22%) 404 (16%)
IV 766 (13%) 87 (3.5%)
LengthofStay 6 (5, 6) 5 (5, 6)
Experience 17 (14, 20) 19 (16, 21)
1 Median (IQR); n (%)

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

  1. 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 effectifs
table(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’exemple
tidy(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 logistique
ModGlmer<-glmer(remission ~  LengthofStay + (1|DID), data = BDD, family = binomial)
summary(ModGlmer)

#### Récupération des logit estimés
probabilities <- predict(ModGlmer, type = "response")
LogOdds<-log(probabilities/(1-probabilities))
#### Test de Box-Tidwell
library(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

Code
#CancerStage : => p-value < 2.2e-16
#LengthofStay : => p-value < 2e-16

Etape 3 : calcul des VIFS

Code
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 variable
library(gtsummary)
ModGlmerMulti %>% 
tbl_regression( exponentiate = TRUE) %>%
  add_global_p
Characteristic OR1 95% CI1 p-value
CancerStage <0.001
I
II 0.66 0.57, 0.76
III 0.37 0.30, 0.44
IV 0.10 0.07, 0.13
LengthofStay 0.89 0.83, 0.95 <0.001
1 OR = Odds Ratio, CI = Confidence Interval
Code
# Représentation graphique
library(finalfit)
dependent='remission'
explanatory=c("CancerStage" , "LengthofStay")      
BDD |>
  or_plot(dependent, explanatory, table_text_size=3, title_text_size=14,
plot_opts=list(xlab("OR, 95% CI"), theme(axis.title = element_text(size=))))

Etape 7 - Visualisation graphique des effets du modèle
Code
library(ggeffects)
cowplot::plot_grid(plotlist = plot(ggeffect(ModGlmerMulti)))

Etape 8 – Calcul d’un pseudo R² pour les modèles glmer

# R2m= R²marginal = variance expliquée par les effets fixes

# R2C = R² conditionnel = variance expliquée par les effets fixes et aléatoires

Code
library(MuMIn)                                                 
r.squaredGLMM(ModGlmerMulti)
                   R2m       R2c
theoretical 0.06750995 0.5855338
delta       0.05680732 0.4927067

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èle
Random <- ranef(ModGlmerMulti)$DID

# QQplot
p <- ggplot(Random, aes(sample = Random$`(Intercept)`))
p + stat_qq() + stat_qq_line()+theme_bw()