différences entre aov et lme

[caroline domerg]

Bonjour,

Dans le cadre d'un split-plot pour le relevé d'un rendement quelconque, nous avons les facteurs suivants:

B : 3 blocks
I : 2 types de parcelle, une irriguée et l'autre pas
V : 2 types de sous-parcelle correspondant à deux variétés.
N : 3 niveaux d'azote

Chaque block est divisé en 2 parcelles (une irriguée et une pluviale) . Chacune de ces parcelles est redivisée en deux sous-parcelles (variétés différentes) et chaque sous-parcelle est partitionnée en 3 niveaux d'azote.

Dans Venables & Ripley (Modern Applied Statistics with S-plus, 3rd edition), un modèle multistrates mais pour trois facteurs est ajusté pour le même type de données:
aov(Y~N*V+Error(B/V), qr=T)

nous l'adaptons à nos données:
aov(Y~N*V*I+Error(B/V/I))

D'après Pinheiro & Bates (Mixed-effect models in S and S-plus) et comme nous l'avons vu dans le message "Re: lme and lmer syntax de Ronaldo Reis-Jr, Wed 26 Oct 2005", nous ajustons également le modèle mixte suivant:
lme(Y~N*V*I, random~1|B/V/I)

Sur un Y genéré aléatoirement, on a tout juste 1 facteur avec le meme test F pour les 2 modèles proposés, pour les autres certains ont un ordre de grandeur équivalent et d'autres sont complétement différents (alors qu'on retrouve des valeurs similaires dans l'exemple présenté dans Venables & Ripley et Pinheiro & Bates)

Y-a-t-il une erreur dans un de nos deux modèles ou une explication à ces différences?

Merci d'avance.

Caroline Domerg and Frédéric Chiroleu
UMR 53 PVBMT (Peuplements Vegetaux et Bio-agresseurs en Milieu Tropical)
CIRAD
Pôle de Protection des Plantes (3P) - Saint-Pierre
Ile de la Réunion

[Renaud Lancelot]

* Si votre plan d'observation est équilibré, les deux modèles devraient fournir les mêmes résultats. Si ce n'est pas le cas, il est normal qu'il y ait des différences dans les coefficients ajustés, et il faut alors considérer les résultats de lme (ceux de aov sont alors biaisés).

* La fct lme fournit par défaut des estimateurs REML alors que aov fournit (je crois !...) des estimateurs ML. Utuliser l'argument "method" pour obtenir des estimateurs ML:

Code : Tout sélectionner

fm <- lme(Y ~ N * V * I, random ~ 1 | B / V / I, method = "ML")

* La table d'analyse de variance fournie par anova:

Code : Tout sélectionner

fm <- lme(Y ~ N * V * I, random ~ 1 | B / V / I)
anova(fm)

donne les résultats d'une anova séquentielle. Pas sûr que ce soit le cas avec aov (je n'utilise pas cette fct car "mes" plans d'obs st tjs déséquilibrés ;-)). Si vous voulez les résultats d'une anova marginale, il faut utiliser l'argument "type":

Code : Tout sélectionner

anova(fm, type = "m")


Good luck !

Renaud

[caroline domerg]

Merci Renaud, on a effectivement un plan équilibré (pour notre simulation, en attendant les données) et entre la méthode REML et ML, les valeurs des tests F ne bougent pas...

Une erreur aussi dans mon message mais pas dans les modèles, il s'agit de B/I/V et pas de B/V/I.

L'anova marginale sur notre modèle mixte ne nous permet pas non plus de retrouver les résultats de l'aov.

Sans vouloir nous décharger de la tâche, voici les lignes de commande que nous utilisons puisqu'apparemment il y a quelque chose qui nous échappe.

Code : Tout sélectionner

library("nlme")
   
Y1<-rnorm(36)
bl<-factor(rep(1:3, times=12))
ir<-factor(rep(1:2, each=18))
va<-factor(rep(1:2,each=3, times=6)) 
az<-factor(rep(1:3, each=3, times=4))
   
m.aov<-aov(Y1~ir*va*az + Error(bl/ir/va))
m.lme<-lme(Y1~ir*va*az, random=~1|bl/ir/va)
   
summary(m.aov)
anova(m.lme, type="m")
anova(m.lme, type="s")


Encore merci!!

Caroline

[Renaud Lancelot]

Pas le temps de regarder cela en détail tout de suite. Cependant une remarque de fond: I, V et N sont typiquement des effets fixes. Pourquoi vouloir les considérer comme aléatoires ? Pour ma part, je me contenterais d'un effet aléatoire "Bloc". Je me simplifierais la vie en numérotant les blocs de 1 à n (contrairement à ce que tu as fait dans le code). Avec lme, le modèle serait alors:

