Mixed Model (lme dans nlme ou lme4)

Questions sur les fonctions statistiques de R

Modérateur : Groupe des modérateurs

Karine Poitrineau
Messages : 18
Enregistré le : 30 Mai 2006, 19:37

Mixed Model (lme dans nlme ou lme4)

Messagepar Karine Poitrineau » 31 Mai 2006, 13:46

Bonjour à tous,

J'ai quelques problèmes à comprendre quel code utiliser pour un modèle mixte avec la fonction lme (library nlme ou lme4).

Dans mon modèle mixte j'ai deux populations, dans ces populations plusieurs souches (nested dans populations), et pour chaque souches plusieurs répétitions (4 normalement s'il n'y a pas de données manquantes).
Au départ je souhaitais faire simplement une ANOVA hiérarchique pour une variable (continue) que j'ai mesurée, mais il s'avère que celle-ci est corrélée à une autre variable (continue aussi), sur l'ensemble de mon expérience.
Je souhaite donc intégrer cette variable dans mon analyse. A noter que cette variable n'est pas significativement différente entre populations, mais l'est entre les souches. Pour cette variable, je ne sais pas trop comment gérer la chose: si je ne me trompe pas, je devrais la considérer comme nested dans la variable "souche"?

Remarque: Comme je fais une analyse d'héritabilité de mon caractère, j'aurais besoin de connaître la variance existant entre souches pour la comparer à la variance intrasouches (en tenant compte de la variabilité de mon second facteur).
J'ai une variable binomiale (type survie-mort) pour un certain nombre d'individus de chaque répétition, ce qui me donne des pourcentages (c'est ma variable dépendante, que j'ai normalisée avec la formule arcsin). Je pourrais apparemment utiliser au lieu d'un modèle linéaire un "logistic mixed model" (ce qui évite la transformation des données), mais si j'ai bien compris la documentation que j'ai trouvée, ça ne me permet pas de faire d'analyse d'héritabilité...

Bref, j'espère avoir été claire, et que mon problème n'est pas trop compliqué (en même temps, l'idée me paraît simple, mais je ne trouve pas de cas ressemblant assez au mien pour être certaine de la façon de procéder). Je vais me procurer le Livre de Pinheiro et Bates (mixed-effects models in S and S-Plus, mais il va falloir que j'attende le transfert inter-BU ;) ).

Merci d'avance
K.ontente de pouvoir discuter avec d'autres utilisateurs de R

Renaud Lancelot
Messages : 2484
Enregistré le : 16 Déc 2004, 08:01
Contact :

Messagepar Renaud Lancelot » 31 Mai 2006, 14:57

Bonjour Karine,

J'ai quelques problèmes à comprendre quel code utiliser pour un modèle mixte avec la fonction lme (library nlme ou lme4).


Bienvenue au club ;-) !

Dans mon modèle mixte j'ai deux populations, dans ces populations plusieurs souches (nested dans populations), et pour chaque souches plusieurs répétitions (4 normalement s'il n'y a pas de données manquantes).


Si vous n'avez pas de données manquantes, vous pouvez utiliser la fonction aov en spécifiant un terme d'erreur (effet aléatoire). Voir l'aide de ?aov et les exemples à la fin.

