Introduction aux GLMM avec données de proportion

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) !

Contexte

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 :

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.

Caractéristiques des données et approche statistique

Données bornée

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.

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.

A propose des GLMM

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.

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 :

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

Ce qu'il faut retenir

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

Trois éléments me semblent importants :

  • C’est que le GLMM va nous permette de comparer les niveaux moyens des pourcentages observés pour les traitements A et B.
  • 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.

  • 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})}} \]

Les données

Importation et structure

library(tidyverse)
library(here)
library(knitr)

mydata <- read.csv2(here::here("data", "Cellules.csv"))
head(mydata)

    ##   vivant mort boite  trt
    ## 1    138   24    B1 TrtA
    ## 2     94   14    B1 TrtA
    ## 3     46    8    B1 TrtA
    ## 4     38    3    B1 TrtA
    ## 5     78    3    B1 TrtA
    ## 6     91    2    B1 TrtA 

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.

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`.

mydata <- mydata %>% 
        mutate(freq_viv =(vivant/(vivant+mort)))

head(mydata)

    ##   vivant mort boite  trt  freq_viv
    ## 1    138   24    B1 TrtA 0.8518519
    ## 2     94   14    B1 TrtA 0.8703704
    ## 3     46    8    B1 TrtA 0.8518519
    ## 4     38    3    B1 TrtA 0.9268293
    ## 5     78    3    B1 TrtA 0.9629630
    ## 6     91    2    B1 TrtA 0.9784946 

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

mydata$freq_viv <- round(mydata$freq_viv,3)
head(mydata)

    ##   vivant mort boite  trt freq_viv
    ## 1    138   24    B1 TrtA    0.852
    ## 2     94   14    B1 TrtA    0.870
    ## 3     46    8    B1 TrtA    0.852
    ## 4     38    3    B1 TrtA    0.927
    ## 5     78    3    B1 TrtA    0.963
    ## 6     91    2    B1 TrtA    0.978 

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).

ggplot(mydata,aes(x=boite,y=freq_viv, fill=trt, colour=trt) )+
        geom_jitter(width=0.25) +
        geom_boxplot(outlier.alpha = 0, alpha=0.5) 

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 :

mydata %>% 
        group_by(trt, boite) %>% 
        summarise(nb=n(),
                  avg=mean(freq_viv))

    ## # A tibble: 8 x 4
    ## # Groups:   trt [2]
    ##   trt   boite    nb   avg
    ##      
    ## 1 TrtA  B1        8 0.922
    ## 2 TrtA  B2       15 0.906
    ## 3 TrtA  B3       18 0.944
    ## 4 TrtA  B4        6 0.876
    ## 5 TrtB  B5       19 0.836
    ## 6 TrtB  B6       24 0.808
    ## 7 TrtB  B7       15 0.735
    ## 8 TrtB  B8       12 0.744 

Ajustement du GLMM

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 :

mydata$y <- cbind(mydata$vivant, mydata$mort)
head(mydata)

    ##   vivant mort boite  trt freq_viv y.1 y.2
    ## 1    138   24    B1 TrtA    0.852 138  24
    ## 2     94   14    B1 TrtA    0.870  94  14
    ## 3     46    8    B1 TrtA    0.852  46   8
    ## 4     38    3    B1 TrtA    0.927  38   3
    ## 5     78    3    B1 TrtA    0.963  78   3
    ## 6     91    2    B1 TrtA    0.978  91   2

head(mydata$y)

    ##      [,1] [,2]
    ## [1,]  138   24
    ## [2,]   94   14
    ## [3,]   46    8
    ## [4,]   38    3
    ## [5,]   78    3
    ## [6,]   91    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`.
library(lme4)
glmer1 <- glmer(y~trt+(1|boite), family=binomial, data=mydata)
summary(glmer1)
 
    ## Generalized linear mixed model fit by maximum likelihood (Laplace
    ##   Approximation) [glmerMod]
    ##  Family: binomial  ( logit )
    ## Formula: y ~ trt + (1 | boite)
    ##    Data: mydata
    ## 
    ##      AIC      BIC   logLik deviance df.resid 
    ##    525.4    533.6   -259.7    519.4      114 
    ## 
    ## Scaled residuals: 
    ##     Min      1Q  Median      3Q     Max 
    ## -4.8299 -0.5462  0.1386  0.8382  2.3740 
    ## 
    ## Random effects:
    ##  Groups Name        Variance Std.Dev.
    ##  boite  (Intercept) 0.01527  0.1236  
    ## Number of obs: 117, groups:  boite, 8
    ## 
    ## Fixed effects:
    ##             Estimate Std. Error z value Pr(>|z|)    
    ## (Intercept)   2.2376     0.1090  20.537  < 2e-16 ***
    ## trtTrtB      -0.8035     0.1349  -5.957 2.58e-09 ***
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ## 
    ## Correlation of Fixed Effects:
    ##         (Intr)
    ## trtTrtB -0.808 

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.

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)`

surdisp <- 519.4/114
 
surdisp
## [1] 4.55614 

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 :

overdisp_fun <- function(model) {
        rdf <- df.residual(model)
        rp <- residuals(model,type="pearson")
        Pearson.chisq <- sum(rp^2)
        prat <- Pearson.chisq/rdf
        pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
        c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}

overdisp_fun(glmer1)
##        chisq        ratio          rdf            p 
    ## 1.941786e+02 1.703321e+00 1.140000e+02 4.216640e-06 

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.

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 :

glmer2 <- glmer(y~trt+(1|rep), family=quasibinomial, data=mydata)

Error in lme4::glFormula(formula = y ~ trt + (1 | rep), data = mydata, : 
"quasi" families cannot be used in glmer 

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`:

