Anova en mesure répétée et données manquantes

[Caroline Rondel]

Bonjour,

J'ai réalisée une expérience avec 6 traitements différents (contrôle compris). Chaque traitement est répliquée 3 fois. Les mesures des différents paramètres ont été prises à des pas de temps variables durant l'expérience.
Malheureusement, il nous manque les valeurs de certains paramètres à différentes dates pour des traitements.

Je veux savoir s'il existe des différences entre les traitements et donc réaliser un ANOVA en mesures répétées.

Quels fonctions vaut-il mieux utiliser dans mon cas?
La fonction lm ou lme?

Et comment ces différentes fonctions vont-elles se comporter par rapport aux données manquantes?

Egalement, pouvez vous me dire si la correction de seuil de probabilité en fonction du nombre de paramètres testés peut se faire automatiquement sous R?

Caroline Rondel
Université Montpellier II, Laboratoire Ecosystèmes Lagunaires,
IRD, UR 167 CyRoCo
Place Eugène Bataillon
Case 093; F- 34095 Montpellier cedex 05

[Renaud Lancelot]

Bonjour Caroline,

Pour éviter les erreurs d'interprétation (de ma part), l'idéal serait de nous montrer le jeu de données, ou une portion suffisante pour que nous comprenions bien ce dont vous parlez.

Si les répétitions concernent des unités différentes, vous êtes dans un cas d'analyse de variance "classique" avec un plan d'expérience déséquilibré. Vous pouvez utilisez la fonction lm. En utilisant des contrastes appropriés, vous pourrez comparer les traitements au traitement de référence à l'aide du test de t fourni en sortie standard de la méthode summary pour lm. Vous pouvez également utiliser le test de Wald pour tester des combinaisons linéaires des paramètres. Voir la fonction linear.hypothesis du package car, ou la fonction wald.test du package aod, qui peuvent prendre en compte le nbre de ddl des résidus.

Si les répétitions sont faites sur les mêmes unités, il faut tenir compte de l'éventuelle source de variation due à des réponses se ressemblant pour un même individu. Différentes méthodes sont possibles:
* si vous ne vous intéressez qu'aux moyennes (comparaison des traitements), vous pouvez utiliser des modèles marginaux tels que les moindres carrés généralisés (package nlme, fonction gls) ou les modèles GEE (fonction geese du package geepack),
* si en plus des moyennes, vous vous intéressez aux effets des individus, vous pouvez utiliser des modèles à effets mixtes comme ceux ajustés par la fonction lme du package nlme, ou lmer des packages Matrix et lme4.
Ces modèles permettent d'analyser des plans d'expérience déséquilibrés.

Une question importante est de savoir si les données manquantes se produisent de manière indépendante du processus étudié:
* si ce n'est pas le cas, les résultats sont biaisés,
* si c'est la cas, il faut quand même faire attention à l'interprétation des paramètres estimés (conditionnelle aux données manquantes).

Pour une présentation générale des méthodes adaptées aux données répétées, voir par exemple:

Diggle PJ, Liang K-Y and Zeger SL, 1994. Analysis of longitudinal data, 1st ed. Clarendon Press, Oxford.

Il y a une édition plus récente (2002 ?). Cet ouvrage discute également du pb des données manquantes et présente des méthodes de diagnostic de la dépendance entre données manquante et réponses.

Les modèles linéaires mixtes et les modèles gls sont présentés simplement et avec beaucoup d'exemples dans le bouquin accompagnant le package nlme:

Pinheiro JC and Bates DM, 2000. Mixed-effect models in S and S-Plus, ed. Springer, New York.

Bien cordialement,

Renaud

[Caroline Rondel]

Bonjour,

Merci pour ces explications.
Je vous expose ci-après le plan d'expérience et un exemple du jeu de donnée que je possède.

D'après ce que vous me dites il me semble que je devrais plutôt utiliser la fonction lm à moins que j'ai mal compris mal la notion d'unitée?
Pourriez-vous également me dire s'il s'il existe une fonction corrigeant le seuil de probabilité vu que je vais tester un très grand nombres de paramètres (environ 50)?


Il s'agit d'une expérience en mésocosme.
Mon plan d'expérience au départ est un croisement entre l'effet de jeune poisson et 2 autres stades de développement du même poisson:

ZP (+/-) xFish (0, GL, OM)

Des analyses grâce à super Anova avaient montré très peu d'effet d'interaction entre les groupes et me font penser qu'il serait plus simple de choisir une anova à 1 facteur plutôt que de tenir compte de l'interaction. A tort peut être?

Les données manquantes sont dû à la perte de certains mésocosmes en fin d'expérience et donc commun à l'ensemble des paramètres.

Voici un exemple de tableau de données, le paramètre choisi est le rapport N/P (X correpondant aux valeurs manquantes):

