Introduction aux GLMM avec données de proportion

GLMM logiciel R pourcentages

 

Pour faire suite au tutoriel sur les GLM avec données de comptage, et pour répondre aux demandes de certains d’entre vous, je vous propose ici une introduction aux GLMM avec données de proportion, sous la forme d’un petit tutoriel. Les GLMM (pour Generalized Linear Mixed Models) sont des modèles linéaires généralisés à effets mixtes. Ils sont employés pour analyser des réponses qui ne sont pas numériques continues non bornées – c’est le cas des données de comptages, des réponses binaires (homme/femme par exemple) ou encore, comme ici des proportions, (ça c’est pour la partie GLM), et lorsque les données ne sont pas indépendantes (ça c’est pour la partie mixte) !

1. Contexte
1.1 Le plan expérimental

Les données analysées dans ce tutoriel sont issues d’une expérimentation qui vise à évaluer les effets de deux traitements sur la survie ou la mort de cellules. Pour chaque traitement (traitement A ou traitement B), 4 boites de cultures sont utilisées (4 boites avec des cellules ayant reçu le traitement A et 4 boites avec des cellules ayant reçu le traitement B).

Pour chacune de ces boites, un certain nombre de photos sont prises (entre 6 et 24) sur une petite zone de la boite, et pour chaque photo, le nombre de cellules mortes et vivantes sont observées.

Le plan expérimental peut être résumé par ce schéma :

GLMM avec pourcentage

1.2 La question posée

La question qui est posée est “Est ce que les proportions (ou pourcentages) de cellules vivantes sont différentes selon le traitement reçu ?” Autrement formulé, “est ce qu’un des traitements est associé à une proportion plus importante de cellules vivantes ?”

D’un point de vu pratique, on cherche donc à faire une comparaison de deux pourcentages.

 

2. Caractéristiques des données et approche statistique

2.1 Données bornées

Compte tenu de la question posée, on va donc s’intéresser à la proportion de cellules vivantes. Cette proportion est une variable numérique continue, mais qui est un peu particulière, car elle est bornée entre 0 et 1. C’est cette spécificité qui rend l’emploi de modèles linéaires classiques caduque, et nous oriente vers l’utilisation d’un modèle linéaire généralisé (GLM pour Generalized Linear Model).

Il existe grossièrement trois types de GLM:

  • les GLM avec lien log et distribution d’erreur de poisson lorsque les réponses observées sont de type comptage,
  • les GLM avec lien logit et distribution d’erreur binomiale lorsque les réponses observées sont binaires (oui/non, mâle/femelle),
  • les GLM avec lien logit et distribution d’erreur binomiale lorsque les réponses observées sont des proportions. C’est ce dernier modèle qui nous intéresse.

2.2 Données non indépendantes

Néanmoins, les données possèdent une autre caractéristique : elles ne sont pas indépendantes. En effet, le plan expérimental nous laisse penser que, pour un traitement donné, les pourcentages de cellules vivantes observées dans une même boite sont plus liés entre eux (corrélés) que les pourcentages provenant de deux boites différentes.

Pour comprendre ce concept, on peut imaginer, par exemple, qu’une des boites à un défaut de fabrication qui entraîne une sur-mortalité des cellules. Dans ce cas, toutes les observations réalisées sur cette boite sont impactées. Et donc les pourcentages de cellules vivantes observées
sur les images de cette boite sont plus semblables entre eux que ceux observées sur les images de deux boites différentes.

C’est cette seconde caractéristique, la non-indépendance, qui nous dirige finalement vers l’utilisation d’un GLMM (Generalized Linear Mixed Model ou modèle linéaire généralisé mixtes). Pour faire simple Un GLMM (Generalized Linear Mixed Model ou modèle linéaire généralisé mixtes) est un GLM avec une fonctionnalité supplémentaire qui lui permet de prendre en considération la non indépendance des données. Ici, on va donc utiliser un GLMM avec un lien logit et une distribution binomiale des erreurs.

Remarque très importante :

Utiliser un GLMM, et non pas un GLM, en cas de non indépendance des données ne relève pas du détail ! C’est CAPITAL ! Sans cette prise en compte de la corrélation entre certaines données, la déviance résiduelle du modèle employé pour analyser les données sera biaisée. Cela conduira à une estimation, elle aussi biaisée, de l’erreur standard des paramètres, et en bout de chaîne à une p value erronée.

 

3. A propos des GLMM
3.1 Pourquoi “mixte” ?

Un GLMM est dit “mixte”, car il comporte au moins un effet dit “fixe” (la variable dont on souhaite évaluer l’effet, ici le traitement) et au moins un effet dit “aléatoire” (la variable de regroupement, ici la boite). Les effets aléatoires ne sont pas évalués, ils servent seulement à indiquer au modèle que les données ne sont pas indépendantes pour une boite donnée. C’est ce qui permet à la déviance résiduelle d’être bien estimée, et ainsi à l’erreur standard des paramètres de ne pas être biaisée, et aux final d’obtenir des résultats fiables. En pratique, cela passe par l‘addition d’une quantité dans la matrice de variance- covariance des résidus.