library(MASS)

glmmPQL1 <- glmmPQL(y~trt, random=~1|boite, family=quasibinomial, data=mydata)

summary(glmmPQL1)

    ## Linear mixed-effects model fit by maximum likelihood
    ##  Data: mydata 
    ##   AIC BIC logLik
    ##    NA  NA     NA
    ## 
    ## Random effects:
    ##  Formula: ~1 | boite
    ##         (Intercept) Residual
    ## StdDev:  0.08275886 1.320594
    ## 
    ## Variance function:
    ##  Structure: fixed weights
    ##  Formula: ~invwt 
    ## Fixed effects: y ~ trt 
    ##                  Value Std.Error  DF  t-value p-value
    ## (Intercept)  2.2338383 0.1243801 109 17.95977  0.0000
    ## trtTrtB     -0.7982217 0.1467616   6 -5.43890  0.0016
    ##  Correlation: 
    ##         (Intr)
    ## trtTrtB -0.847
    ## 
    ## Standardized Within-Group Residuals:
    ##        Min         Q1        Med         Q3        Max 
    ## -3.6902089 -0.4189964  0.1369458  0.6488932  1.8356178 
    ## 
    ## Number of Observations: 117
    ## Number of Groups: 8 

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é:

cc <- coef(summary(glmer1))
phi <- overdisp_fun(glmer1)["ratio"]
cc <- within(as.data.frame(cc),
    {   `Std. Error` <- `Std. Error`*sqrt(phi)
        `z value` <- Estimate/`Std. Error` `Pr(>|z|)` <- 2*pnorm(abs(`z value`), lower.tail=FALSE) }) printCoefmat(cc,digits=3) ## Estimate Std. Error z value Pr(>|z|)    
    ## (Intercept)    2.238      0.142   15.74   <2e-16 ***
    ## trtTrtB       -0.804      0.176   -4.56    5e-06 ***
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

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 )

Résultats

summary(glmmPQL1)

    ## Linear mixed-effects model fit by maximum likelihood
    ##  Data: mydata 
    ##   AIC BIC logLik
    ##    NA  NA     NA
    ## 
    ## Random effects:
    ##  Formula: ~1 | boite
    ##         (Intercept) Residual
    ## StdDev:  0.08275886 1.320594
    ## 
    ## Variance function:
    ##  Structure: fixed weights
    ##  Formula: ~invwt 
    ## Fixed effects: y ~ trt 
    ##                  Value Std.Error  DF  t-value p-value
    ## (Intercept)  2.2338383 0.1243801 109 17.95977  0.0000
    ## trtTrtB     -0.7982217 0.1467616   6 -5.43890  0.0016
    ##  Correlation: 
    ##         (Intr)
    ## trtTrtB -0.847
    ## 
    ## Standardized Within-Group Residuals:
    ##        Min         Q1        Med         Q3        Max 
    ## -3.6902089 -0.4189964  0.1369458  0.6488932  1.8356178 
    ## 
    ## Number of Observations: 117
    ## Number of Groups: 8 

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

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 :

p_trtA = exp(2.23)/(1+exp(2.23)) 
p_trtA
## [1] 0.9029114

p_trtB=exp(2.23-0.80)/(1+exp(2.23-0.80)) 
p_trtB
## [1] 0.8069013 

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 :

