Modèle mixte: construction des effets aléatoires

[Sébastien Nusslé]

Bonjour,

J'essaie de faire des analyses de variance avec un modèle à effets mixtes.

Le design est relativement simple, nous avons croisé des poissons (12 mâles et 4 femelles, toutes les combinaisons possibles) et élevé des œufs dans trois différents traitements (stress petit, moyen et élevé). Nos variables réponses sont: la survie de chaque œuf (variable binaire) ainsi que la date d'éclosion (variable normale).

Nous aimerions mesurer l'effet mâle (effet génétique) et l'effet femelle (effet génétique + investissement maternel) et regarder l'effet du traitement sur ces effets.

Premièrement: est-ce que la proportion de variance expliquée par l'effet mâle est différente selon les traitements (autrement dit est-ce l'effet génétique paternel est plus important si le stress est élevé).

Deuxièmement: est-ce que ce sont les mêmes mâles qui ont un effet positif (respectivement négatif) dans tous les traitements (on aimerait savoir si certains mâles ont un meilleur succès que les autres et si ces mâles sont les mêmes dans chaque traitement).

Problème: la variance des réponses est très différente selon les traitements de plus le nombre d'oeufs utilisés par combinaison parentale varie de 10 (traitement low stress) à 4 (traitement high stress).

Ma toute première question est la suivante:

- Est-ce que cela fait sens d'effectuer l'analyse en une seule fois ?

lmer(survie ~ traitement + (traitement | male) + (traitement | femelle), family=binomial)

- Ou faut-il analyser chaque traitement séparément ?
stress1 <- lmer(survie ~ 1 + (1 | male) + (1 | femelle), subset(data, stress=="low", family=binomial)
stress2 <- ...

La seconde question est: comment choisir (au niveau conceptuel) la matrice variance-covariance des effets aléatoires et comment le faire au niveau pratique.

Calcul d'une valeur moyenne pour chaque parent dans chaque stress différent:
lmer(survie ~ traitement + (traitement | male) + (traitement | femelle), family=binomial)
Calcul d'une valeur moyenne, puis d'une valeur pour chaque niveau de stress:
lmer(survie ~ traitement + (1 | male) + (traitement -1| male) + (1 | femelle) + (traitement-1 | femelle), family=binomial)
Calcul d'une valeur pour chaque parent dans chaque niveau de stress (idem que 1) mais avec des valeurs indépendantes les unes des autres:
???

Merci d'avance pour vos réponses !

[Renaud Lancelot]

Bonjour,

J'essaie de faire des analyses de variance avec un modèle à effets mixtes.

Le design est relativement simple, nous avons croisé des poissons (12 mâles et 4 femelles, toutes les combinaisons possibles) et élevé des œufs dans trois différents traitements (stress petit, moyen et élevé). Nos variables réponses sont: la survie de chaque œuf (variable binaire) ainsi que la date d'éclosion (variable normale).

Nous aimerions mesurer l'effet mâle (effet génétique) et l'effet femelle (effet génétique + investissement maternel) et regarder l'effet du traitement sur ces effets.

Premièrement: est-ce que la proportion de variance expliquée par l'effet mâle est différente selon les traitements (autrement dit est-ce l'effet génétique paternel est plus important si le stress est élevé).

La réponse étant binaire, on est dans le cadre du modèle linéaire généralisé à effets aléatoires (GLMM) ==> utiliser glmer avec une famille binomiale et un lien logit. Cela pose des pbs pour la mesure de l'héritabilité car la variance résiduelle et la variance des effets aléatoires ne sont pas sur la même échelle. Question récemment abordée sur ce forum, sans succès apparemment.

Je pense qu'il faudra utiliser des méthodes d'estimation MCMC pour construire des intervalles de confiance des quantités d'intérêt. A ma connaissance, pas de solution évidente avec R (la fct mcmcsamp du package lmer n'est pas opérationnelle pour les modèles GLMM). Il faudra probablement aller vers (Win)BUGS ou autres (MLwiN ?).

Deuxièmement: est-ce que ce sont les mêmes mâles qui ont un effet positif (respectivement négatif) dans tous les traitements (on aimerait savoir si certains mâles ont un meilleur succès que les autres et si ces mâles sont les mêmes dans chaque traitement).

Il faut extraire les effets aléatoires et les comparer aux traitements. Possible avec glmer.

Problème: la variance des réponses est très différente selon les traitements de plus le nombre d'oeufs utilisés par combinaison parentale varie de 10 (traitement low stress) à 4 (traitement high stress).