3.2 Complexité

Les GLMM sont des approches plutôt complexes. Pour vous en convaincre, vous pouvez consulter cette FAQ;  rédigée, entre autres, par Ben Bolker, un des grand spécialiste des GLMM.

Voici ce qu’il écrit en introduction :

GLMM avec R

Je ne vais donc pas ici rentrer dans les détails de cette approche, mais uniquement me concentrer sur les aspects pratiques.

3.3 Ce qu’il faut retenir :

Trois éléments me semblent importants :

  1. c’est que le GLMM va nous permette de comparer les niveaux moyens des pourcentages observés pour les traitements A et B.

  2. La fonction de lien utilisée est la fonction logit, définie par :
    $$ logit(p) = ln(\frac{p}{q}) = \sum_{j=1}^{p} \beta_j\;X_{ij} $$

Pour plus d’informations sur les fonctions de lien, je vous renvoie vers le tutoriel sur les GLM avec données de comptage.

  1. Ceci implique que pour estimer les proportions de cellules vivantes et mortes à partir des paramètres estimés par le modèle (le GLMM), il sera nécessaire d’employer la fonction inverse :

$$ p = \frac{ \exp^{(\sum_{j=1}^{p} \beta_j\;X_{ij})}}{1+exp^{(\sum_{j=1}^{p} \beta_j\;X_{ij})}} $$

 

4. Les données

4.1 Importation et structure

 

 

Le jeu de données est au format tidy avec :

  • une colonne “vivant” : qui contient le nombre de cellules vivantes observées sur une image
  • une colonne “mort” : qui contient le nombre de cellules mortes observées sur une image
  • une colonne “boite” qui renseigne sur la boite dans laquelle l’observation à eu lieu
  • une colonne “trt” qui renseigne sur le traitement reçu par les cellules.

4.2 Création d’une variable de proportion (freq_viv)

Les données ne contiennent pas de variable correspondant à la proportion de cellules vivantes. Cette variable sera utile pour représenter visuellement les données, nous allons donc la créer à l’aide de la fonction mutate() du package dplyr qui appartient au super package tidyverse.

 

On peut alors arrondir ces proportions, en utilisant la fonction round() comme ceci :

 

5. Exploration des données

Pour avoir une première idée des données, le plus parlant est de réaliser une représentation graphique. Je choisis d’utiliser des boxplots en superposant les données. Pour cela, j’utilise le package ggplot2 et ses couches geom_jitter() et geom_boxplot(). L’argument width=0.25 permet aux points d’être contenus dans les boites. L’argument alpha=0.5 permet d’avoir une couleur dans les boites plus claire que celle des points, et l’argument outlier.alpha=0 permet aux outliers de ne pas être représentés deux fois (une fois par la couche geom_jitter et une fois par la couche geom_boxplot).

 

GLMM boxplot

A première vue, le pourcentage de cellules vivantes est plus important en présence du traitement A qu’en présence du traitement B. L’analyse statistique, et utilisation du GLMM, nous permettra d’évaluer si cette différence est significative ou non.

On peut également résumer les données sous la forme d‘une table, en indiquant le nombre d’images prises par boite et le pourcentage moyen de cellules vivantes :

 

6. Ajustement du GLMM
6.1 création d’une variable combinée

Pour ajuster le GLMM, il est nécessaire de créer une variable combinée (nombre de cellules vivantes, nombre de cellules mortes). C’est cette nouvelle variable qui sera considérée comme la réponse dans le modèle :

 

6.2 Ajustement du modèle

Pour cela, on utilise la fonction glmer() du package lme4, en spécifiant :

  • le traitement (variable trt) comme effet fixe,
  • la boite comme effet aléatoire, avec la syntaxe (1|boite),
  • la fonction de lien logit et la distribution binomial des erreurs avec family=binomial.

 

 

La fonction summary() renvoit plusieurs éléments. Les plus importants sont:

  • des éléments de qualité du modèle (AIC, BIC, deviance, etc..),
  • des éléments concernant les effets aléatoires,
  • des éléments concernant les effets fixes,,
  • la première ligne concerne le traitement A,
  • la seconde ligne le traitement B, mais exprimé relativement au traitement A,
  • la première colonne correspond aux estimations des paramètres,
  • la seconde aux estimations des erreurs standard des paramètres,
  • la troisième correspond à la statistique du test d’égalité à 0 des paramètres,
  • la quatrième correspond à la p value du test.

Remarques:

  • Les éléments des effets fixes sont exprimés dans l’échelle logit. Pour les estimer dans l’échelle des proportions, il sera nécessaire d’utiliser la fonction inverse du logit, décrite précédemment.
  • Le fait que le coefficient relatif au traitement B soit <0 signifie que le niveau moyen des proportions de cellules vivantes est plus faible en présence du traitement B, par rapport au
    traitement A. Pour savoir si nous pouvons avoir confiance dans la valeur de la p value, nous devons d’abord évaluer s’il existe une surdispersion.