Traitement Mésocosme t0 t3 t6 t9 t12 t15 t18 t24
GL 2 13,61 9,69 8,59 10,31 10,22 10,39 9,61 9,35
GL 10 15,61 8,34 9,29 10,32 13,86 10,05 13,56 11,17
GL 18 14,08 8,95 8,92 9,44 11,29 10,18 13,79 12,28
GL + ZP 7 13,37 10,38 9,14 11,21 10,64 13,42 X X
GL + ZP 11 11,48 8,29 8,89 10,15 13,92 12,29 15,56 9,49
GL + ZP 15 15,16 9,26 9,72 9,83 11,53 11,93 13,56 11,83
OM 4 15,15 9,49 8,65 10,12 11,26 11,89 10,34 9,56
OM 9 16,77 9,44 8,46 10,29 14,15 10,15 11,96 8,82
OM 13 14,48 8,84 9,26 10,05 12,02 13,88 13,43 10,94
ZP 3 12,58 8,24 8,55 11,38 11,34 12,89 11,77 8,32
ZP 6 14,86 7,93 8,77 11,88 11,71 10,86 13,82 11,29
ZP 14 13,42 10,51 10,63 10,34 11,80 12,49 9,05 11,69
ZP + OM 1 10,87 10,61 7,01 13,27 12,96 13,46 7,36 X
ZP + OM 8 14,92 8,60 8,65 10,45 9,87 12,55 13,62 9,50
ZP + OM 17 14,72 9,49 9,24 9,06 11,92 10,54 12,53 8,68
Sans 5 13,77 8,86 8,85 12,41 12,18 13,41 11,54 9,60
Sans 12 15,28 9,43 9,35 11,55 14,34 12,23 13,32 14,70
Sans 16 16,19 8,52 9,66 10,17 11,49 12,13 12,61 10,47

[Renaud Lancelot]

Bonjour Caroline,

J'ai légèrement édité votre jeu de données pour qu'il soit exploitable sous R et je l'ai placé à ftp://ftp.cirad.fr/pub/group-r/groupe-r/Data/rondel.txt. De cette manière, les lecteurs intéressés pourront intervenir et mettre en application leurs suggestions.

Je confirme qu'il s'agit bien d'une anova classique. Cependant, le facteur temps semble jouer un rôle important. Je ne suis pas sûr que vous puissiez le laisser de côté:

Code : Tout sélectionner

> # récupère les données
> rondel <- read.table("d:/analyses/travail/data/rondel.txt", header = TRUE)
>
> # fichier
> rondel
   Traitement Mesocosme    t0    t3    t6    t9   t12   t15   t18   t24
1          GL         2 13.61  9.69  8.59 10.31 10.22 10.39  9.61  9.35
2          GL        10 15.61  8.34  9.29 10.32 13.86 10.05 13.56 11.17
3          GL        18 14.08  8.95  8.92  9.44 11.29 10.18 13.79 12.28
4       GL.ZP         7 13.37 10.38  9.14 11.21 10.64 13.42    NA    NA
5       GL.ZP        11 11.48  8.29  8.89 10.15 13.92 12.29 15.56  9.49
6       GL.ZP        15 15.16  9.26  9.72  9.83 11.53 11.93 13.56 11.83
7          OM         4 15.15  9.49  8.65 10.12 11.26 11.89 10.34  9.56
8          OM         9 16.77  9.44  8.46 10.29 14.15 10.15 11.96  8.82
9          OM        13 14.48  8.84  9.26 10.05 12.02 13.88 13.43 10.94
10         ZP         3 12.58  8.24  8.55 11.38 11.34 12.89 11.77  8.32
11         ZP         6 14.86  7.93  8.77 11.88 11.71 10.86 13.82 11.29
12         ZP        14 13.42 10.51 10.63 10.34 11.80 12.49  9.05 11.69
13      ZP.OM         1 10.87 10.61  7.01 13.27 12.96 13.46  7.36    NA
14      ZP.OM         8 14.92  8.60  8.65 10.45  9.87 12.55 13.62  9.50
15      ZP.OM        17 14.72  9.49  9.24  9.06 11.92 10.54 12.53  8.68
16       Sans         5 13.77  8.86  8.85 12.41 12.18 13.41 11.54  9.60
17       Sans        12 15.28  9.43  9.35 11.55 14.34 12.23 13.32 14.70
18       Sans        16 16.19  8.52  9.66 10.17 11.49 12.13 12.61 10.47
>
> # mise en forme pour l'analyse
> rondel2 <- reshape(data = rondel,
+                    varying = list(names(rondel)[3:10]),
+                    idvar = c("Traitement", "Mesocosme"),
+                    direction = "long", v.names = "NP")
> rownames(rondel2) <- seq(nrow(rondel2))
>
> # début du fichier
> head(rondel2)
  Traitement Mesocosme time    NP
1         GL         2    1 13.61
2         GL        10    1 15.61
3         GL        18    1 14.08
4      GL.ZP         7    1 13.37
5      GL.ZP        11    1 11.48
6      GL.ZP        15    1 15.16
>
> # récupère les vrais temps
> tps <- as.numeric(substring(names(rondel)[3:10], 2, nchar(names(rondel)[3:10])))
> rondel2$Time <- rondel2$time
> for(i in seq(length(tps)))
+   rondel2$Time[rondel2$time == i] <- tps[i]
>
> # élimine les lignes avec données manquantes
> rondel3 <- na.omit(rondel2)
> rondel4 <- rondel3[order(rondel3$Mesocosme, rondel3$Traitement, rondel3$Time), ]
> head(rondel4)
    Traitement Mesocosme time    NP Time
