[RESOLU] Mesures répétées, lme et random effects

[Simon NAVEL]

Tout d'abord bonjour à tous.

J'éprouve d'énormes difficultés à analyser les résultas obtenus d'un design exprimental un peu complexe.
Je me suis dirigé (suite à des discussions) vers les lme, sans avoir eu de formations sur ce domaine là de la Statistique.
Je m'excuse d'ailleur à l'avance si toutefois ma question vous semble idiote.


Le design est le suivant:

- 18 unités expérimentales (des colonnes verticales cylindriques, alimentées en eau depuis le haut) définies chaucune par:
1 porosité (facteur a 3 modalités)
1 traitement faunistique (variable "GAMMARE" à2 modalités: avec et sans faune)
-> (Donc des triplicats de colonne pour chaque combinaison).

- Pour chacune de ces colonnes: nous réalisons des prélevements d'eau à 3 horizons (=profondeurs fixées). "Horizon" est donc niché dans "Colonne".

- Cette opération est répétée 1 fois par semaine pendant 5 semaines. (1 unique prélèvement par semaine, pas de réplicats).



Nous mesurons à la fois des teneurs en oxygène, nitrate, ammonium et phosphate. Cependant dans un premier temps je souhaiterai les analyser séparément.

En prenant l'exemple du nitrate (NO3), la question que je me pose est simple: les variables "porosite" et "faune" ont-elles un impact significatif sur les teneurs mesurées?



Après discussion avec des collègues, il m'a été conseillé d'étudier le livre de Pinheiro & Bates : "Mixed effects Model in S and S-Plus".
Je l'ai lu et relu (j'avoue ne pas avoir tout bien compris) mais je reste coincé dans ma façon de coder les random effects.

J'aimerais regarder l'impact des 2 variables explicatives (POROSITE et FAUNE), après m'être affranchi de:
-la part de variabilité expliquée par les variations temporelles
-la variablité due à l'horizon étudié au sein de la colonne.
qui potentiellement estompe les effets des facteurs explicatifs.


J'ai aboutit à un modèle de la sorte:

Code : Tout sélectionner

lme2=lme(fixed=NO3~POROSITE*GAMMARE,data=analyseNO3,random=list(COLONNE=~1,HORIZON=~1))

Seulement je suis bloqué à cette étape, et ne vois pas comment faire en sorte d"enlever" la part de variablité résiduelle due aux variablités dans le temps .

R se fige quand je lui envoie un modele du type:

Code : Tout sélectionner

lme1=lme(fixed=NO3~POROSITE*GAMMARE,data=analyseNO3,random=list(COLONNE=~DATE,HORIZON=~DATE))

Malgré la longueur du message, ma requête est au final assez courte: quelqu'un saurait-il comment faire en sorte de me "séparer" de cette variablité des mesures répétées dans le temps? (pourquoi R ne veut pa appliquer le modele dernierement employé?)[/code]

En espérant ne pas dire de betises plus grosses que moi et en m'excusant encore si toutefois ce devait être le cas.

Simon

[Pierre Bady]

bonjour,

J'éprouve d'énormes difficultés à analyser les résultas obtenus d'un design exprimental un peu complexe.
Je me suis dirigé (suite à des discussions) vers les lme, sans avoir eu de formations sur ce domaine là de la Statistique.
Je m'excuse d'ailleur à l'avance si toutefois ma question vous semble idiote.

si vous ne maîtrisez pas le sujet, il est souvent plus prudent et plus efficace d’aller voir un biostateux (ou un stateux, à Lyon, ce n’est pas cela qui manque) qui maîtrise le sujet. Souvent, cela permet de gagner du temps, d’apprendre plein de chose … et des fois, ça permet même de mettre en place des collaborations.

Après discussion avec des collègues, il m'a été conseillé d'étudier le livre de Pinheiro & Bates : "Mixed effects Model in S and S-Plus". Je l'ai lu et relu (j'avoue ne pas avoir tout bien compris) mais je reste coincé dans ma façon de coder les random effects (pourquoi R ne veut pa appliquer le modele dernierement employé?)

R fait bêtement ce qu’on lui demande, ses erreurs sont nos erreurs ...(la plupart du temps)
La formulation que vous utilisez n’est pas correcte.

