GLMM: application/ groupes non équilibrés

Je réalise une analyse stat pour laquelle j'ai 2 VI inter (groupe et sexe) ainsi que 1 VI intra (shift). Toutes mes VI ont 2 modalités et sont codées en -0.5, +0.5.
Les 2 VD que j'étudie sont: les réponses correctes et le temps de réponse.
(Mon nombre d'observations par participant varie entre 14 et 22)
Dans mon groupe patient, j'ai 53 hommes et 47 femmes.
Dans mon groupe contrôle, j'ai environ 60 hommes et 115 femmes.

J'ai choisi de réaliser 2 GLMM (mais je ne suis pas experte du tout et j'aurais des questions).
Pour les réponses correctes:

Code:

glm2A <- glmer (CR2 ~ GroupC * ShiftC * SexC + (1 + ShiftC ||ID) + (1| stim),
               data = ESTdfWO1,
               family = binomial (link = "logit"),
               control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))


Pour le temps de réponse (j'ai préféré ce type de modèle à un LMM avec transformation du RT en log, ce qui est souvent fait, suite à la lecture de cet article: https://www.frontiersin.org/articles/10.3389/fpsyg.2015.01171/full; mais toute remarque à ce sujet est bienvenue):

Code:

glmRT2 <- glmer (R2 ~ GroupC * ShiftC * SexC + (1 + ShiftC |ID) + (1| stim),
              data = ESTdfWO2,
              family = inverse.gaussian (link = "identity"))


L'un de mes encadrants m'indique qu'il n'est pas bien d'avoir des groupes déséquilibrés, et qu'il faudrait que j'enleve des particpants femmes contrôle. Je comprends que cela soit le cas pour les ANOVA (surtout dans mon cas où il y a une interaction); néanmoins, je ne suis pas certaine que cela pose problème avec les GLMM.

J'ai essayé de lire sur le sujet et je ne trouve pas d'information claire. Ainsi, j'aurais aimé savoir:

1. Les groupes inégaux sont ils problématiques pour les GLMM? Connaitriez vous un papier que je puisse lire qui traite le sujet?
   
2. Si la réponse est oui: quel type d'analyse stat devrais-je utiliser à la place? Quelle est la meilleure stratégie à adopter dans mon cas?
   
3. Si la réponse est non c'est ok: y a t'il des conditions d'applications à respecter pour les GLMM? (j'ai aussi du mal à trouver quelque chose de clair à ce sujet)
   
4. Pour investiguer l'interaction, en générale j'utilise le package emmeans et pairs (dans R). Savez-vous s'il est ok d'utiliser ce type de tests post-hoc dans le cas de groupes inégaux (notamment en spécifiant un ajustement "Dunnett" qui semble plus ok en cas de groupes inégaux...)?
   

Merci d'avance à ceux qui m'apporteront des réponses car j'ai beaucoup de mal à en trouver...

Rapidement,

- Non, les groupes inégaux ne sont pas à ce point problématiques. Il ne faut pas tout de même aller vers trop de déséquilibre, ce qui n'est vraiment le problème dans votre cas je pense.

- A mon avis, il est complètement inutile d'aller jusqu'à l'interaction de niveau 3. Ca allourdi le modèle inutilement. Ca augmente sensiblement le temps de calcul, et ça donne des résultats généralement non (ou très difficilement) interprétables. Je remplacerait "GroupC * ShiftC * SexC" pas "(GroupC + ShiftC + SexC)^2".

- Que représente la variable "CR2" ? Il serait intéressant de savoir comment vous codez ceci.

- Vous être dans le cas binomiale, or vous parlez de "les réponses correctes et le temps de réponse". Rien de binomial ici. Je soupçonne une erreur ici.

- Pourquoi coder -0.5, +0.5 vos variables binaires ? Quelle est l’intérêt ici ? Pourquoi ne pas les coder "male/female" et "groupe1/groupe2", etc? Ca serait largement plus clair, je pense.

HTH, Eric.

Merci beaucoup pour votre réponse.