13       ZP.OM         1    1 10.87    0
31       ZP.OM         1    2 10.61    3
49       ZP.OM         1    3  7.01    6
67       ZP.OM         1    4 13.27    9
85       ZP.OM         1    5 12.96   12
103      ZP.OM         1    6 13.46   15
>
> # exploration graphique
> library(lattice)
> trellis.device(color = FALSE)
> xyplot(NP ~ Time | Traitement, data = rondel4,
+   panel = function(x, y){
+     panel.grid()
+     panel.xyplot(x, y)
+     panel.loess(x, y)


Ce dernier graphique montre un effet non linéaire du temps. Sans rentrer dans des modèles non linéaires pouvant rendre compte de l'explication biologique sous-jacente, on peut simplement se débarraser de cette variation en utilisant des splines (par exemple).

NB: vrai ici car il ne semble pas y avoir d'interaction temps * traitement: ttes les courbes ont à peu près le même profil.

Code : Tout sélectionner

> library(splines)
>
> # avant l'analyse, on recode le facteur traitement:
> levels(rondel4$Traitement) <- c("Sans", "GL", "OM", "Zp", "GL.ZP", "ZP.OM")
>
> fm1 <- lm(NP ~ bs(Time, df = 5) + Traitement, data = rondel4)
> par(mfrow = c(2, 2))
> plot(fm1)
> summary(fm1)

Call:
lm(formula = NP ~ bs(Time, df = 5) + Traitement, data = rondel4)

Residuals:
    Min      1Q  Median      3Q     Max
-4.3292 -0.8202 -0.1068  0.8457  3.6836

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)   
(Intercept)        13.9996     0.4088  34.247  < 2e-16 ***
bs(Time, df = 5)1  -6.4741     0.7336  -8.826 6.30e-15 ***
bs(Time, df = 5)2  -4.8350     0.8207  -5.891 3.09e-08 ***
bs(Time, df = 5)3  -0.4866     1.0046  -0.484   0.6289   
bs(Time, df = 5)4  -3.0678     1.0890  -2.817   0.0056 **
bs(Time, df = 5)5  -3.7815     0.4676  -8.088 3.72e-13 ***
TraitementGL        0.4627     0.4018   1.152   0.2516   
TraitementOM        0.2708     0.3925   0.690   0.4915   
TraitementZp        0.7983     0.3925   2.034   0.0440 * 
TraitementGL.ZP     0.1342     0.3925   0.342   0.7331   
TraitementZP.OM    -0.1218     0.3970  -0.307   0.7595   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.36 on 130 degrees of freedom
Multiple R-Squared: 0.6179,     Adjusted R-squared: 0.5885
F-statistic: 21.02 on 10 and 130 DF,  p-value: < 2.2e-16


Dans ce cas, seul le ttt "Zp" semble avoir un effet significatif sur NP.

Pour ma part, j'aurais conservé les deux variables Zp et Fish sous leur forme initiale, quitte à éliminer les interactions, voire les effets principaux s'ils s'avéraient sans influence sur la réponse.

Je ne comprends pas bien ce que vous voulez dire par "corriger le seuil de probabilité" ? Si vous avez un pb de sélection de modèles induisant des comparaisons multiples, vous pouvez utiliser le critère d'information d'Akaike. Il y a une procédure de sélection pas-à-pas disponible dans le package MASS: fonction stepAIC.

Cordialement,

Renaud

[Caroline Rondel]

Bonjour,

Vu que mes connaissances statistiques sont encores balbutiantes, je me permet de vous poser de nouvelles questions peut être naïve.

Pourquoi éliminez-vous les lignes contenant les données manquantes? En simulant des valeurs (par exemple en utilisant la moyenne pour le groupe) le test n'aurait-il pas plus de puissance?

En quoi exactement la variation non linéaire pose t'elle problème pour l'analyse? Il me semblait qu'en traitant le temps comme un effet aléatoire cela prenait en compte la variation intra-sujet?

Je comprends bien l'idée de ne pas mettre de l'effet d'interaction des facteurs ZP et Fish dans le modèle structurel du design ZPxFish puisque dans mon cas celui-ci s'avère négligeable. Cependant, quels informations peut m'apporter un modèle structurel où ne figure ni l'interaction ni les effets principaux? Il ne resterait alors que l'effet temps?

NB:
Ma question concernant la correction de seuil de probalité vient dû fait que lorsqu'on réalise des comparaisons de moyenne plusieurs fois dans une même expérience, on effectue, par exemple, une correction de bonferroni du seuil Alpha.
Dans mon cas il ne s'agit pas de multiplier les comparaisons de moyennes (pas encore) mais de répéter de nombreuses fois un même type d'ANOVA sur des paramètres différents d'une même expérience. Je me demandais alors si le risque d'erreur Beta n'augmentais pas et s'il existait une fonction sous R qui pouvait corriger l'erreur que je fais.

En vous remerciant de toutes les informations que vous m'avez déjà apporté.

Cordialement,[/code]

[Renaud Lancelot]

Bonjour,

Avant de répondre à vos questions, je reviens sur ce que je disais dans mon précédent message. Vos observations sont répétées sur les mêmes unités (les mésocosmes) et en toute rigueur, il faudrait utiliser des méthodes appropriés pour analyser les données. Toutefois, le jeu que vous présentez permet de faire comme si les observations étaient indépendantes. Exemple avec un modèle linéaire mixte:

Code : Tout sélectionner

> library(nlme)
> # modèle linéaire mixte
> fm2 <- lme(NP ~ bs(Time, df = 5) + Traitement, random = ~ 1 | Mesocosme, data = rondel4)
>
> # variance de l'effet aléatoire
> VarCorr(fm2)
Mesocosme = pdLogChol(1) 
            Variance     StdDev     
(Intercept) 2.263515e-08 0.0001504498
Residual    1.849050e+00 1.3597976704


La variance de l'effet aléatoire est très faible, et très inférieure à la variance résiduelle.

Réponses à vos questions:

1) Le temps représente une source de variation importante dans vos expériences et ce serait une erreur de ne pas en tenir compte, surtout si (comme ici) le plan d'expérience est déséquilibré. Il y a peu de données manquantes dans votre exemple, mais si elles étaient plus nombreuses et réparties différemment dans les différents mésocosmes, la comparaison des réponses entre mésocosmes serait biaisée sans prise en compte du temps. En d'autres termes, si vous incluez le temps dans le modèle, la comparaison statistiques des différents traitements se fait pour un temps fixé, le même pour tous les traitements.