glmmPQL2 <- glmmPQL(y~trt-1, random=~1|boite, family=quasibinomial, data=mydata)
summary(glmmPQL2)

    ## Linear mixed-effects model fit by maximum likelihood
    ##  Data: mydata 
    ##   AIC BIC logLik
    ##    NA  NA     NA
    ## 
    ## Random effects:
    ##  Formula: ~1 | boite
    ##         (Intercept) Residual
    ## StdDev:  0.08275886 1.320594
    ## 
    ## Variance function:
    ##  Structure: fixed weights
    ##  Formula: ~invwt 
    ## Fixed effects: y ~ trt - 1 
    ##            Value  Std.Error DF  t-value p-value
    ## trtTrtA 2.233838 0.12438015  6 17.95977       0
    ## trtTrtB 1.435617 0.07790088  6 18.42876       0
    ##  Correlation: 
    ##         trtTrA
    ## trtTrtB 0     
    ## 
    ## Standardized Within-Group Residuals:
    ##        Min         Q1        Med         Q3        Max 
    ## -3.6902089 -0.4189964  0.1369458  0.6488932  1.8356178 
    ## 
    ## Number of Observations: 117
    ## Number of Groups: 8

# quantile 97.5 selon distribution de student
T <-qt(0.975,6)
T
## [1] 2.446912

binfA = 2.23-2.45*0.12 
bsupA = 2.23 +2.45*0.12
IC_infA = exp(binfA)/(1+exp(binfA))
IC_supA = exp(bsupA)/(1+exp(bsupA))

IC_infA
## [1] 0.873912

IC_supA
## [1] 0.9258073

binfB = 1.44-2.45*0.08 
bsupB = 1.44 +2.45*0.08
IC_infB = exp(binfB)/(1+exp(binfB))
IC_supB = exp(bsupB)/(1+exp(bsupB))

IC_infB
## [1] 0.7762595
IC_supB
## [1] 0.8369899 

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]

Visualisation des résultats

colA<- "#FF7F00"
colB <- "#1E90FF"
mycol=c(colA, colB)

c("#FF7F00", "#1E90FF", "#FFFFFF")
## [1] "#FF7F00" "#1E90FF" "#FFFFFF"

ggplot(mydata,aes(x=boite,y=freq_viv, fill=trt, colour=trt) )+
        geom_jitter(width=0.25)+
        geom_boxplot(outlier.alpha = 0, alpha=0.5)+
        scale_y_continuous(limits=c(0,1.5))+
        scale_colour_manual(values=mycol)+
        scale_fill_manual(values=mycol)+
        geom_segment(x=0.75, xend=4.25, y=1.1, yend=1.1,col=colA)+
        # geom_segment(x=0.75, xend=0.75, y=1.05, yend=1.1,col=colA)+
        # geom_segment(x=4.25, xend=4.25, y=1.05, yend=1.1,col=colA)+

        geom_segment(x=4.75, xend=8.25, y=1.1, yend=1.1,col=colB)+
        # geom_segment(x=4.75, xend=4.75, y=1.05, yend=1.1,col=colB)+
        # geom_segment(x=8.25, xend=8.25, y=1.05, yend=1.1,col=colB)+
        
        geom_segment(x=2.25, xend=6.25, y=1.25, yend=1.25,col="black")+
        geom_segment(x=2.25, xend=2.25, y=1.20, yend=1.25,col="black")+
        geom_segment(x=6.25, xend=6.25, y=1.20, yend=1.25,col="black")+
        annotate("text", x = 4, y = 1.3, label = c("p < 0.001")) 

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.

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édit photo : qimono

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