L'hétérogénéité de la variance est structurelle dans les modèles GLM(M). Par ailleurs, les modèles utilisés s'accomodent de plans déséquilibrés.

Ma toute première question est la suivante:

- Est-ce que cela fait sens d'effectuer l'analyse en une seule fois ?

lmer(survie ~ traitement + (traitement | male) + (traitement | femelle), family=binomial)

Ouch ! Je ne mettrais pas le traitement dans les effets aléatoires, et je n'utiliserais pas le sexe comme facteur de regroupement.

Je ferais:

Code : Tout sélectionner

glmer(survie ~ traitement * sexe + (1 | pere) + (1 | mere), family = binomial)

et je regarderais si les effets aléatoires du pere ou de la mere sont liés au traitement (graphique, test).

Renaud

[Sébastien Nusslé]

Tout d'abord, merci d'avoir pris un peu de temps pour me répondre, ensuite toute mes excuses pour mon manque de clarté. Je vais essayer de clarifier mon problème.

Ma toute première question est la suivante:

- Est-ce que cela fait sens d'effectuer l'analyse en une seule fois ?

lmer(survie ~ traitement + (traitement | male) + (traitement | femelle), family=binomial)

Ouch ! Je ne mettrais pas le traitement dans les effets aléatoires, et je n'utiliserais pas le sexe comme facteur de regroupement.