2) Données manquantes. Comme indiqué dans ma première réponse, le traitement des données manquantes est un pb difficile. Si vous avez "beaucoup" de données manquantes, la première des choses est de s'assurer de leur statut vis-à-vis de la réponse. Je vous renvoie à la référence citée (Diggle et al.) pour une discussion plus précise. Vous pouvez effectivement remplacer les données manquantes par des valeurs arbitraires. En choisissant la moyenne, vous diminuez artificiellement la variance et risquez de biaiser les tests. Il est préférable d'utilise des techniques d'imputation multiples. Voir par exemple:

Schafer J-L, 1997. Analysis of incomplete multivariate data, ed. Chapman & Hall / CRC, London, 448 p.

Il y a des packages R qui font cela: voir par exemple mix.

3) Je préfèrerais laisser les variables en l'état, estimer les différents modèles plausibles et les comparer avec un test du rapport des vraisemblances ou s'ils sont nombreux, avec le critère d'information d'Akaike. Si une ou plusieurs variables vous semblent importantes, vous pouvez les imposez dans les modèles même si elles ne sont pas significatives. Vous aurez ainsi des estimation conditionnelles à chacune des combinaisons de ces variables.

4) Seuil de probabilité. Deux possibilités:

* faire une analyse de données multivariées: d'abord ACP ou ACM pour étudier les relations entre les réponses, et ensuite une ACPVI ou AFCVI pour étudier l'effet des variables explicatives sur le tableau des variables à expliquer. C'est un autre monde. Les méthodes sont disponibles dans le package ade4, accompagné d'une documentation abondante. Voir http://pbil.univ-lyon1.fr/ADE-4/ADE-4F.html.

* faire une analyse de variance multivariée. Voir les fonctions manova et summary.manova:

Code : Tout sélectionner

> example(summary.manova)

smmry.> tear <- c(6.5, 6.2, 5.8, 6.5, 6.5, 6.9, 7.2, 6.9,
    6.1, 6.3, 6.7, 6.6, 7.2, 7.1, 6.8, 7.1, 7, 7.2, 7.5, 7.6)

smmry.> gloss <- c(9.5, 9.9, 9.6, 9.6, 9.2, 9.1, 10, 9.9,
    9.5, 9.4, 9.1, 9.3, 8.3, 8.4, 8.5, 9.2, 8.8, 9.7, 10.1, 9.2)

smmry.> opacity <- c(4.4, 6.4, 3, 4.1, 0.8, 5.7, 2, 3.9, 1.9,
    5.7, 2.8, 4.1, 3.8, 1.6, 3.4, 8.4, 5.2, 6.9, 2.7, 1.9)

smmry.> Y <- cbind(tear, gloss, opacity)

smmry.> rate <- factor(gl(2, 10), labels = c("Low", "High"))

smmry.> additive <- factor(gl(2, 5, len = 20), labels = c("Low",
    "High"))

smmry.> fit <- manova(Y ~ rate * additive)

smmry.> summary.aov(fit)
 Response tear :
              Df  Sum Sq Mean Sq F value   Pr(>F)   
rate           1 1.74050 1.74050 15.7868 0.001092 **
additive       1 0.76050 0.76050  6.8980 0.018330 *
rate:additive  1 0.00050 0.00050  0.0045 0.947143   
Residuals     16 1.76400 0.11025                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Response gloss :
              Df  Sum Sq Mean Sq F value  Pr(>F) 
rate           1 1.30050 1.30050  7.9178 0.01248 *
additive       1 0.61250 0.61250  3.7291 0.07139 .
rate:additive  1 0.54450 0.54450  3.3151 0.08740 .
Residuals     16 2.62800 0.16425                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Response opacity :
              Df Sum Sq Mean Sq F value Pr(>F)
rate           1  0.421   0.421  0.1036 0.7517
additive       1  4.901   4.901  1.2077 0.2881
rate:additive  1  3.960   3.960  0.9760 0.3379
Residuals     16 64.924   4.058               


smmry.> summary(fit, test = "Wilks")
              Df  Wilks approx F num Df den Df   Pr(>F)   
rate           1 0.3819   7.5543      3     14 0.003034 **
additive       1 0.5230   4.2556      3     14 0.024745 *
rate:additive  1 0.7771   1.3385      3     14 0.301782   
Residuals     16                                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Cordialement,

Renaud

[Renaud Lancelot]

Bonjour,

Je vous remercie beaucoup pour le temps que vous m'avez accordé et vos réponses.