• auriez vous pare hasard une piste de lecture permettant d'expliquer pourquoi les groupes inégaux ne sont pas problématiques?
  
  
• concernant l'interaction de niveau 3, c'est le modèle le plus complet (et j'ai pré-enregistré cette étude ainsi, avec des hypothèses derrières). néanmoins, je fais de la comparaison de modèles et il est possible que ce modèle soit ensuite réduit (je pensais à cela : GroupC * ShiftC + SexC). Par contre, n'étant pas experte, pouvez vous m'indiquer ce que signifierai le modèle que vous proposez? (je ne comprends le "^2")
  
  
• CR2 est l'exactitude ou non de la réponse. c'est codé en 0 et 1. Donc j'analyse d'abord cela, d'où le modèle binomial, non? le second modèle a ensuite pour but d'analyser le temps de réponse sur les réponses correctes uniquement (R2 correspond au temps de réponse en ms)
  
  
• Je recode les variables en -0.5, +0.5 car les paramètres sont plus faciles à interpréter ainsi. On nous avait rapidement évoqué cela en cours de stat et j'ai pris l'habitude de le faire. Après, ce que j'ai pu constater (sur le peu d'analyse stat que j'ai fait dans ma vie), c'est que lorsque je ne le faisais pas, j'avais plus facilement des pb de convergence et j'obtenais des résultats différents (ce que je ne comprends pas d'ailleurs car pour moi, ne pas recoder devrait changer uniquement la valeur des paramètres mais pas la significativité des résultats...)
  

EternalSunshine a écrit:auriez vous pare hasard une piste de lecture permettant d'expliquer pourquoi les groupes inégaux ne sont pas problématiques?

Non, je n'ai pas de piste à vous fournir. L'idée vient surtout que vous êtes en fait dans un modèle de régression, qui - lui - n'a pas de contrainte sur l'équilibre entre les modalités.

EternalSunshine a écrit:concernant l'interaction de niveau 3, c'est le modèle le plus complet (et j'ai pré-enregistré cette étude ainsi, avec des hypothèses derrières). néanmoins, je fais de la comparaison de modèles et il est possible que ce modèle soit ensuite réduit (je pensais à cela : GroupC * ShiftC + SexC). Par contre, n'étant pas experte, pouvez vous m'indiquer ce que signifierai le modèle que vous proposez? (je ne comprends le "^2")

La syntaxe que je vous ai proposée est un raccourci - dans R - pour avoir tous les effets principaux et toutes les interactions deux-à-deux. Très pratique.

EternalSunshine a écrit:CR2 est l'exactitude ou non de la réponse. c'est codé en 0 et 1. Donc j'analyse d'abord cela, d'où le modèle binomial, non? le second modèle a ensuite pour but d'analyser le temps de réponse sur les réponses correctes uniquement (R2 correspond au temps de réponse en ms)

Ok alors pour le premier modèle. Ca me parait correct. En revanche, pour le second modèle, vous analyser un temps, qui n'est donc pas binomial. Or votre syntaxe indique vous voulez ajuster un modèle binomial à vos données. Il y a problème, je pense.

EternalSunshine a écrit:Je recode les variables en -0.5, +0.5 car les paramètres sont plus faciles à interpréter ainsi. On nous avait rapidement évoqué cela en cours de stat et j'ai pris l'habitude de le faire. Après, ce que j'ai pu constater (sur le peu d'analyse stat que j'ai fait dans ma vie), c'est que lorsque je ne le faisais pas, j'avais plus facilement des pb de convergence et j'obtenais des résultats différents (ce que je ne comprends pas d'ailleurs car pour moi, ne pas recoder devrait changer uniquement la valeur des paramètres mais pas la significativité des résultats...)

C'est vous qui voyez, mais je trouve largement plus simple et facile à lire de coder vos facteurs qualitatifs, y compris binaires, par des textes "male", "female", etc.

HTH, Eric.

Je vous remercie à nouveau beaucoup pour vos réponses!

Très bien pour le point 1.

Merci pour la syntaxe! Je ne la connaissais pas donc je suis bien contente de la connaitre :-))

Concernant le second modèle, j'ai en effet fait une erreur de copié-collé, et j'avais copié 2 fois le même modèle sans m'en rendre compte. je comprends mieux votre incompréhension et j'ai corrigé.

Pour le dernier point, j'ai vraiment du mal à comprendre ce que fait R lorsqu'on ne recode pas et pourquoi je n'obtiens pas les mêmes résultats (au niveau des p-value).

Pare exemple, si j'ai ce modèle avec donc les VI recodées en -0.5, +0.5:

J'obtiens cela:

Code:

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: CR1 ~ GroupC * SexC + (1 | ID) + (1 | stim)
   Data: ESTdfWO

     AIC      BIC   logLik deviance df.resid
  4168.1   4208.3  -2078.0   4156.1     6044

Scaled residuals:
    Min      1Q  Median      3Q     Max
-9.2828  0.1554  0.2506  0.3732  1.6291

Random effects:
 Groups Name        Variance Std.Dev.
 ID     (Intercept) 0.3504   0.5919  
 stim   (Intercept) 1.3464   1.1603  
Number of obs: 6050, groups:  ID, 275; stim, 22

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   2.2775     0.2567   8.873   <2e-16 ***
GroupC       -0.2629     0.1130  -2.328   0.0199 *  
SexC          0.2859     0.1130   2.531   0.0114 *  
GroupC:SexC  -0.1961     0.2260  -0.868   0.3854    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) GroupC SexC  
GroupC       0.040              
SexC        -0.005  0.152      
GroupC:SexC  0.031 -0.043  0.196

Si je ne recode pas mes variables, avec le même modèle j'obtiens cela:

Code:

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: CR1 ~ group * sex + (1 | ID) + (1 | stim)
   Data: ESTdfWO

     AIC      BIC   logLik deviance df.resid
  4168.1   4208.3  -2078.0   4156.1     6044

Scaled residuals:
    Min      1Q  Median      3Q     Max
-9.2829  0.1554  0.2506  0.3732  1.6291

Random effects:
 Groups Name        Variance Std.Dev.
 ID     (Intercept) 0.3504   0.5919  
 stim   (Intercept) 1.3464   1.1603  
Number of obs: 6050, groups:  ID, 275; stim, 22

Fixed effects:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)    2.2400     0.2816   7.954  1.8e-15 ***
groupTD        0.3610     0.1562   2.310   0.0209 *  
sexM          -0.1879     0.1748  -1.075   0.2824    
groupTD:sexM  -0.1961     0.2259  -0.868   0.3853    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) gropTD sexM  
groupTD     -0.375              
sexM        -0.340  0.610      
groupTD:sxM  0.259 -0.692 -0.773

