dispositif enclos/exclos avec mesures répétées dans le temps

[Christophe Baltzinger]

Bonjour,

Je dois analyser des données avec des mesures répétées régulièrement tous les ans sur un dispositif relativement structuré avec 10 sites contenant chacun 2 répétitions enclos/exclos (pour tester l'effet de l'herbivorie des cerfs sur la végétation : diversité/recouvrement).
Je pensais analyser ces données selon les "repeated measures analysis", je voulais savoir si sur ces modèles d'analyses il était possible de prendre en compte l'appariement enclos1/exclos1 du site i en analysant par exemple comme variable la différence pour cette variable entre enclos et exclos (ce qui risque de donner des valeurs positives et/ou négatives)

merci pour tout éclairement

[Renaud Lancelot]

La question est trop vague pour répondre précisément. Il faudrait un jeu de données en exemple. Cependant, pas de doute a priori: il est possible d'analyser ces données dans un cadre de mesures répétées. Les méthodes sont nombreuses et leur choix dépend des objectifs. Par exemple, le seul intérêt de l'analyse est-il de mettre en évidence l'effet d'herbivorie ou veut-on également mettre en évidence et quantifier des variations liées au site ? Dans le premier cas, des modèles marginaux sont suffisants. Dans le second il faudrait envisager des modèles mixtes. La nature de la variable à expliquer est également importante: si elle est quantitative, cela oriente vers le modèle linéaire, sinon, vers le modèle linéaire généralisé. Dans les deux cas, des extensions permettent d'ajuster des modèles marginaux ou mixtes.

Renaud

[Christophe Baltzinger]