Code : Tout sélectionner

Bl <- paste(ir, va, az, sep = "_")
m.lme <- lme(Y1 ~ ir*va*az, random = ~ 1 | Bl)


Je n'ai pas vérifié. Je regarderai ce soir.

A+

Renaud

[Renaud Lancelot]

Les modèles sont trop complexes par rapport aux données:

Code : Tout sélectionner

> library("nlme")
> set.seed(123)
> Data <- data.frame(
+   Y1 = rnorm(36),
+   bl = factor(rep(1:3, times = 12)),
+   ir = factor(rep(1:2, each = 18)),
+   va = factor(rep(1:2, each = 3, times = 6)),
+   az = factor(rep(1:3, each = 3, times = 4)))
> Data$bloc <- with(Data, factor(paste(ir, va, az, sep = "")))
>
> m1 <- aov(Y1 ~ ir * va * az + Error(bl / ir / va), data = Data)
> m2 <- lme(Y1 ~ ir * va * az, random = ~ 1 | bl / ir / va, data = Data)
>
> m1

Call:
aov(formula = Y1 ~ ir * va * az + Error(bl/ir/va), data = Data)

Grand Mean: 0.05560446

Stratum 1: bl

Terms:
                Residuals
Sum of Squares   2.019690
Deg. of Freedom         2

Residual standard error: 1.004910

Stratum 2: bl:ir

Terms:
                       ir Residuals
Sum of Squares  0.2855233 1.4188413
Deg. of Freedom         1         2

Residual standard error: 0.8422711
5 out of 6 effects not estimable
Estimated effects are balanced

Stratum 3: bl:ir:va

Terms:
                       va     ir:va Residuals
Sum of Squares  1.0369324 0.1304834 1.0626924
Deg. of Freedom         1         1         4

Residual standard error: 0.5154349
4 out of 6 effects not estimable
Estimated effects may be unbalanced

Stratum 4: Within

Terms:
                       az     ir:az     va:az  ir:va:az Residuals
Sum of Squares   0.114185  1.249827  1.915827  2.315214 19.306785
Deg. of Freedom         2         2         2         2        16

Residual standard error: 1.098487
Estimated effects may be unbalanced


Noter les avertissements sur les effets non estimables.

Code : Tout sélectionner

> summary(m2)
Linear mixed-effects model fit by REML
 Data: Data
       AIC      BIC    logLik
  113.0993 131.9481 -40.54964

Random effects:
 Formula: ~1 | bl
        (Intercept)
StdDev:   0.0402903

 Formula: ~1 | ir %in% bl
         (Intercept)
StdDev: 0.0001522510

 Formula: ~1 | va %in% ir %in% bl
         (Intercept)  Residual
StdDev: 7.529249e-05 0.9951769

Fixed effects: Y1 ~ ir * va * az
                 Value Std.Error DF    t-value p-value
(Intercept)  0.2560184 0.5750363 16  0.4452212  0.6621
ir2         -0.5357715 0.8125585  2 -0.6593635  0.5774
va2          0.1233928 0.8125585  4  0.1518572  0.8867
az2         -0.2708140 0.8125585 16 -0.3332856  0.7432
az3         -0.7530177 0.8125585 16 -0.9267242  0.3678
ir2:va2      0.2460439 1.1491313  4  0.2141130  0.8409
ir2:az2      0.8927399 1.1491313 16  0.7768824  0.4486
ir2:az3      0.5414556 1.1491313 16  0.4711869  0.6439
va2:az2      0.5296899 1.1491313 16  0.4609481  0.6510
va2:az3      0.4796553 1.1491313 16  0.4174069  0.6819
ir2:va2:az2 -1.8989230 1.6251170 16 -1.1684838  0.2597
ir2:va2:az3  0.4383414 1.6251170 16  0.2697291  0.7908
[...]
Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max
-2.08535315 -0.55331868 -0.04829478  0.46672070  1.68519531

Number of Observations: 36
Number of Groups:
                bl         ir %in% bl va %in% ir %in% bl
                 3                  6                 12
> VarCorr(m2)
            Variance     StdDev     
bl =        pdLogChol(1)             
(Intercept) 1.623309e-03 4.029030e-02
ir =        pdLogChol(1)             
(Intercept) 2.318036e-08 1.522510e-04
va =        pdLogChol(1)             
(Intercept) 5.668959e-09 7.529249e-05
Residual    9.903770e-01 9.951769e-01