Pourquoi je n'obtiens pas les mêmes résultats? (là ca va dans le sens de plus d'effets significatifs, et je précise que ce n'est pas pour cela que je préfère l'autre approche, car c'est variable; mais j'aimerais bien comprendre pourquoi je ne trouve pas les mêmes effets)

Ca montre bien que votre recodage en -0.5 et 0.5 est foireux (selon moi).

Je le rappelle : Vous être dans une régression, en fait. Si vous ne recodez pas vos variables binaires, une modalité (par exemple ici sexe femelle) va prendre en fait la valeur 0 et le modèle va tester l'effet de l'ajout de l'autre modalité, qui va prendre la valeur de 1. Dans vos sorties, vous avez bien comme unique coefficient un truc intitulé 'SexM" qui correspond bien à ceci.

Si vous recodez, Vous allez faire passer au mieux votre regression entre des points ayant pour abscisses -0.5 ou 0.5 et le modèle va vous sortir un pente globale (qui - cette fois-ci - va s'intituler directement "SexC"). Les points (observations) ne sont plus aux mêmes endroits, et il est normal que le modèle et ne soit pas le même et les résultats également.

Ce raisonnement est le même pour la variable groupC, et pour l'interaction, etc.

J'espère que c'est plus clair ainsi.

HTH, Eric.

Merci à nouveau d'avoir pris le temps de me répondre.

En fait ce que je comprends (et qui me semble corroboré par ce qui est écrit sur cette page: http://talklab.psy.gla.ac.uk/tvw/catpred/), c'est que pour les effets principaux, mon modèle de régression va tester l'effet de ma VI quand l'autre VI est à 0. Du coup, je pense plutôt que c'est le fait de ne pas recoder qui est foireux, en tout cas par rapport à ce que je souhaite tester.

Je m'explique:

Avec le deviation coding (donc -0.5 +0.5), lorsque mon modèle estime la pente liée à l'effet Sex pour Group = 0, cela correspond à une valeur moyenne de Group. Ainsi, je teste l'effet principal de la variable Sex, indépendamment de la variable Group. C'est ce que je souhaite. Idem, lorsque je teste l'effet de Group (la pente pour l'effet Group), c'est indépendamment du fait que les participants soient des hommes ou des femmes, et c'est ce que je souhaite savoir.

Si à l'inverse je laisse mes variables telles qu'elles sont, lorsque mon modèle estime la pente liée à l'effetSex pour Group = 0, cela correspond à mon groupe de personnes autistes (puisque R code par ordre alphabétique et que j'ai un groupe ASD qui correspond à Autisme et un groupe TD correspondant à Typically Developped). Du coup, je ne teste plus un effet principal mais l'effet simple du Sex chez les personnes autistes... De même, l'effet du Group est en réalité l'effet simple de Group chez les femmes.

Pour vérifier si cela était juste, j'ai fait 2 autres recodages de la variable sexe.

1/ au lieu d'avoir Male à -0.5 et Female à 0.5, j'ai mis Female à -0.5 et Male à 0.5 pour faire à tourner le glm. Cela ne doit donc pas changer les paramètres de mon glm0, et en effet, ca ne change pas (à part que la pente change de signe ce qui est logique). (les résultats de ce modèle se réfèrent au premier que j'avais posté en dessus)

Code:

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: CR1 ~ GroupC * SexC2 + (1 | ID) + (1 | stim)
   Data: ESTdfWObis

     AIC      BIC   logLik deviance df.resid
  4168.1   4208.3  -2078.0   4156.1     6044

Scaled residuals:
    Min      1Q  Median      3Q     Max
-9.2829  0.1554  0.2506  0.3732  1.6291

Random effects:
 Groups Name        Variance Std.Dev.
 ID     (Intercept) 0.3504   0.5919  
 stim   (Intercept) 1.3463   1.1603  
Number of obs: 6050, groups:  ID, 275; stim, 22

Fixed effects:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)    2.2776     0.2567   8.874   <2e-16 ***
GroupC        -0.2629     0.1130  -2.328   0.0199 *  
SexC2         -0.2859     0.1130  -2.531   0.0114 *  
GroupC:SexC2   0.1961     0.2259   0.868   0.3853    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) GroupC SexC2
GroupC       0.040              
SexC2        0.005 -0.152      
GroupC:SxC2 -0.031  0.043  0.196

2/ Pour ce qui est du modèle non recodé, j'ai renommé Female en NM (Non Male), afin que ce soit la variable Male qui soit à 0 (au lieu de Female). ici, je devrais donc avoir un changement de mes paramètres pour le Groupe car celui ci correspondra à l'effet simple du groupe (la pente liée à l'effet Group) chez les Male (au lieu de Female). Et en effet, l'effet du groupe change alors que les effets du sexe et de l'interaction ne changent pas (les résultats de ce modèle se réfèrent au deuxième que j'avais posté en dessus).

Code:

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: CR1 ~ group * sex2 + (1 | ID) + (1 | stim)
   Data: ESTdfWObis

     AIC      BIC   logLik deviance df.resid
  4168.1   4208.3  -2078.0   4156.1     6044

Scaled residuals:
    Min      1Q  Median      3Q     Max
-9.2829  0.1554  0.2506  0.3732  1.6291

Random effects:
 Groups Name        Variance Std.Dev.
 ID     (Intercept) 0.3504   0.5919  
 stim   (Intercept) 1.3464   1.1603  
Number of obs: 6050, groups:  ID, 275; stim, 22

Fixed effects:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)      2.0522     0.2765   7.423 1.15e-13 ***
groupTD          0.1649     0.1632   1.010    0.312    
sex2NM           0.1879     0.1748   1.075    0.282    
groupTD:sex2NM   0.1961     0.2259   0.868    0.385    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) gropTD sex2NM
groupTD     -0.307              
sex2NM      -0.286  0.487      
grpTD:sx2NM  0.225 -0.722 -0.773