Si il y a des données manquantes, il vaut mieux utiliser lme (package nlme) ou lmer (packages Matrix / lme4) car le modèle sous-jacent est adapté aux plans déséquilibrés (cf l'article de Laird et Ware). Il n'en demeure pas moins qu'il est indispensable de vérifier si les données manquantes ou perdus de vue se produisent de manière complètement aléatoire par rapport à la variable dépendante. Vous trouverez une méthode d'évaluation dans:

Diggle, P.J., Liang, K.-Y., Zeger, S.L., 1994. Analysis of longitudinal data, 1st Edition. Clarendon Press, Oxford, 253 p.

(il y a une édition plus récente).

Au départ je souhaitais faire simplement une ANOVA hiérarchique pour une variable (continue) que j'ai mesurée, mais il s'avère que celle-ci est corrélée à une autre variable (continue aussi), sur l'ensemble de mon expérience.
Je souhaite donc intégrer cette variable dans mon analyse. A noter que cette variable n'est pas significativement différente entre populations, mais l'est entre les souches. Pour cette variable, je ne sais pas trop comment gérer la chose: si je ne me trompe pas, je devrais la considérer comme nested dans la variable "souche"?


Il faudrait savoir si la variation aléatoire liée à la souche est suffisante pour rendre compte de l'effet de la seconde variable sans que celle-ci ne soit incluse dans les effets fixes du modèle. Cela peut se vérifier en regardant la corrélation des résidus de ce modèle avec cette variable. Une autre manière d'aborder le problème serait de faire le contraire: subsiste-t-il une variation liée à la souche si on inclut la seconde variable dans le modèle ? En d'autres termes, l'effet de la souche est-il attribuable à l'effet de cette seconde variable (ou vice versa) ? Si oui, et si l'objectif du modèle est de faire de la prévision, je laisserais tomber la variable la plus coûteuse à observer: souche ou seconde variable. Sinon...

Remarque: Comme je fais une analyse d'héritabilité de mon caractère, j'aurais besoin de connaître la variance existant entre souches pour la comparer à la variance intrasouches (en tenant compte de la variabilité de mon second facteur).


C'est un modèle des composantes de la variance. Matthieu Lesnoff a fait une fiche sur ce sujet, disponible sur ce forum: exemple traité avec lme.

J'ai une variable binomiale (type survie-mort) pour un certain nombre d'individus de chaque répétition, ce qui me donne des pourcentages (c'est ma variable dépendante, que j'ai normalisée avec la formule arcsin). Je pourrais apparemment utiliser au lieu d'un modèle linéaire un "logistic mixed model" (ce qui évite la transformation des données), mais si j'ai bien compris la documentation que j'ai trouvée, ça ne me permet pas de faire d'analyse d'héritabilité...


NB: la "normalisation" se fait avec la transformation arcsin(sqrt(p)) et non arcsin(p). Il faut voir quelle est la valeur de p, et l'ordre de grandeur du dénominateur. Dans certains cas, un modèle linéaire sur les probas peut être suffisant (probas pas trop hautes ou basses, grand dénominateur).

Avec un modèle linéaire généralisé mixte, le problème est que la variance résiduelle et la variance aux niveaux supérieurs (souche, population) ne sont pas sur la même échelle, ce qui rend les comparaisons difficiles. Il faudrait faire la biblio pour voir les évolutions récentes. Il y a peut-être des solutions avec des approches bayésiennes, elles-mêmes d'ailleurs disponibles dans lme4 (voir fonction mcmcsamp).

Bref, j'espère avoir été claire, et que mon problème n'est pas trop compliqué (en même temps, l'idée me paraît simple, mais je ne trouve pas de cas ressemblant assez au mien pour être certaine de la façon de procéder). Je vais me procurer le Livre de Pinheiro et Bates (mixed-effects models in S and S-Plus, mais il va falloir que j'attende le transfert inter-BU ;) ).


C'est un bon point de départ. La fonction lmer est plus puissante que lme en cas de structure aléatoire complexe, mais elle est moins bien carrossée et son aide est rachitique. Voir toutefois le package mlmRev qui traite des exemples d'analyses hiérarchiques avec lmer.

Amicalement,

Renaud

Karine Poitrineau
Messages : 18
Enregistré le : 30 Mai 2006, 19:37

Messagepar Karine Poitrineau » 31 Mai 2006, 16:32

Bonjour,
Merci beaucoup pour votre réponse :)

Si vous n'avez pas de données manquantes, vous pouvez utiliser la fonction aov en spécifiant un terme d'erreur (effet aléatoire). Voir l'aide de ?aov et les exemples à la fin.

Malheureusement j'ai des données manquante (rien n'est parfait, surtout en biologie!)...