J'ai compris qu'en prenant le temps comme un effet aléatoire, je ne pouvais plus alors tester son effet sur le paramètre. Et comme vous l'avez noté il y'a une importance biologique à cette disparité du temps puisqu'il y'a eu une période avec ajout de nutriment et une période sans nutriment. Il serait donc dommage de ne pas le prendre en compte.

Ma dernière question est : comment faire si je dois tester des paramètres où il semble qu'il y'est une interaction entre le temps et les traitements. Je devrais donc à ce moment le traiter en tant qu'effet aléatoire.

Malheureusement, si j'utilise la fonction lme, je ne peux plus tester l'effet du temps.
J'avais envisagé de faire une Anova de type III avec le package "car". Mais, je ne retrouve pas les mêmes résultats que vous :

Code : Tout sélectionner

fm3<-lm(NP~Time*Traitement,data=rondel4)
library(car)
Anova(fm3,type="III")

Anova Table (Type III tests)

Response: NP
                Sum Sq  Df  F value Pr(>F)   
(Intercept)     940.99   1 200.2358 <2e-16 ***
Traitement        5.60   5   0.2384 0.9449   
Time              0.08   1   0.0172 0.8957   
Traitement:Time   9.92   5   0.4220 0.8327   
Residuals       606.22 129                   
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1


Et pour la fonction lm

Code : Tout sélectionner

anova(fm1)

Analysis of Variance Table

Response: NP
                  Df Sum Sq Mean Sq F value Pr(>F)   
bs(Time, df = 5)   5 375.54   75.11 40.6197 <2e-16 ***
Traitement         5  13.21    2.64  1.4291  0.218   
Residuals        130 240.38    1.85                   
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1


Par exemple, dans le jeu de donnée suivant issu de la même expérience, je ne pense pas pouvoir utiliser le temps comme un facteur fixe:

Code : Tout sélectionner