6.3 Evaluation de la surdispersion

Il y a surdispersion (overdispersion en anglais) lorsque la variance des réponses observée est supérieure à la variance théorique, ici définie par la loi binomiale. Dans cette situation, l’erreur standard des paramètres est sous-estimée. Ceci peut conduire à une p value excessivement faible, et donc aboutir à une conclusion erronée. Il est donc important d’évaluer la présence d’une éventuelle surdispersion.

De manière grossière, la surdispersion peut s’évaluer par le rapport de la déviance sur le nombre de degrés de liberté, qui figurent sur la sortie de la fonction summary(glmer1)

 

Une surdispersion est mise en évidence lorsque le ratio est supérieur à 1 (autour de 1.5 par exemple). Le ratio étant ici de l’ordre de 5, on conclue à la présence d’une surdispersion.

Une autre approche, sous la forme d’une fonction, est proposée par Ben Bolker, ici :

 

La p-value du test est < 0.05, on conclut également à la présence d’une surdispersion.

En cas de surdispersion, il est nécessaire d’utiliser une autre distribution des erreurs, telle que la structure “quasibinomial”. Celle-ci va conduire à une augmentation des erreurs standard des paramètres du modèle.

6.4 Utilisation de la loi quasi binomiale

Il n’est malheureusement pas possible d’utiliser la distribution “quasiobinomial” au sein de la fonction glmer(), puisque cela génère un message erreur :

Dans cette situation, deux solutions sont envisageables. La première consiste à employer la fonction glmmPQL() du package MASS. Dans ce cas, l’effet aléatoire doit être défini avec cette syntaxe : random=~1|boite:

On peut, voir que les erreurs standard des paramètres ont augmenté (par rapport à celles du modèle glmer1).

Ben Bolker propose également une autre méthode qui consiste à ajuster la table des coefficients, c’est-à-dire en multipliant l’erreur type par la racine carrée du facteur de dispersion et en recalculant les valeurs Z et p en conséquence. Le code suivant est donné:

 

Les erreurs standard sont alors encore plus élevées qu’avec le modèle glmmPQL. A présent que la surdispersion a été prise en considération, on peut analyser les résultats (ici par exemple ceux du modèle glmmPQL1 )

7. Résultats

 

Comme expliqué précédemment, le fait que le coefficient relatif au traitement B soit < 0 signifie que le niveau moyen des proportions de cellules vivantes est plus faible en présence du traitement B, par rapport au traitement A.

La p value relative au traitement B correspond à la comparaison des niveaux moyens des proportions de cellules vivantes dans les deux traitements. Sa valeur < 0.05 nous indique que le niveau moyen des proportions de cellules vivantes en présence de traitement B est significativement plus faible qu’en présence de traitement A.

7.1 Estimation des pourcentages de cellules vivantes

Comme expliqué précédemment, pour calculer ces pourcentages, il est nécessaire d’utiliser la fonction inverse du logit :

 

7.2 Estimation des intervalles de confiance des pourcentages

L’estimation des intervalles de confiance est plus difficile, car selon la paramétrisation du modèle, l’erreur standard de la seconde ligne est en réalité une erreur standard de différence.

Un moyen de ne pas se tromper est d’utiliser une autre paramétrisation du modèle (en ajoutant -1 à la fin des effets fixes). Dans cette paramétrisation,  la seconde ligne n’est plus exprimée relativement à la première ligne. L’erreur standard est donc bien celle correspondant au niveau moyen prédit de proportion de cellule vivante en présence de traitement B. On peut alors s’en servir pour construire l’intervalle de confiance à 95% de cette proportion :

 

Ainsi, les proportions de cellules vivantes pour les traitement A et B, avec leurs intervalles de confiance à 95%,  sont respectivement :

  • 0.90 [0.87 ; 0.93] et
  • 0.81 [0.78 ; 0.84]

8. Visualisation des résultats

 

GLMM résulat

9 Pour aller plus loin

Si le sujet des GLMM vous intéresse, je vous recommande ces deux publications :

 

Voilà, j’espère que cet article permettra de vous familiariser avec les GLMM, et de comprendre dans quelles situations les employer.

Et si cet article vous a plu, ou vous a été utile, n’oubliez pas de le partager 😉

 

Crédit photo : qimono

Crédit data : Rémi G 😉 (Merci !)

 

Poursuivez votre lecture :

Partager l'article
  •  
  •  
  •  
  •  
  •  
    27
    Partages
  • 27
  •  
  •  
  •  
  •  

6 commentaires

  1. EL GAOUBI Fodil Répondre

    Bonjour,

    C’est un document important et par conséquence, je vous remercie vivement pour ce don gratuit de pleine utilité et de dévouement.

    Avec mes encouragements pour Claire.

    Fodil EL GAOUBI

Laisser un commentaire

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