Si il y a des données manquantes, il vaut mieux utiliser lme (package nlme) ou lmer (packages Matrix / lme4) car le modèle sous-jacent est adapté aux plans déséquilibrés (cf l'article de Laird et Ware). Il n'en demeure pas moins qu'il est indispensable de vérifier si les données manquantes ou perdus de vue se produisent de manière complètement aléatoire par rapport à la variable dépendante.

Ah, les choses se compliquent alors!
Il est bien probable que j'aie plus de donnée manquantes lorsque ma seconde variable est faible pour une souche : cette seconde variable est un degré d'infestation par des parasitoïdes, alors que ma variable dépendante est un taux de résistance aux parasitoïdes. Hors, dans certains cas je n'ai pas eu d'infestation (donc pas de résistance possible), ce qui entraîne évidemment des données manquantes.

Il faudrait savoir si la variation aléatoire liée à la souche est suffisante pour rendre compte de l'effet de la seconde variable sans que celle-ci ne soit incluse dans les effets fixes du modèle.

En fait, dans une ANOVA hiérarchique "simple" (sans inclure ma seconde variable) je n'ai pas d'effet visible ni de la population, ni de la souche. Donc la différence de degré d'infestation entre souche ne se répercute pas visiblement sur les différences de résistance entre les souches (malgré une corrélation globale significative).
C'est d'ailleurs pour ceci que je serais intéressée par introduire le degré d'infestation dans mon analyse, car je suppose que sa variabilité intra-souche pourrait masquer celle de la résistance entre souches (je ne sais pas si je suis claire et si mon raisonnement est correct... je commence même de plus en plus à me dire qu'il est peut-être bien foireux :-( :oops: ).

Une autre manière d'aborder le problème serait de faire le contraire: subsiste-t-il une variation liée à la souche si on inclut la seconde variable dans le modèle ? En d'autres termes, l'effet de la souche est-il attribuable à l'effet de cette seconde variable (ou vice versa) ?

Cf avant: je me posais la question inverse : "pourrait-on observer une variation liée à la souche qu'on n'observe pas autrement, si on introduit l'effet de la seconde variable?"

Si oui, et si l'objectif du modèle est de faire de la prévision, je laisserais tomber la variable la plus coûteuse à observer: souche ou seconde variable. Sinon...

Mon objectif n'est pas de faire un modèle prédictif, mais de mettre en évidence les effets (ou l'absence d'effet) de la population (apparemment nul), du degré d'infestation et de la souche.

NB: la "normalisation" se fait avec la transformation arcsin(sqrt(p)) et non arcsin(p). Il faut voir quelle est la valeur de p, et l'ordre de grandeur du dénominateur. Dans certains cas, un modèle linéaire sur les probas peut être suffisant (probas pas trop hautes ou basses, grand dénominateur).

J'ai effectivement fait ce traitement (arcsin transformation = angular transformation), car j'ai des valeurs de p qui peuvent être très variables, notamment très faibles.

En tout cas merci de vous pencher sur mon problème (j'apprends des choses :) ).

K.

Renaud Lancelot
Messages : 2484
Enregistré le : 16 Déc 2004, 08:01
Contact :

Messagepar Renaud Lancelot » 31 Mai 2006, 17:10

Karine Poitrineau a écrit :Bonjour,
Merci beaucoup pour votre réponse :)

Si vous n'avez pas de données manquantes, vous pouvez utiliser la fonction aov en spécifiant un terme d'erreur (effet aléatoire). Voir l'aide de ?aov et les exemples à la fin.

Malheureusement j'ai des données manquante (rien n'est parfait, surtout en biologie!)...

Si il y a des données manquantes, il vaut mieux utiliser lme (package nlme) ou lmer (packages Matrix / lme4) car le modèle sous-jacent est adapté aux plans déséquilibrés (cf l'article de Laird et Ware). Il n'en demeure pas moins qu'il est indispensable de vérifier si les données manquantes ou perdus de vue se produisent de manière complètement aléatoire par rapport à la variable dépendante.

Ah, les choses se compliquent alors!


Oui... Vous trouverez la présentation d'un exemple de diagnostic des données manquantes dans ma thèse, disponible en ligne.

Il est bien probable que j'aie plus de donnée manquantes lorsque ma seconde variable est faible pour une souche : cette seconde variable est un degré d'infestation par des parasitoïdes, alors que ma variable dépendante est un taux de résistance aux parasitoïdes. Hors, dans certains cas je n'ai pas eu d'infestation (donc pas de résistance possible), ce qui entraîne évidemment des données manquantes.


Il faut probablement reformuler les questions posées, et peut-être séparer les analyses: probabilité d'infestation par les parasitoïdes, puis résistance chez les infestés. A discuter avec vos collègues biologistes.

Il faudrait savoir si la variation aléatoire liée à la souche est suffisante pour rendre compte de l'effet de la seconde variable sans que celle-ci ne soit incluse dans les effets fixes du modèle.

En fait, dans une ANOVA hiérarchique "simple" (sans inclure ma seconde variable) je n'ai pas d'effet visible ni de la population, ni de la souche. Donc la différence de degré d'infestation entre souche ne se répercute pas visiblement sur les différences de résistance entre les souches (malgré une corrélation globale significative).
C'est d'ailleurs pour ceci que je serais intéressée par introduire le degré d'infestation dans mon analyse, car je suppose que sa variabilité intra-souche pourrait masquer celle de la résistance entre souches (je ne sais pas si je suis claire et si mon raisonnement est correct... je commence même de plus en plus à me dire qu'il est peut-être bien foireux :-( :oops: ).


Oui, je pense qu'il y a un problème de logique. Une manière d'y voir + clair est de faire un pré-modèle conceptuel d'analyse: schéma résumant les hypothèses de causalité entre variables dépendantes et explicatives (voir les bouquins de Frontier sur l'échantillonnage). Même sans aller jusque là, je ne crois pas que vous puissiez mélanger dans la même analyse des observations sur des populations non comparables. Cf ci-dessus.