- Les facteurs males et femelles ici ne sont pas les sexe mais les identités des parents.
- Inclure le traitement dans les effet mixtes permet de spécifier une structure. Avec (1 | pere) l'effet aléatoire se fait sur l'intercept et donc calcule un effet aléatoire (ranef()) moyen pour chaque mâle, indépendamment du traitement.
- Avec (traitement | male), c'est l'équivalent d'un effet aléatoire sur la pente, sauf qu'ici la variable traitement est catégorielle. On calcule d'abord les effets aléatoires sur l'intercept (donc le premier niveau de stress avec les "dummy variables") et ensuite les effets individuels pour chacun des deux autres niveau de stress.
- Le problème c'est que par défaut la fonction met une structure à mes variables aléatoires et que les résultats ne sont pas du tout pareils si je change cette structure. De plus je la change en bricolant, mais je ne sais pas comment faire pour être rigoureux.
- Ce que j'aimerais faire c'est calculer un effet aléatoire niché dans le traitement pour voir si l'effet individuel d'un père diffère selon le niveau de stress. On s'attend en effet à ce que d'autres gènes soient activés lors de stress important, ce qui donnerait un avantage (que l'on résume par les ranef) à d'autres mâles que ceux qui sont habituellement favorisés.

Je ferais:

Code : Tout sélectionner

glmer(survie ~ traitement * sexe + (1 | pere) + (1 | mere), family = binomial)

et je regarderais si les effets aléatoires du pere ou de la mere sont liés au traitement (graphique, test).

Le problème ici c'est qu'un seul effet par mâle est calculé.

La réponse étant binaire, on est dans le cadre du modèle linéaire généralisé à effets aléatoires (GLMM) ==> utiliser glmer avec une famille binomiale et un lien logit.

Je n'ai pas encore compris la logique de glmer(), la fonction lmer a toujours pu, en spécifiant l'argument "family", faire des modèles linéaires généralisés. J'ai toujours fait mes régressions logistiques avec lmer... Est-ce que Douglas Bates a retiré cette possibilité de lmer pour la mettre dans glmer ? Quel est l'avantage d'utiliser glmer par rapport à lmer si les fonctions habituelles (mcmcsamp par exemple) ne sont pas encore implémentées ? Je n'ai rien vu dans le R news depuis la présentation de la fonction lmer en 2005.

Cela pose des pbs pour la mesure de l'héritabilité car la variance résiduelle et la variance des effets aléatoires ne sont pas sur la même échelle. Question récemment abordée sur ce forum, sans succès apparemment.

J'attends aussi avec impatience une solution, pour l'instant j'essaie de bricoler avec des permutations, mais cela ne permet pas de tout faire.

Problème: la variance des réponses est très différente selon les traitements de plus le nombre d'oeufs utilisés par combinaison parentale varie de 10 (traitement low stress) à 4 (traitement high stress).

L'hétérogénéité de la variance est structurelle dans les modèles GLM(M). Par ailleurs, les modèles utilisés s'accomodent de plans déséquilibrés.

Le problème, c'est que j'utilise une régression logistique, donc avec une distribution binomiale, et la variance devrait être maximale quand la probabilité de décès est de 0.5 est minimale quand on s'approche des extrêmes. Dans cette expérience ma variance augmente avec la mortalité. J'avais aussi noté que lmer n'avait pas trop de problèmes avec des plans déséquilibrés, mais je me demandais si la précision des estimations était diminuée et surtout si 4 oeufs pour une combinaison parentale était suffisante pour estimer précisément une probabilité de décès.

Encore merci

Sébastien

[Renaud Lancelot]

Je n'ai pas encore compris la logique de glmer(), la fonction lmer a toujours pu, en spécifiant l'argument "family", faire des modèles linéaires généralisés. J'ai toujours fait mes régressions logistiques avec lmer... Est-ce que Douglas Bates a retiré cette possibilité de lmer pour la mettre dans glmer ? Quel est l'avantage d'utiliser glmer par rapport à lmer si les fonctions habituelles (mcmcsamp par exemple) ne sont pas encore implémentées ? Je n'ai rien vu dans le R news depuis la présentation de la fonction lmer en 2005.

glmer a été introduit il y a qque temps, avec nlmer et gnlmer. Ce sont juste des interfaces (regarder le code de glmer) et pas de différence (pour le moment) avec lmer. Pas sûr que ce soit tjs le cas à l'avenir ==> autant utiliser les fonctions correspondant aux modèles qu'on veut effectivement ajuster.

Renaud

[Sébastien Nusslé]

Je me suis concentré sur les données avec le jour d'éclosion, car elles sont plus ou moins normales et les méthodes statistiques sont plus légères.

J'ai donc testé deux structures d'effets mixtes qui donnent des résultats assez différents, mais je ne suis pas certain de comprendre ce que fait la fonction.

J'ai l'impression que dans la première analyse, les effets aléatoires sont indépendants, alors qu'il ne le sont pas dans la seconde. Je me base sur la variance des effets aléatoires qui sont plus où moins égales dans la seconde analyse alors qu'ils sont liés à la variance des observations dans la première.

Deux questions:
Est-ce que quelqu'un sait comment la fonction traite la structure des effets aléatoires ?

Quel est votre avis sur la pertinence de l'un où l'autre des modèles quant aux questions posées, a savoir est-ce que les effets génétiques (paternels) sont différents selon les traitements.

Merci d'avance

Code : Tout sélectionner

Linear mixed model fit by REML
Formula: days_until_hatching ~ treatment + (treatment - 1 | Female) +      (treatment - 1 | Male)
   Data: data1
  AIC  BIC logLik deviance REMLdev
 4294 4368  -2131     4266    4262
Random effects:
 Groups   Name            Variance Std.Dev. Corr         
 Male     treatmenthigh   12.66886 3.5593               
          treatmentlow     0.52781 0.7265   0.185       
          treatmentmedium  1.00801 1.0040   0.528  0.932
 Female   treatmenthigh    4.69515 2.1668               
          treatmentlow     1.44648 1.2027   -0.017       
          treatmentmedium  2.00711 1.4167    0.215  0.973
 Residual                 19.28760 4.3918               
Number of obs: 728, groups: Male, 12; Female, 4
 
Fixed effects:
                Estimate Std. Error t value
(Intercept)       61.261      1.581   38.75
treatmentlow       4.529      1.699    2.66
treatmentmedium    4.552      1.600    2.85
 
Correlation of Fixed Effects:
            (Intr) trtmntl
treatmentlw -0.920       
treatmntmdm -0.860  0.963

Code : Tout sélectionner

Linear mixed model fit by REML
Formula: days_until_hatching ~ treatment + (treatment | Female) + (treatment |      Male)
   Data: data1
  AIC  BIC logLik deviance REMLdev
 4294 4368  -2131     4266    4262
Random effects:
 Groups   Name            Variance Std.Dev. Corr         
 Male     (Intercept)     12.6690  3.5594               
          treatmentlow    12.2406  3.4987   -0.979       
          treatmentmedium  9.9009  3.1466   -0.963  0.998
 Female   (Intercept)      4.6950  2.1668               
          treatmentlow     6.2278  2.4955   -0.876       
          treatmentmedium  5.3806  2.3196   -0.803  0.991
 Residual                 19.2876  4.3918               
Number of obs: 728, groups: Male, 12; Female, 4
 
Fixed effects:
                Estimate Std. Error t value
(Intercept)       61.261      1.581   38.75
treatmentlow       4.529      1.699    2.66
treatmentmedium    4.552      1.600    2.85
 
Correlation of Fixed Effects:
            (Intr) trtmntl
treatmentlw -0.920       
treatmntmdm -0.860  0.963