date   mesocosme   traitement   PL   Fish   PETITHERB   GRANDHERB   CARNIVORES
0   M5   sans   0   0   148,409   276,310   33,167
1   M5   sans   0   0   153,586   213,725   106,418
3   M5   sans   0   0   114,077   290,413   14,842
6   M5   sans   0   0   1876,367   496,445   406,202
9   M5   sans   0   0   444,119   374,859   255,015
12   M5   sans   0   0   74,246   936,131   289,371
15   M5   sans   0   0   31,446   1028,866   260,572
18   M5   sans   0   0   79,205   872,444   227,063
24   M5   sans   0   0   65,040   168,601   47,686
0   M12   sans   0   0   122,334   366,861   48,796
1   M12   sans   0   0   380,844   892,988   100,015
3   M12   sans   0   0   199,398   1086,968   10,149
6   M12   sans   0   0   143,873   513,053   609,254
9   M12   sans   0   0   482,000   692,561   87,331
12   M12   sans   0   0   81,003   824,556   162,766
15   M12   sans   0   0   37,320   351,674   186,040
18   M12   sans   0   0   67,663   927,863   939,887
24   M12   sans   0   0   15,615   508,334   511,207
0   M16   sans   0   0   125,781   374,700   181,691
1   M16   sans   0   0   222,605   404,226   73,937
3   M16   sans   0   0   206,599   744,915   140,930
6   M16   sans   0   0   198,991   466,546   274,761
9   M16   sans   0   0   599,589   453,345   167,144
12   M16   sans   0   0   70,659   655,184   302,109
15   M16   sans   0   0   90,609   460,652   151,127
18   M16   sans   0   0   60,505   604,614   152,056
24   M16   sans   0   0   20,860   651,255   456,134
0   M3   PL   1   0   100,918   191,279   28,661
1   M3   PL   1   0   196,624   363,385   48,167
3   M3   PL   1   0   65,768   161,607   16,591
6   M3   PL   1   0   359,668   141,929   284,688
9   M3   PL   1   0   574,427   1926,323   230,478
12   M3   PL   1   0   101,038   1189,067   231,276
15   M3   PL   1   0   84,853   826,451   177,339
18   M3   PL   1   0   36,699   1039,361   393,822
24   M3   PL   1   0   22,196   96,654   1507,883
0   M6   PL   1   0   131,711   255,623   30,262
1   M6   PL   1   0   227,646   10,755   17,684
3   M6   PL   1   0   23,355   55,565   13,336
6   M6   PL   1   0   508,744   22,136   352,743
9   M6   PL   1   0   814,420   171,144   176,626
12   M6   PL   1   0   177,797   35,775   123,084
15   M6   PL   1   0   232,469   200,265   186,190
18   M6   PL   1   0   157,781   448,633   265,461
24   M6   PL   1   0   20,754   291,963   181,048
0   M14   PL   1   0   157,254   209,459   35,338
1   M14   PL   1   0   246,851   397,175   54,762
3   M14   PL   1   0   124,147   542,973   100,420
6   M14   PL   1   0   323,729   630,835   280,376
9   M14   PL   1   0   197,966   535,368   223,113
12   M14   PL   1   0   75,401   982,440   333,628
15   M14   PL   1   0   84,771   559,138   426,787
18   M14   PL   1   0   56,654   453,909   358,582
24   M14   PL   1   0   36,496   295,767   548,778
0   M2   AL   0   AL   75,299   230,483   25,633
1   M2   AL   0   AL   177,677   273,294   93,258
3   M2   AL   0   AL   226,794   126,936   19,070
6   M2   AL   0   AL   429,533   29,167   153,245
9   M2   AL   0   AL   386,026   403,713   17,835
12   M2   AL   0   AL   37,168   56,412   12,770
15   M2   AL   0   AL   35,610   107,244   32,253
18   M2   AL   0   AL   49,526   171,298   37,108
24   M2   AL   0   AL   108,557   240,782   42,071
0   M10   AL   0   AL   129,051   230,227   28,088
1   M10   AL   0   AL   104,351   249,185   23,539
3   M10   AL   0   AL   116,162   336,998   160,687
6   M10   AL   0   AL   124,079   277,406   231,953
9   M10   AL   0   AL   377,216   0,696   5,133
12   M10   AL   0   AL   82,326   0,906   6,192
15   M10   AL   0   AL   139,130   0,988   7,955
18   M10   AL   0   AL   151,351   6,915   3,204
24   M10   AL   0   AL   21,096   56,837   5,049
0   M18   AL   0   AL   146,319   351,636   26,901
1   M18   AL   0   AL   107,477   452,827   49,428
3   M18   AL   0   AL   176,059   791,504   137,144
6   M18   AL   0   AL   167,233   586,895   69,718
9   M18   AL   0   AL   32,050   2,283   3,071
12   M18   AL   0   AL   46,145   7,424   26,612
15   M18   AL   0   AL   85,574   56,796   13,467
18   M18   AL   0   AL   73,891   5,562   6,484
24   M18   AL   0   AL   146,717   18,106   8,333
0   M7   ALPL   1   AL   158,828   226,809   27,117
1   M7   ALPL   1   AL   210,239   374,825   37,588
3   M7   ALPL   1   AL   106,705   75,014   17,875
6   M7   ALPL   1   AL   226,634   272,837   227,495
9   M7   ALPL   1   AL   488,365   248,581   50,307
12   M7   ALPL   1   AL   135,209   10,040   36,372
15   M7   ALPL   1   AL   226,287   29,999   54,115
18   M7   ALPL   1   AL   NA   NA   NA
24   M7   ALPL   1   AL   NA   NA   NA
0   M11   ALPL   1   AL   83,295   234,830   25,969
1   M11   ALPL   1   AL   286,855   606,755   122,982
3   M11   ALPL   1   AL   89,287   537,672   112,250
6   M11   ALPL   1   AL   103,065   42,567   214,532
9   M11   ALPL   1   AL   91,670   4,048   7,632
12   M11   ALPL   1   AL   103,252   3,739   15,128
15   M11   ALPL   1   AL   102,481   4,874   11,326
18   M11   ALPL   1   AL   123,381   7,934   6,224
24   M11   ALPL   1   AL   4,321   69,506   11,698
0   M15   ALPL   1   AL   103,260   256,096   33,267
1   M15   ALPL   1   AL   163,747   517,678   62,070
3   M15   ALPL   1   AL   118,003   621,330   139,005
6   M15   ALPL   1   AL   46,728   395,608   95,893
9   M15   ALPL   1   AL   39,530   18,052   12,226
12   M15   ALPL   1   AL   39,668   15,753   7,572
15   M15   ALPL   1   AL   137,215   308,100   54,476
18   M15   ALPL   1   AL   133,505   19,565   13,874
24   M15   ALPL   1   AL   46,715   50,115   104,120
0   M4   OM   0   OM   89,916   285,824   30,236
1   M4   OM   0   OM   107,875   462,271   80,417
3   M4   OM   0   OM   31,522   140,439   18,182
6   M4   OM   0   OM   145,951   1121,387   388,205
9   M4   OM   0   OM   176,016   412,593   521,575
12   M4   OM   0   OM   98,640   314,691   406,856
15   M4   OM   0   OM   104,945   197,609   205,408
18   M4   OM   0   OM   30,953   408,122   1023,753
24   M4   OM   0   OM   34,301   179,438   231,997
0   M9   OM   0   OM   NA   NA   NA
1   M9   OM   0   OM   91,583   309,544   41,160
3   M9   OM   0   OM   95,833   550,832   10,223
6   M9   OM   0   OM   276,576   467,613   105,774
9   M9   OM   0   OM   342,041   109,657   29,192
12   M9   OM   0   OM   339,872   23,687   133,417
15   M9   OM   0   OM   446,049   22,989   90,680
18   M9   OM   0   OM   388,985   40,306   144,185
24   M9   OM   0   OM   56,980   71,036   378,069
0   M13   OM   0   OM   120,285   213,760   70,870
1   M13   OM   0   OM   128,423   315,255   61,798
3   M13   OM   0   OM   137,810   736,002   209,851
6   M13   OM   0   OM   133,540   603,777   292,508
9   M13   OM   0   OM   166,406   4468,520   329,410
12   M13   OM   0   OM   150,953   267,074   152,249
15   M13   OM   0   OM   157,113   166,164   123,806
18   M13   OM   0   OM   163,314   470,811   554,357
24   M13   OM   0   OM   7,275   140,815   194,198
0   M1   PLOM   1   OM   139,696   393,474   27,742
1   M1   PLOM   1   OM   160,609   262,255   45,150
3   M1   PLOM   1   OM   19,613   302,152   18,139
6   M1   PLOM   1   OM   123,617   1987,730   353,624
9   M1   PLOM   1   OM   46,399   943,421   121,490
12   M1   PLOM   1   OM   61,940   1031,544   104,752
15   M1   PLOM   1   OM   93,734   1074,952   67,440
18   M1   PLOM   1   OM   46,264   536,972   236,963
24   M1   PLOM   1   OM   NA   NA   NA
0   M8   PLOM   1   OM   129,373   284,907   28,964
1   M8   PLOM   1   OM   159,159   416,556   50,974
3   M8   PLOM   1   OM   93,810   696,399   25,254
6   M8   PLOM   1   OM   324,836   927,455   541,007
9   M8   PLOM   1   OM   197,709   549,647   200,049
12   M8   PLOM   1   OM   115,324   187,850   241,379
15   M8   PLOM   1   OM   151,188   29,661   300,188
18   M8   PLOM   1   OM   103,133   45,595   333,985
24   M8   PLOM   1   OM   87,757   430,743   65,627
0   M17   PLOM   1   OM   106,597   302,356   35,955
1   M17   PLOM   1   OM   149,141   501,711   43,071
3   M17   PLOM   1   OM   83,811   857,553   121,024
6   M17   PLOM   1   OM   107,365   1222,996   445,976
9   M17   PLOM   1   OM   112,977   378,703   180,375
12   M17   PLOM   1   OM   95,837   163,410   267,944
15   M17   PLOM   1   OM   141,899   982,256   160,975
18   M17   PLOM   1   OM   50,576   123,574   198,350
24   M17   PLOM   1   OM   67,076   257,275   88,957