- avec lme, on peut considérer l’autocorrélation temporelle via des matrices corrélations spécifiques (ex. corAR1, viewtopic.php?t=319&highlight=lme )

- on pourrait également envisager une approche avec des modèles marginaux (gee, cf par exemple viewtopic.php?p=518 ).

- il y a également des alternatives robustes (cf les packages sandwich, lmtest, …).

- et enfin, une approche baysienne de la question pourrait être envisagée.

Mais bon, dans tous les cas, ce n’est pas simple à mettre en œuvre …

Personnellement*, je vous conseille d’aller voir des gens compétents dans le domaine (pourquoi pas les gens qui vous ont conseillé le bouquin ?). De plus, une fois le modèle fait … il reste le plus important à faire (e.g. analyser la structure des résidus, les points leviers, les valeurs influentes, etc …un modèle ne se résume pas à quelques étoiles).



hth


pierre





liens utiles :
viewtopic.php?t=245
http://www.ats.ucla.edu/stat/books/default.htm (section Longitudinal Data and Repeated Measures)
http://cran.r-project.org/doc/contrib/F ... models.pdf
* enfin, c’est juste mon avis, ça vaut ce que ça vaut.

[Simon NAVEL]

Merci pour ces recommandations

Simon

[Simon NAVEL]

J'ai simplifier mon jeu de données.
L'effet profondeur a été supprimé, je ne résonne désormais que sur la différences des 2 relevés les plus en profondeur.

Voici le détail du jeu de données si cela peux éclaircir mes propos:
J'ai 18 mesocosmes (COLONNE) caractérisés par
* 1 variable explicative POROSITE "fixe" à 3 valeurs (1, 2 ou 3)
* 1 variable explicative GAMMARE "fixe" à 2 valeurs (1 ou 0)
(soit 3 mesocosmes "identiques" pour chaque croisement possible de ces 2 variables)
Dans ces mesocosmes je réalise une mesure "02" par semaine (soit donc 18 valeurs par semaine)
Enfin je répète la mesure 4 fois, sur les mêmes colonnes (variable DATE de T1 à T4).
Ce qui fait 72 valeurs en tout.

Code : Tout sélectionner

  DATE COLONNE POROSITE GAMMARE     O2
1   T1       1        1       1 1.0731
2   T1       2        1       0 1.0512
3   T1       3        1       1 1.0731
4   T1       4        1       0 0.9855
5   T1       5        1       1 1.0950
6   T1       6        1       0 1.0293

Mon problème est le suivant:
Dans le cadre d'analyses de mesures répétées, j'aimerais pouvoir spécifier dans R, le fait que toutes les interactions avec le random affect (DATE) sont elles aussi "random" (que se soit en utilisant l'approche aov avec Error() ou l'approche lme).


J'ai d'abord adopté l'approche suivante sous R:

Code : Tout sélectionner

summary(aov(O2~as.factor(POROSITE)*as.factor(GAMMARE)+Error(DATE))

En faisant une anova répétée sur un logiciel "presse bouton" (Statistica), j'obtiens exactement les mêmes valeurs pour ce qui est des SumSq et ddl pour les effets fixes, qu'avec l'approche précédente sous R. Cependant le calcul du F est faussé par le SumSq et ddl résiduel qui sont très différents entres les 2 logiciels.
Dans R, le F est calculé avec un résidus calculé avec 63ddl (72 - 1 (intercept) -2 (effet POROSITE) -1 (effet GAMMARE) -2(interaction)).
Dans Statistica, j'ai un résidus avec 12ddl (qui correspondent à 2 ddl pour chaque croisement d'effet fixe possible (soit 2x6)..ici le résidus correspond donc bien aux simples variabilités "intra groupe")

Cette différence de ddl viend du fait que R n'a pas pris comme random effects (et donc exclu du residu servant à calculer les F) les interactions entre DATE et les 2 variables explicatives fixes (DATE:P; DATE:G; DATE:G:P), ce que Statistica fait automatiquement (a priori les interactions avec un effet alétoire sous censées être elles même aléatoires).