Noter les valeurs très faibles (proches de 0) de plusieurs variances des effets aléatoires. L'algorithme semble plus performant, mais ce n'est qu'une apparence:

Code : Tout sélectionner

> intervals(m2)
Erreur dans intervals.lme(m2) : Cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance


Les problèmes subsistent avec des modèles aléatoires plus simples:

Code : Tout sélectionner

> m3 <- aov(Y1 ~ ir * va * az + Error(bloc), data = Data)
> m4 <- lme(Y1 ~ ir * va * az, random = ~ 1 | bloc, data = Data)
> m3

Call:
aov(formula = Y1 ~ ir * va * az + Error(bloc), data = Data)

Grand Mean: 0.05560446

Stratum 1: bloc

Terms:
                       ir        va        az     ir:va     ir:az     va:az
Sum of Squares  0.2855233 1.0369324 0.1141853 0.1304834 1.2498272 1.9158266
Deg. of Freedom         1         1         2         1         2         2
                 ir:va:az
Sum of Squares  2.3152141
Deg. of Freedom         2

Estimated effects may be unbalanced

Stratum 2: Within

Terms:
                Residuals
Sum of Squares   23.80801
Deg. of Freedom        24

Residual standard error: 0.9959921
> summary(m3)

Error: bloc
         Df  Sum Sq Mean Sq
ir        1 0.28552 0.28552
va        1 1.03693 1.03693
az        2 0.11419 0.05709
ir:va     1 0.13048 0.13048
ir:az     2 1.24983 0.62491
va:az     2 1.91583 0.95791
ir:va:az  2 2.31521 1.15761

Error: Within
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals 24 23.808   0.992               
> summary(m4)
Linear mixed-effects model fit by REML
 Data: Data
       AIC      BIC    logLik
  109.0996 125.5924 -40.54982

Random effects:
 Formula: ~1 | bloc
        (Intercept)  Residual
StdDev:    1.533430 0.9959921

Fixed effects: Y1 ~ ir * va * az
                 Value Std.Error DF    t-value p-value
(Intercept)  0.2560184  1.637704 24  0.1563276  0.8771
ir2         -0.5357715  2.316063  0 -0.2313285     NaN
va2          0.1233928  2.316063  0  0.0532770     NaN
az2         -0.2708140  2.316063  0 -0.1169286     NaN
az3         -0.7530177  2.316063  0 -0.3251283     NaN
ir2:va2      0.2460439  3.275408  0  0.0751186     NaN
ir2:az2      0.8927399  3.275408  0  0.2725584     NaN
ir2:az3      0.5414556  3.275408  0  0.1653093     NaN
va2:az2      0.5296899  3.275408  0  0.1617172     NaN
va2:az3      0.4796553  3.275408  0  0.1464414     NaN
ir2:va2:az2 -1.8989230  4.632127  0 -0.4099462     NaN
ir2:va2:az3  0.4383414  4.632127  0  0.0946307     NaN
[...]
Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max
-2.08100636 -0.54993282 -0.04561536  0.46985780  1.68762809

Number of Observations: 36
Number of Groups: 12
Warning message:
production de NaN in: pt(q, df, lower.tail, log.p)
> VarCorr(m4)
bloc = pdLogChol(1)
            Variance  StdDev   
(Intercept) 2.3514082 1.5334302
Residual    0.9920003 0.9959921
> intervals(m4)
Erreur dans intervals.lme(m4) : Cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance
De plus : Warning message:
production de NaN in: qt(p, df, lower.tail, log.p)


Idem avec les modèles des composantes de la variance:

Code : Tout sélectionner

> m4 <- aov(Y1 ~ 1 + Error(bl / ir / va), data = Data)
> m5 <- lme(Y1 ~ 1, random = ~ 1 | bl / ir / va, data = Data, method = "ML")
>
> m4

Call:
aov(formula = Y1 ~ 1 + Error(bl/ir/va), data = Data)

Grand Mean: 0.05560446

Stratum 1: bl

Terms:
                Residuals
Sum of Squares   2.019690
Deg. of Freedom         2

Residual standard error: 1.004910

Stratum 2: bl:ir

Terms:
                Residuals
Sum of Squares   1.704365
Deg. of Freedom         3

Residual standard error: 0.7537384

Stratum 3: bl:ir:va

Terms:
                Residuals
Sum of Squares   2.230108
Deg. of Freedom         6

Residual standard error: 0.6096595

Stratum 4: Within

Terms:
                Residuals