Merci pour ce début de réponse, alors je précise :
le dispositif est composé de 10 sites comprenant chacun 2 couples enclos/exclos : soit au total 20 enclos et 20 exclos.
Les variables à analyser sont quantitatives (richesse totale, recouvrement par strates de hauteur).
L'objectif est effectivement d'étudier l'effet des cervidés au cours du temps, mais il faut contrôler pour l'effet site d'une part et prendre en compte l'appariement des couples enclos/exclos d'autre part (d'où l'idée d'analyser les différences par variable entre enclos et exclos, ce qui a l'intérêt de régler la question de l'appariement et de l'exclusion des cerfs)
les années de mesure sont 1997 à 2001 chaque année et une fois en 2005 ( donc si on considère l'ensemble des campagnes elles ne sont pas également espacées dans le temps, et il faut donc prendre cette spécificité en compte ).

D'avance merci

[Renaud Lancelot]

J'utiliserais un modèle linéaire mixte:

1. je traiterais la variabilité inter-site comme un effet aléatoire,
2. je considérerais les autres effets comme fixes:
- année traitée en facteur
- strate
- enclos (deux modalités)

En supposant les données rassemblées dans un data.frame mydata, le modèle pourrait s'écrire:

Code : Tout sélectionner

library(lme4)
fm <- lmer(recouvrement ~ an + strate + enclos + (1 | site), data = mydata)
summary(mydata)


Avec ce modèle, vous estimez les moyennes de pop à l'aide des effets fixes. L'effet aléatoire permet de rendre compte de la variabilité inter-site autour de cette moyenne: nul en moyenne et distribué selon une loi normale dont la variance est estimée.
Vous pouvez tester l'existence d'interactions entre les effets fixes à l'aide de tests du rapport de vraisemblance, et/ou comparer les modèles à l'aide de critères d'information.

Attention:

- si vous n'avez que 20 mesures / an, il faudra rester avec des modèles très simples: le modèle ci-dessus a déjà bcp de paramètres.

- si la richesse est un comptage avec des valeurs modestes (quelques unités), l'hypothèse de distribution normale des résidus est probablement fausse et il faudra peut-être utilisé un modèle linéaire généralisé mixte.

Renaud

[Christophe Baltzinger]

Merci pour ces conseils Renaud.
Tout n'est pas encore très clair pour moi.
les données répétées dans le temps le sont tous les ans de 1997 à 2001, puis 4 ans après en 2005, j'ai l'impression que dans ce cas il faudrait plutôt prendre la variable an comme continue plutôt que facteur, puisqu'alors l'année 2005 serait considérée comme une année "2002".
Qu'entends tu par effet "enclos" : la protection contre les cerfs ?

Voici la sortie obtenue sur S+, en limitant aux années 1997 à 2001, 4 valeurs manquantes et prise en compte de l'interaction entre t=an et protection !

> fm<-lme(richesse~factor(t)*factor(protection), random=~site, na.action=na.exclude,data=hg0)
Warning messages:
SINGULAR CONVERGENCE. in: ms( ~ - logLik(lmeSt, lmePars), start = list(lmePars = c(coef(lmeSt))), control = list(rel.tolerance = controlvals$msTol, maxiter =
controlvals$msMaxIter, scale = controlvals$ ....
> summary(fm)
Linear mixed-effects model fit by REML
Data: hg0
AIC BIC logLik
1035.308 1080.468 -503.6538

Random effects:
Formula: ~ site | 1
Structure: General positive-definite
StdDev Corr
(Intercept) 0.5908897 (Inter
site 0.2564936 0
Residual 3.1021711

Fixed effects: richesse ~ factor(t) * factor(protection)
Value Std.Error DF t-value p-value
(Intercept) 10.77456 1.245479 186 8.650942 <.0001
factor(t)1 0.12500 0.346833 186 0.360404 0.7190
factor(t)2 0.37500 0.200244 186 1.872713 0.0627
factor(t)3 0.51250 0.141594 186 3.619501 0.0004
factor(t)4 0.57504 0.114522 186 5.021167 <.0001
factor(protection) -0.32667 0.221781 186 -1.472928 0.1425
factor(t)1factor(protection) 0.15000 0.346833 186 0.432484 0.6659
factor(t)2factor(protection) 0.05000 0.200244 186 0.249695 0.8031
factor(t)3factor(protection) 0.02500 0.141594 186 0.176561 0.8600
factor(t)4factor(protection) -0.12667 0.114449 186 -1.106751 0.2698

Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.000475 -0.6361214 -0.1107242 0.473544 4.441023

Number of Observations: 196
Number of Groups: 1

Merci pour tout éclairement

[Renaud Lancelot]

Non, je ne crois pas qu'il faut prendre l'année comme une variable quantitative sauf si une exploration graphique montre qu'il y a bien une tendance linéaire (croissance ou décroissance continue dans le temps).

Le modèle ne converge pas car la formulation des effets aléatoires est probablement inappropriée: aboutit à une structure aléatoire comportant n * (n + 1) / 2 paramètres, où n est le nb de modalités de site. Je ne pense pas que ce soit cela que vous voulez, mais plutôt un seul effet aléatoire lié au site (fluctuation "aléatoire" de la moyenne de richesse selon le site). Avec lme, il faut écrire:

fm1 <- lme(richesse ~ factor(t) * factor(protection), random = ~ 1 | site, na.action = na.exclude, data = hg0)

Si vous avez encore des pbs de convergence, essayer le modèle additif

fm2 <- lme(richesse ~ factor(t) * factor(protection), random = ~ 1 | site, na.action = na.exclude, data = hg0)

Les modèles équivalents avec le package lme4 et la fct lmer sont:

library(lme4)
fm3 <- lmer(richesse ~ factor(t) * factor(protection) + (1 | site), data = hg0)
fm4 <- lmer(richesse ~ factor(t) + factor(protection) + (1 | site), data = hg0)

Il vous faudrait lire un bon bouquin sur les modèles mixtes. Si vous utilisez lme, voir

Pinheiro, J.C., Bates, D.M., 2000. Mixed-effect models in S and S-Plus. Springer, New York (USA), 598 p.

Renaud

[Noémie Stroh]

Bonjour,
je travaille avec Christophe Baltzinger, sur les mêmes données ; on voudrait donc formaliser un modèle pour ananlyser nos données de recouvrement et de richesse spécifique. Le dispositif expérrimental est le suivant :

10 sites d'étude sur lesquels 2 dispositifs d'enclos / exclos sont installés ; on a donc 2 réplications des mesures sur chaque site. Les mesures sont faites sur 6 années : 1997 à 2001 puis en 2005.

Donc on voudrait utilisé un modèlelinéaire mixte en prenant :
- l'effet traitement "enclos/exclos" (A) et l'effet "temps" (B) comme effets fixes
- l'effet "site" comme effet aléatoire (b)
Il faudrait également prendre en compte une corrélation temporelle dans nos données car les mesures sont faites chaque année sur les mêmes 10 sites.

Donc le modèle serait du type :

Yijkt = variable (recouvrement ou richesse) au site i , pour le traitement j , à la réplication k , à l'année t
= u + Aj + Bt + bi + Eijkt

i = 1 à 10 (sites)
j = 1, 2 (enclos ou esclos)
k = 1, 2 (dispositif 1 ou 2 du site)
t = 1,...5, 8 (1997à 2001, 2005)

Donc, plusieurs questions :
-d'abord est-ce cohérent ? et est-ce que les mesures répliquées 2 fois sur chaque sites sont biens prises en compte?
- si oui, comment le formaliser sous R et surtout comment formaliser la corrélation temporelle (dans le livre de Pinheiro sur les modèles mixtes sur S et S+, il est question de l'argument "correlation" pour spécifié les corrélations dans lme et des objets "corStruct", je comprend le principe, mais je vois mal comment cela se goupille..).

Somme nous sur la bonne voie?

Merci d'avance pour tout aide

Noémie

[Renaud Lancelot]

Il y a sur ce forum un exemple traité (code, données et publi) de construction d'un modèle mixte sur des données de croissance. Le pb se rapproche bcp du vôtre. Voir http://forums.cirad.fr/logiciel-R/viewtopic.php?t=245

Dans le cas de vos données, je ne vois pas bien si vous avez "seulement" 2 observations par an et par site, ou si les observations sont répétées chaque année pour chaque catégorie enclos / exclos. Si vous avez bcp de répétitions par site, vous pouvez vous permettre de coder l'année en facteur et de mettre cette variable dans les effets fixes. Le modèle peut être ajusté avec bcp de fcts différentes dans plusieurs packages R, par exemple lme dans le package nlme:

Code : Tout sélectionner

library(nlme)
fm1 <- lme(recouvrement ~ A + B + T, random = ~ 1 | S, data = mydata)

où mydata est le data.frame contenant les données et s est la variable décrivant le site.

* Attention au codage de T (doit être un facteur). En passant, ce serait une très mauvaise idée d'utiliser les noms t et s dans le script car ce sont des noms de fonction R ==> risque de conflit.
* dans ce cas, l'autocorrélation des observations est prise en compte par l'effet aléatoire et il n'y a pas lieu de modéliser l'autocorrélation des résidus, sauf à rechercher une structure spatiale.

Sinon (2 observations / site / an), vous n'aurez pas assez de données pour estimer correctement ce modèle et il faudra une autre stratégie. Je vois deux possibilités:

1. Mettre l'année en effet aléatoire, i.e., vous aurez deux effets aléatoires "croisés" (crossed random effects). La solution la plus simple sous R est d'utiliser le package lme4, le modèle pouvant s'écrire:

Code : Tout sélectionner

library(lme4)
fm2 <- lmer(recouvrement ~ A + B + (1 | t) + (1 | s), data = mydata)


Dans ce cas, pas besoin de spécifier une structure d'autocorrélation des résidus.

2. Si l'exploration graphique révèle une structure de corrélation simple dans les données, vous pouvez ignorer l'année dans les effets fixes et aléatoires, et représenter l'autocorrélation des résidus avec un modèle explicite. Je n'ai pas le temps de rentrer dans les détails: consulter la référence indiquée avec l'analyse des données de croissance. Il vous faudra utiliser le package nlme:

Code : Tout sélectionner

library(nlme)
fm3 <- lme(recouvrement ~ A + B,
                 random = ~ 1 | S,
                 correlation = corAR1(form = ~ 1 | Time),
                 data = mydata)

pour modéliser l'autocorrélation des observations successives par une modèle auto-régressif d'ordre 1. Attention, il faudrait probablement coder l'année sous forme de valeurs entières commençant à 1 (1997 = 1, 1998 = 2,...).

Il faut du temps pour comprendre ces modèles et l'ensemble des concepts et méthodes qui sont mises en oeuvre dans ce cas. Le bouquin de Pinheiro et Bates est très bien fait mais il vous faudra y consacrer qques semaines de lecture, exercices, etc.

Attention, les packages nlme et lme4 ne doivent pas être chargés au cours de la même session (conflits) ==> démarrer deux sessions parallèles ou successives pour ajuster les modèles.

Bon courage,

Renaud