Après consultation auprès d'un "spécialiste modèle mixte" sur Lyon, nous avons conclut qu'effectivement je devais n'avoir que 12 ddl résiduels pour calculer les F des effets fixes. Nous avons aboutit à la formulation suivante, en utilisant les lme (le but étant donc d'avoir l'équivalent sous R de ce qu'on obtient sous Statistica).

Code : Tout sélectionner

lme1=lme(fixed=O2~as.factor(POROSITE)*as.factor(GAMMARE),
    random=~as.factor(POROSITE)*as.factor(GAMMARE)|as.factor(DATE),
    method="ML")

Ainsi je retirerai des ddl dans le residus servant a calculer les F et pvalue des variables explicatives (à savoir tous les ddl correspondant aux interactions avec la variable DATE).

Seulement R m'affiche le message d'erreur suivant:

Code : Tout sélectionner

Erreur dans lme.formula(fixed = O2 ~ as.factor(POROSITE) * as.factor(GAMMARE),  :
  nlminb problem, convergence error code = 1
  message = iteration limit reached without convergence (9)

J'ai parcouru les forums et j'ai compris que cela devait venir des données, et non du code. J'ai essayé d'utiliser la méthode optim au lieu de nlmimb, avec également un nombre bien plus important d'iterations, mais rien n'y fait.
Dès que je mets un random plus complexe qu'un simple effet DATE (~1|DATE), je n'arrive plus à obtenir la convergence

Le jeu de données m'empêche t-il tout emploi du lme pour spécifier des effets aléatoires complexes, ou y a t'il une solution à laquelle je n'ai pas pensé...me permettant d'enfin intégrer les interactions avec le temps comme étant également "aléatoires"?


En m'excusant pour la longueur du texte et en vous remerciant par avance pour votre aide.

Simon NAVEL

[Eric Pagot]

Pour moi, le facteur d'erreur serait plutôt la colonne. Je travaille avec systat et la présentation des données se ferait de la façon suivante :
Colonne Porosité Gammare MesureO2(1) MesureO2(2) MesureO2(3) MesureO2(4).
La donnée répétée est la mesureO2 (4 fois) et les deux facteurs sont Porosité et Gammare. L'unité statistitique est la colonne (il y en a 18).
Avec Systat :
ANOVA
DEPEND MesureO2(1) ...(4)/ REPEAT = 4, NAMES = DATE
CATEGORY Porosité / Gammare
ESTIMATE
Avec R, il faut introduire le facteur date qui caractérise la répétition (Date1...4).
Le modèle (qui a l'air d'être orthogonal) serait : summary(aov(O2~as.factor(DATE)*as.factor(POROSITE)*as.factor(GAMMARE)+Error(COLONNE)). En principe, le résultat donne les variations inter et intra.

[Simon NAVEL]

Merci Eric.

J'avais pensé fut un temps à analyser les données comme cela.
Cependant, je ne tiens pas à mettre l'effet DATE comme fixe, je n'ai pas de tendance au cours du temps qui justifie de le mettre en fixe, de plus ce n'est pas une variable "d'intérêt" pour moi (pris avec des pincettes).
Un modèle comme vous me le suggérer me laisse 47 ddl résiduels pour calculer les F et pvalue des effets fixes.
Or je tiens à n'en conserver que 12.
J'explique pourquoi (de manière sûrement peu conventionnelle), avec les ddl associés:

J'ai d'une part les variabilté "d'interet" ... pris avec des pincettes (non lié à la date, considérée comme inter-colonne):
*intercept* 1ddl
*effet porosite 2ddl
*effet gammare 1ddl
*interaction porosite-gammare 2ddl
*résidus au sein des traitements croisés 12ddl

d'autres part les variabiltés que je souhaite être considérées comme aléatoires (liées à la date).
*effet date 3ddl
*interaction date-porosite 6ddl
*interaction date-gammare 3ddl
*interaction date-porosite-gammare 6ddl
*résidus liés à l'effet date 36ddl

L'idée derrière tout cela, est de "s'affranchir" des variabilités dues à l'effet DATE, puisque je n'ai aucune tendance (j'ai en quelque sorte une population de date). Je souhaite donc faire en sorte que R résonne sur une moyenne des 4 dates confondues, pour chaque colonne. Le résidus "d'intérêt" correspond donc aux variabilités au sein de chacun des croisements de traitement porosite-faune possible (soit la variabilité entre 3 colonnes). Le ddl résiduel, serait donc 6 (nombre de croisements possibles porosite-gammares) x 2 (3colonnes par croisement - l'intercept par croisement), soit 12 ddl, ce que j'obtiens en faisant bêtement l'anova répétée sous statistica. Ce ddl résiduel de 12 est absolument important, il correspond à la variabilité "vraie" inter colonne (au sein des groupes), une fois toute implication de la DATE "ecartée" de l'analyse. C'est à partir de ce résidus à 12 ddl que je souhaite calculer les F et pvalue des effets fixes.

Votre démarche, tout à fait légitime également, raisonne sur la COLONNE plutôt que sur la DATE. C'est tout à fait correct mais ça ne répond pas à ce que j'attends véritablement dans la façon de gérer les résidus.

J'espere que mon explication permettra de mieux comprendre mon problème.
Je vous remercie Eric pour votre aide.

Simon

[Eric Pagot]

J'entrevois quelque peu votre raisonnement, mais de mon point de vue, il est biaisé. A priori, un effet fixe est un effet qui est imposé par l'expérimentateur (par exemple porosité). Donc, pour moi, l'effet date a bien été décidé et que l'on aie une tendance n'a rien à voir dans l'histoire.
Vouloir regrouper toutes les données en une population de dates signifierait que les mesures sont indépendantes. Or, si je comprend bien, pour une colonne donnée, il y a 4 mesures, ce qui signifie que les données ne sont pas indépendantes (liées au sujet). Le test d'ANOVA en mesure répétée permet d'augmenter la puissance. Autrement, la solution est de faire la moyenne des 4 valeurs aux 4 dates et de faire une ANOVA classique avec deux groupes (ce qui diminue le nombre de données) et donnera bien les 12 ddl.

[Simon NAVEL]

Merci Eric pour votre aide.

Je comprend tout à fait vos interrogations quand à ma gestion de l'effet DATE. Cependant l'approche que vous me suggérez ne correspond pas vraiment à ce que je recherche. En effet les 47 ddl résiduels fournis avec l'analyse que vous me suggérez sont très peu interprétables. Pour être peut-être plus clair, j'essaye en fait de refaire l'anova répétée réalisée sous Statistica, qui d'un point de vue des ddl me sort des valeurs cohérentes. Je suis sensible à vos commentaires et tout à fait d'accord avec votre raisonnement ( unité expérimentale = colonne..).
Cependant la gestion des résidus ne correspond pas à ce que j'attends.
Peut-être que l'anova répétée se gère d'une toute autre manière sous R (méthode autre que lme, aov et autres), auquel cas je m'excuse vraiment de n'avoir pas été plus précis plus tôt).
Peut être cette précision apportera quelque chose de nouveau.
Le Stateux que j'ai consulté, approuve "le raisonnement" de Statistica et le nombre de ddl utilisé. Il ne me manque en fait qu'à essayer de retranscrire le même raisonnement sous R.

Dans tous les cas je vous remercie Eric pour votre aide précieuse.


Simon

[Simon NAVEL]

Bien bien bien,

J'ai enfin résolu mon probleme d'anova répétée.
J'ai utilisé Rcmdr et le plugin que Eric a proposé

http://forums.cirad.fr/logiciel-R/viewtopic.php?p=4031

Après m'être familiarisé avec l'interface et après avoir recodé mes données, j'ai pu réalisé l'anova répétée en mode "presse-bouton".
J'obtiens les bons ddl et des valeurs de F et pvalue correspondantes à ce que j'ai pu obtenir sous Statistica.

L'avantage du Rcmdr est de me permettre de disposer du script ayant permit cette analyse. Je vais donc pouvoir étudié un peu la manière employée ici et ainsi mieux comprendre la façon de coder.

Merci pour l'aide apportée.

Simon

[Eric Pagot]

Je suis bien content que mon plugin ait pu servir. Pour trouver les codes il faut reprendre le fichier présent dans le sous-répertoire R. J'ai nommé la fonction "repeatedAnova".

[Simon NAVEL]

Bien recu,
merci beaucoup Eric