Sum of Squares   24.90184
Deg. of Freedom        24

Residual standard error: 1.018615
> summary(m4)

Error: bl
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals  2 2.0197  1.0098               

Error: bl:ir
          Df  Sum Sq Mean Sq F value Pr(>F)
Residuals  3 1.70436 0.56812               

Error: bl:ir:va
          Df  Sum Sq Mean Sq F value Pr(>F)
Residuals  6 2.23011 0.37168               

Error: Within
          Df  Sum Sq Mean Sq F value Pr(>F)
Residuals 24 24.9018  1.0376               
> VarCorr(m5)
            Variance     StdDev     
bl =        pdLogChol(1)             
(Intercept) 1.651793e-09 4.064226e-05
ir =        pdLogChol(1)             
(Intercept) 4.983046e-11 7.059070e-06
va =        pdLogChol(1)             
(Intercept) 6.647528e-12 2.578280e-06
Residual    8.571111e-01 9.258030e-01
> intervals(m5)
Erreur dans intervals.lme(m5) : Cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance


En bref: il faudrait mieux définir ce que vous voulez analyser: effets fixes, composantes de la variance et choisir un modèle approprié. D'après ce que vous décrivez dans le jeu de données simulées, je pense qu'il faut partir d'un modèle additif des effets fixes, tester l'existence d'interactions entre ces effets fixes et n'utiliser un modèle mixte (effet aléatoire lié aux blocs) uniquement si cela s'avère nécessaire (forte corrélation intra-bloc).

Amicalement,

Renaud

[caroline domerg]

On vient de récupérer les données réelles et il s'avère que lorsqu'on compare les deux modèles suivants :

Code : Tout sélectionner

maov.06<- aov(Y~I*V*N + Error(B/I/V/N), data=spli06)
mlme.06<-lme(Y~I*V*N , random=~1| B/I/V/N, data=spli06)


on obtient les mêmes valeurs des tests F sauf pour le facteur irrigation mais avec une p-value du même ordre...

Code : Tout sélectionner


summary(maov.06)

Error: B
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals  2 98.676  49.338               

Error: B:I
          Df  Sum Sq Mean Sq F value Pr(>F)
I          1 1542.26 1542.26  3.9973 0.1836
Residuals  2  771.65  385.82               

Error: B:I:V
          Df  Sum Sq Mean Sq F value  Pr(>F) 
V          1 1724.05 1724.05 10.5974 0.03122 *
I:V        1    1.75    1.75  0.0107 0.92246 
Residuals  4  650.74  162.69                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Error: B:I:V:N
          Df  Sum Sq Mean Sq  F value    Pr(>F)   
N          2 20720.5 10360.2 138.8727 7.748e-11 ***
I:N        2   442.5   221.3   2.9659  0.080238 . 
V:N        2  1047.4   523.7   7.0199  0.006477 **
I:V:N      2    73.7    36.9   0.4943  0.619031   
Residuals 16  1193.6    74.6                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(mlme.06)
            numDF denDF   F-value p-value
(Intercept)     1    16 1285.1790  <.0001
I               1     2    7.0883  0.1169
V               1     4   10.5975  0.0312
N               2    16  138.8718  <.0001
I:V             1     4    0.0107  0.9225
I:N             2    16    2.9658  0.0802
V:N             2    16    7.0199  0.0065
I:V:N           2    16    0.4943  0.6190



Etonnant que l'on obtienne des résultats si différents entre les modèles ajustés sur les données simulées (aucune ifluence des facteurs sur réponse et/ou mauvais choix de paramètres pour générér Y)

Comme vous le disiez dans le dernier message, on va interpréter le modèle lm avec interactions entre effets fixes et erreur liée aux blocks.

De plus, j'ai trouvé une référence
http://www3.imperial.ac.uk/portal/pls/p ... 171923.PDF
dans laquelle ce type de données (split plot à 4 facteurs) était modélisé de la même façon.

Merci beaucoup pour le temps que vous avez consacré à notre problème.

Caroline

[Renaud Lancelot]

Etonnant que l'on obtienne des résultats si différents entre les modèles ajustés sur les données simulées (aucune ifluence des facteurs sur réponse et/ou mauvais choix de paramètres pour générér Y)

Pas étonnant: compte tenu de la manière dont vous avez généré les données, il n'y avait aucune corrélation entre Y et les effets fixes, et aucune corrélation intra-bloc ==> très difficile ou impossible d'estimer les effets fxes et aléatoires (d'où composantes de la variance très faibles, etc.).

Merci pour la réf.

Renaud