20 réponses

  1. 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

  2. Bonjour,
    Merci pour cet article
    Lorsque je réalise cette analyse en l’adaptant à mes données j’ai le message d’erreur suivant
    “fixed-effect model matrix is rank deficient so dropping columns / coefficients”
    Je ne le comprend pas Pourriez vous m’expliquer ?

  3. Bonjour et un grand merci pour votre travail et la qualité de vos articles.
    Je souhaite bien comparer des proportions dans le cadre de mes recherches de thèse mais je ne suis pas certains que le GLMM soit bien adapté et je souhaiterai obtenir votre avis.

    L’expérience se présente ainsi :
    Je travaille sur des cultures de microalgues toxiques soumises à un traitement chimique.
    Ces algues produisent 3 toxines différentes.
    On prélève différents échantillons à différents temps après l’application du traitement.
    J’ai donc:
    – un facteur traitement
    – un facteur temps sur des échantillons non indépendants car on prélève plusieurs fois dans la même culture au cours du temps.
    – enfin chaque culture est réalisée en triplicat, ce qui n’est pas l’idéal.

    Je souhaiterai déterminer si le traitement et potentiellement le temps on un effet sur le profil toxinique, autrement dit j’ai exprimé ma quantité de toxine par rapport à la somme des quantités des 3 toxines par cellules.
    A un temps T, pour un traitement 1, j’ai donc des quantités de toxines exprimées comme-ci :
    – Toxine A = 10% du total de toxine par cellule
    – Toxine B = 80% du total de toxine par cellule
    – Toxine C = 10% du total de toxine par cellule

    Dans cette configuration est-ce bien le GLMM qui est adapté pour répondre à ma question ou suis-je à coté de la plaque ?

    Merci encore.

    1. Bonsoir,
      à priori je partirai plutôt sur un modèle linéaire à effet mixte car votre variable réponse est un ratio de quantité qui est elle même une variable numérique continue non bornée. Les GLMM pour données de proportion c’est pour des ratio de données de comptages. Vous trouverez des exemples dans le R book (en cherchant un peu vous le trouverez en pdf 😉
      Bonne continuation

  4. Bonjour, Claire.
    Merci pour cet article très intéressant et très clairement présenté.
    Je voudrais savoir si la définition de l’odds ratio (OR) pour les Glm est valable et applicable pour les Glmm.
    Par avance, merci.
    Cordialement,

  5. Hello Claire 🙂

    J’avais exactement la même problématique que tu présentes au début de cet article, avec un facteur aléatoire de groupe, et un effet fixe dont je voulais calculer des intervalles de confiance, et ton article m’a énormément aidé !

    Les modèles mixtes ça remonte à très loin, et je trouve que la plupart des articles sur les internets ressemblent à mes anciens cours de stats… pas très clairs 😀

    Mais ton article est top et m’a beaucoup aidé !

    J’ai encore quelques questions et je me demande si tu as des livres à recommander pour creuser la question (je vais déjà regarder les 2 articles que tu as proposés) et comprendre un peu mieux les maths et le R utilisés. Par exemple :
    – Pourquoi ne peut-on pas donner les données brutes 0/1 en entrée ? Je remarque que ça donne des résultats proches mais pas identiques. Notamment au niveau de la détection de la surdispersion.
    – Pourquoi et comment la surdispersion ? Je ne comprends pas d’où ça vient. Les données sont forcément issues d’une Bernoulli, et donc d’une binomiale quand on les ajoute, du coup.. pourquoi ?

    Merci 😉

    1. Hello Charles,
      je suis ravie que l’article t’ait aidé.
      Concernant les livres, ce que j’ai lu de plus compréhensible et de plus utile sur les LMM ou GLMM est dans le livre d’Alain Zuur : https://www.amazon.fr/Effects-Models-Extensions-Ecology-Author/dp/B007S7GBTM/ref=sr_1_1?__mk_fr_FR=ÅMÅŽÕÑ&keywords=alain zuur&qid=1578776225&sr=8-1
      Tu verras aux avis sur Amazon, que je ne suis pas la seule à le trouver chouette !

      En ce qui concerne l’utilisation des données brutes, ça veut dire que tu crées par exemple 138 lignes avec une réponse codée 1 et 24 lignes avec une réponse codée 0 (première ligne de mon exemple). Ensuite la fonction de lien est sensée être identique….je n’ai pas de réponse. Est ce que les petites différences sont sur les coeff ou sur les erreurs standard, ou les deux ?
      Pour la surdispersion, ma compréhension est que le modèle calcule les erreurs standard des paramètres en faisant l’hypothèse que les données sont distribuées selon une distribution binomiale. Quand il y a surdispersion c’est que la dispersion réelle (cad observée) est supérieure à ce qui est attendue en théorie avec la distribution binomiale. Il faut donc corriger pour que l’erreur standard des paramètres soit calculée correctement, et donc en bout de chaîne que les tests soient exactes.
      J’espère que ça t’aide 😉

  6. Bonjour Claire,
    j’apprécie bien les chapitres traités dans ton blog: félicitations!
    Cependant, je souhaite apprendre à analyser dans R en utilisant lmer.
    Pouvez-vous m’aider?

    Merci d’avance

  7. Bonjour !
    Merci beaucoup pour ce cours, il est très bien expliqué.

    Par ailleurs, il me semble que les résidus doivent se comporter normalement avec un GLM comme un GLMM, je voulais savoir comment procéder quand ce n’est pas vraiment le cas (transformation des données ? autre approche à utiliser ?).

    Dans un de mes essais, je fonctionne par bloc de répétitions de 2 traitements (comme ton exemple) mais j’effectue mon comptage (vivant/mort) sur l’ensemble du bloc (et pas sur plusieurs zones comme toi avec tes photos). Je me demandais si cela changeais quelque chose au modèle ou si cela restait un GLMM ?

    Merci d’avance pour ta réponse 🙂

  8. Bonjour Claire,
    merci beaucoup pour cet article très bien écrit et très utile.

    Je voulais vous demander si vous connaissiez par hasard une toolbox permettant de faire du Generalized NonLinear Mixed Model (GNLMM) ? J’ai une variable y binaire (0 1) et des prédicteurs organisé selon un modèle non linéaire.
    Je sais que la librairie BRMS dans R est capable de le faire mais c’est de l’estimation Bayesienne (je n’ai rien contre le bayesien, mais je comprends mieux pour l’instant l’estimation par Maximum de Vraisemblance classique).

    Merci d’avance !
    Vincent

Laisser un commentaire

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