Je ne suis pas un très grand spécialiste du modèle mixte. Je pense cependant que votre interprétation est correcte pour chaque coefficient donné dans l'ajustement (et j'imagine que vous fondez votre raisonnement là-dessus). En revanche, si vous construisez votre table d'anova (avec la fonction anova() sur votre objet glmm), les tests que vous récupérez sont globaux, pour chaque variable explicative, et ne correspondent pas au test d'une variable pour un sous-groupe d'une autre.

Essayez, et vous verrez. Et merci de nous tenir informés.

Un dernier point : Pour ce qui est des glmm, la page de Ben Bolker reprend tout ce que l'on doit savoir et est indispensable, ici : https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html.

Eric.

Merci encore une fois pour vos réponses.

En effet, si je réalise l'anova de mon modèle, j'obtiens du coup la même chose pour les deux modèles:

Code:

> Anova(glm0Dum1)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: CR1
          Chisq Df Pr(>Chisq) 
group    5.4588  1    0.01947 *
sex      7.9610  1    0.00478 **
group:sex 0.8025  1    0.37034 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> Anova(glm0Cont1)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: CR1
            Chisq Df Pr(>Chisq) 
GroupC      5.4592  1    0.01947 *
SexC        7.9609  1    0.00478 **
GroupC:SexC 0.8024  1    0.37038 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

En tout cas, merci, cet échange m'a fait avancer je crois. Merci aussi pour le lien que j'avais trouvé dans mes recherches aussi et qui fait partie des sites que j'ai conservé précieusement!