Une autre manière d'aborder le problème serait de faire le contraire: subsiste-t-il une variation liée à la souche si on inclut la seconde variable dans le modèle ? En d'autres termes, l'effet de la souche est-il attribuable à l'effet de cette seconde variable (ou vice versa) ?

Cf avant: je me posais la question inverse : "pourrait-on observer une variation liée à la souche qu'on n'observe pas autrement, si on introduit l'effet de la seconde variable?"


Dit autrement, cette seconde variable définit des strates de populations, et la première chose à faire (à mon avis...) est une analyse intra-strate.

Si oui, et si l'objectif du modèle est de faire de la prévision, je laisserais tomber la variable la plus coûteuse à observer: souche ou seconde variable. Sinon...

Mon objectif n'est pas de faire un modèle prédictif, mais de mettre en évidence les effets (ou l'absence d'effet) de la population (apparemment nul), du degré d'infestation et de la souche.

NB: la "normalisation" se fait avec la transformation arcsin(sqrt(p)) et non arcsin(p). Il faut voir quelle est la valeur de p, et l'ordre de grandeur du dénominateur. Dans certains cas, un modèle linéaire sur les probas peut être suffisant (probas pas trop hautes ou basses, grand dénominateur).

J'ai effectivement fait ce traitement (arcsin transformation = angular transformation), car j'ai des valeurs de p qui peuvent être très variables, notamment très faibles.

Si ces valeurs sont effectivement très faibles, on sort du cadre où la transformation arcsin(sqrt(p)) est pertinente (cela se voit en étudiant les résidus). Il faudrait alors s'orienter vers des modèles binomiaux ou de Poisson. Peut-être pourriez-vous:

(1) dans un premier temps reformuler vos objectifs d'analyse et en déduire le ou les jeux de données à analyser,

(2) faire une bonne analyse exploratoire pour identifier les sources de variations importantes,

(3) utiliser un modèle marginal (plus "facile" que les modèles mixtes) pour estimer l'effet des facteurs importants: modèle bétabinomial pour des proportions (si p pas trop petit) ou binomial négatif sur comptages (si p petit: < 1% par exemple): méthodes disponibles dans le package aod, ou GEE,

(4) modèle mixte pour séparer les sources de variance aléatoire (souche, pop).

En tout cas merci de vous pencher sur mon problème (j'apprends des choses :) ).

K.

Karine Poitrineau
Messages : 18
Enregistré le : 30 Mai 2006, 19:37

Messagepar Karine Poitrineau » 31 Mai 2006, 19:09

Si ces valeurs sont effectivement très faibles, on sort du cadre où la transformation arcsin(sqrt(p)) est pertinente (cela se voit en étudiant les résidus). Il faudrait alors s'orienter vers des modèles binomiaux ou de Poisson.

Plus je retourne le problème, plus je me dis que c'est effectivement un des "os" de mon problème...
Bon je n'ai pas tout expliqué pour ne pas surcharger mes explications, mais dans mon analyse, j'avais trois caractères, qui sont des pourcentages (degré d'infestation, développement des parasitoïdes, et résistance). Ces trois caractères se sont avérés être corrélés entre eux, et les analyses montrent une variation significative des deux premiers caractères entre les souches... mais pas du dernier (ce qui était frustrant). Or, les fameux pourcentages de résistance sont très faibles...
A l'époque (ce sont de vieilles expériences de DEA) quand j'avais présenté mes analyse à mon professeur de statistiques (que vous connaissez ;-) ), il m'avait dit que la méthode d'arcsin-transformation était "grossière", qu'il y avait des méthodes plus élégantes, mais que l'ANOVA était une méthode robuste, donc ça "pouvait passer"...

Bref, je vais laisser tomber l'idée de calculer des héritabilités (qui de toute façon ne signifient pas grand chose car ce sont des "héritabilités au sens large") et essayer de trouver le moyen de faire les analyses de façon un peu plus propre en considérant le caractère binomial de mes données.

Remarque: là où ça se corse, c'est qu'on avait fait une analyse de corrélation canonique entre les trois pourcentages et trois autres caractères de mes lignées, qui eux étaient quantitatifs... certainement pas si facile à refaire...

K.


Retourner vers « Archives : Fonctions statistiques »

Qui est en ligne

Utilisateurs parcourant ce forum : Aucun utilisateur enregistré et 1 invité

cron