zop<-read.table("zoop2.txt",h=T)
library(lattice)
trellis.device(color=FALSE)
xyplot( GRANDHERB~ date | traitement, data = zop,    panel = function(x, y){ panel.grid() ;  panel.xyplot(x, y) ;panel.loess(x, y) })



On voit alors clairement que les profils du temps varient avec les traitements.

En vous remerciant par avance,

[Renaud Lancelot]

Pour les interactions temps * traitement, c'est qque chose de fréquent. Le pb est que votre fct du temps est complexe. Je ne vois pas de solution "facile".

1. Trouver un modèle dynamique décrivant l'évolution de votre système (logistique,...) et utiliser la régression non linéaire pour en estimer les paramètres. Il faudrait discuter de cela avec des spécialistes de la question.

2. Rester sur du modèle linéaire, et discrétiser le temps en pas où la croissance est ~ constante. Pb: vous risquez d'avoir un grand nb de paramètres surtout si vous avez des interactions.

3. Toujours en modèle linéaire, utiliser un polynôme du temps de degré suffisamment élevé pour rendre compte des variations complexes. Pas toujours simple. Dans votre cas, ça marche plutôt bien même avec un polynôme d'ordre 3: cf ci-dessous.

NB: dans chaque cas, le temps est considéré comme un effet fixe. Sauf à montrer que vos mésocosmes évoluent différemment pour les mêmes valeurs des effets fixes, il n'y a pas lieu d'introduire d'effet aléatoire.

Pour les résultats d'anova, la fct "anova" donne la somme des carrés de type I. La somme des carrés de type III est très critiquée par les spécialistes (ne correspond à aucune hypothèse nulle quand le plan est désiquilibré). Si vous voulez mettre en évidence un effet de l'interaction (ou autre variable) il vaut mieux comparer les modèles par un test du F.

Voici un exemple avec vos données dont j'ai placé une copie ici. Adapter la commande read.table() ci-dessous pour le reproduire.

Code : Tout sélectionner

> zop <- read.table("d:/analyses/travail/data/rondel2.txt", header = TRUE)
> library(lattice)
> trellis.device(color = FALSE)
> xyplot(lGRANDHERB ~ date | traitement,
+        data = zop,
+        panel = function(x, y){
+           panel.grid()
+           panel.xyplot(x, y)
+           panel.loess(x, y)
+           })

> xyplot(log(GRANDHERB) ~ date | traitement,
+        data = zop,
+        panel = function(x, y){
+           panel.grid()
+           panel.xyplot(x, y)
+           panel.loess(x, y)
+           })
>
> fm1 <-lm(log(GRANDHERB) ~ poly(date, 3) * traitement, data = zop)
> summary(fm1)

Call:
lm(formula = log(GRANDHERB) ~ poly(date, 3) * traitement, data = zop)

Residuals:
     Min       1Q   Median       3Q      Max
-3.69703 -0.54435  0.02171  0.47869  3.11116

Coefficients:
                              Estimate Std. Error t value Pr(>|t|)   
(Intercept)                     4.0376     0.2267  17.811  < 2e-16 ***
poly(date, 3)1                -11.4944     2.8852  -3.984 0.000111 ***
poly(date, 3)2                 12.0668     2.8852   4.182 5.18e-05 ***
poly(date, 3)3                  2.6939     2.8852   0.934 0.352144   
traitementALPL                  0.1942     0.3288   0.591 0.555718   
traitementOM                    1.4835     0.3243   4.575 1.07e-05 ***
traitementPL                    1.5285     0.3206   4.768 4.79e-06 ***
traitementPLOM                  1.9500     0.3250   6.001 1.73e-08 ***
traitementsans                  2.2174     0.3206   6.917 1.72e-10 ***
poly(date, 3)1:traitementALPL  -0.5047     4.3241  -0.117 0.907250   
poly(date, 3)2:traitementALPL  -2.2540     4.2656  -0.528 0.598092   
poly(date, 3)3:traitementALPL   1.0410     4.2292   0.246 0.805951   
poly(date, 3)1:traitementOM     6.0793     4.1556   1.463 0.145832   
poly(date, 3)2:traitementOM   -13.1288     4.1475  -3.165 0.001917 **
poly(date, 3)3:traitementOM     1.5540     4.1452   0.375 0.708337   
poly(date, 3)1:traitementPL    15.2950     4.0804   3.748 0.000264 ***
poly(date, 3)2:traitementPL   -15.9930     4.0804  -3.920 0.000141 ***
poly(date, 3)3:traitementPL    -5.8663     4.0804  -1.438 0.152849   
poly(date, 3)1:traitementPLOM   8.3068     4.2666   1.947 0.053634 . 
poly(date, 3)2:traitementPLOM -13.3911     4.2636  -3.141 0.002074 **
poly(date, 3)3:traitementPLOM   3.0133     4.1391   0.728 0.467876   
poly(date, 3)1:traitementsans  12.4846     4.0804   3.060 0.002677 **
poly(date, 3)2:traitementsans -14.5869     4.0804  -3.575 0.000488 ***
poly(date, 3)3:traitementsans  -3.3015     4.0804  -0.809 0.419874   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.178 on 134 degrees of freedom
Multiple R-Squared: 0.5405,     Adjusted R-squared: 0.4616
F-statistic: 6.853 on 23 and 134 DF,  p-value: 1.250e-13

> # data.frame pour les prédictions
> New <- expand.grid(date = 0:24, traitement = levels(zop$traitement))
> pred <- predict(fm1, newdata = New, se = TRUE)
> New$fit <- pred$fit
> New$se <- pred$se.fit
>
> graphics.off()
> trellis.device(color = FALSE)
> # on utilise subscripts pour passer l'index des valeurs de chaque panel
> xyplot(log(GRANDHERB) ~ date | traitement,
+        data = zop, subscripts = TRUE,
+        panel = function(x, y, subscripts){
+           panel.grid()
+           panel.xyplot(x, y)
+           panel.loess(x, y)
+           xtraitement <- unique(zop$traitement[subscripts])
+           llines(New$date[New$traitement == xtraitement], New$fit[New$traitement == xtraitement], col = "red")
+           })
>
> fm2 <-lm(log(GRANDHERB) ~ poly(date, 3) + traitement, data = zop)
> anova(fm1, fm2, test = "F")
Analysis of Variance Table

Model 1: log(GRANDHERB) ~ poly(date, 3) * traitement
Model 2: log(GRANDHERB) ~ poly(date, 3) + traitement
  Res.Df     RSS  Df Sum of Sq      F    Pr(>F)   
1    134 185.917                                   
2    149 268.114 -15   -82.197 3.9496 7.252e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


[Marianne Robert]

bonjour a tous

je souhaiterez avoir un peut plus d info les critiques faites aux sommes de carrés des écarts de type 3.

je souhaite réaliser une manova afin de savoir l effet du type de banc (variables independante (IVs) à 2 modalités L ou O) a un effet sur la combinaison des mes variables dépendantes (DVs) que sont TAG, ST et PL

Code : Tout sélectionner

   Type_Banc TAG  ST  PL
23         L 0.2 0.1 2.0
24         L 0.6 0.1 1.7
25         L 0.4 0.1 4.1
26         O 0.9 0.1 8.0
27         O 0.7 0.2 2.5
28         O 0.7 0.2 1.9


neamoins, le plan d expérience ne compte pas le même nombre d'individus dans chaque modalité du type de banc ,
j ai 25 données pour le banc L et 30 pour la modalité O


je pensais qu'une des somution permettant de prendre en compte cela était de passer d'une somme des carrés des écarts de type 1 à 3.

Code : Tout sélectionner

> Y1 <- cbind(list$TAG, list$PL)
> summary(manova(Y1 ~ Sexe+Type_Banc+ Sexe*Type_Banc, data=list),test="Wilks", type="III")
               Df   Wilks approx F num Df den Df  Pr(>F) 
Sexe            3 0.76908  2.10431      6     90 0.06036 .
Type_Banc       1 0.89879  2.53355      2     45 0.09065 .
Sexe:Type_Banc  3 0.77622  2.02547      6     90 0.07029 .
Residuals      46                                         
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>


merci a vous

[Marianne Robert]

erratum sur le dernier code

Code : Tout sélectionner

YY <- lm(Y1 ~Sexe+Type_Banc+ Sexe:Type_Banc, data=list)
Manova(YY,test="Wilks", type="III",multivariate=T)


erreur de copier coller

merci

[Marianne Robert]

Code : Tout sélectionner

Y1 <- cbind(list$TAG, list$PL)

YY <- lm(Y1 ~-1+Sexe+Type_Banc+ Sexe:Type_Banc, data=list)
> Manova(YY,test="Wilks", type="III",multivariate=T)

Type III MANOVA Tests: Wilks test statistic
               Df test stat approx F num Df den Df    Pr(>F)   
Sexe            4    0.2979   9.3613      8     90 2.577e-09 ***
Type_Banc       1    0.8556   3.7985      2     45   0.02990 * 
Sexe:Type_Banc  3    0.7762   2.0255      6     90   0.07029 . 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>


le -1 permet de retirer l interspect qui pour moi n a pas d utilité (pas d interpretation biologique)