Interprétation de la sortie d'une ANOVA à trois facteurs

Bonjour,
Pour une analyse je dois faire une ANOVA avec trois facteurs:
variable: nombre d'individus
facteur 1 : station (14 stations par années)
facteur 2: habitat (4 habitats, dans lesquels sont regroupés les stations; ex: dans l'habitat FR on a les stations FRa, FRb, FRc et FRd)
facteur 3: année d'échantillonnage (sur trois ans)

Utilisation de lm () car j'ai cru comprendre que c'était plus adéquate que aov() , plus tolérant dans le traitement de données . Mon soucis est que dans tous les cas, il me manque une interaction dans les deux dernières tentatives d'ANOVA.
soit il me manque l'interaction "habitats:stations" soit "habitats:ans"

Code:

> summary(ANOV_pocilloT)
    recrues         stations   habitats     annees    
 Min.   :0.000   ax     : 60   BR:116   Min.   :2012  
 1st Qu.:0.000   bx     : 60   FR:234   1st Qu.:2012  
 Median :1.099   by     : 60   MS:115   Median :2013  
 Mean   :1.042   frc    : 60   OS:356   Mean   :2013  
 3rd Qu.:1.609   az     : 59            3rd Qu.:2014  
 Max.   :4.094   brb    : 59            Max.   :2014  
                 (Other):463                          
> stations=as.factor(ANOV_pocilloT$stations)
> habitats=as.factor(ANOV_pocilloT$habitats)
> ans=as.factor(ANOV_pocilloT$annees)
> anova(lm(recrues~habitats+stations+ans,data=ANOV_pocilloT))
Analysis of Variance Table

Response: recrues
           Df Sum Sq Mean Sq F value    Pr(>F)    
habitats    3  85.46 28.4866  60.684 < 2.2e-16 ***
stations   10 218.74 21.8743  46.598 < 2.2e-16 ***
ans         2  41.33 20.6648  44.021 < 2.2e-16 ***
Residuals 805 377.89  0.4694                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> anova(lm(recrues~habitats*ans*stations,data=ANOV_pocilloT))
Analysis of Variance Table

Response: recrues
              Df  Sum Sq Mean Sq F value    Pr(>F)    
habitats       3  85.460 28.4866 76.6131 < 2.2e-16 ***
ans            2  40.598 20.2989 54.5927 < 2.2e-16 ***
stations      10 219.475 21.9475 59.0266 < 2.2e-16 ***
habitats:ans   6  53.157  8.8594 23.8270 < 2.2e-16 ***
ans:stations  20  35.082  1.7541  4.7176 7.183e-11 ***
Residuals    779 289.651  0.3718                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> anova(lm(recrues~habitats:stations+stations:ans+habitats:ans,data=ANOV_pocilloT))
Analysis of Variance Table

Response: recrues
                   Df Sum Sq Mean Sq F value    Pr(>F)    
habitats:stations  13 304.20 23.4002  62.934 < 2.2e-16 ***
stations:ans       28 129.57  4.6275  12.445 < 2.2e-16 ***
Residuals         779 289.65  0.3718                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

mais la c'est un problème que j'ai soumis a un autre forum spécialisé dans le logiciel R^^

et quand j'essaie de voir les différences dans le détails avec tukeyHSD est bien R me dit que ça ne fonctionne pas avec lm().... donc je pense tester un SNK (mais je ne sais pas si c'est mieux qu'un Tukey).

bref je tente avec aov():

Code:

> r=aov(recrues~habitats*ans*stations,data=ANOV_pocilloT)
> summary(r)
              Df Sum Sq Mean Sq F value   Pr(>F)    
habitats       3  85.46  28.487  76.613  < 2e-16 ***
ans            2  40.60  20.299  54.593  < 2e-16 ***
stations      10 219.47  21.947  59.027  < 2e-16 ***
habitats:ans   6  53.16   8.859  23.827  < 2e-16 ***
ans:stations  20  35.08   1.754   4.718 7.18e-11 ***
Residuals    779 289.65   0.372                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> r2=aov(recrues~habitats:stations+stations:ans+habitats:ans,data=ANOV_pocilloT)
> r2
Call:
   aov(formula = recrues ~ habitats:stations + stations:ans + habitats:ans,
    data = ANOV_pocilloT)

Terms:
                habitats:stations stations:ans Residuals
Sum of Squares           304.2029     129.5686  289.6509
Deg. of Freedom                13           28       779

Residual standard error: 0.6097737
49 out of 91 effects not estimable
Estimated effects may be unbalanced


A part le fait qu'il me manque encore une interaction; autant en r qu'en r2; je ne comprend pas le message disant "49 out of 91 effects not estimable Estimated effects may be unbalanced" ????

quelques suggestions serait les bienvenues

Cordialement,
Camille

Vous avez posé cette question également sur le forum francophone sur l'usage de R. Il faut savoir, c'est un problème statistique, ou un problème sur l'usage de R. Je vous ai répondu sur l'autre forum. Je pense que votre jeu de données n'est pas équilibré "unbalanced". Il vous faut lire, je pense, des cours de base sur l'ANOVA.

HTH, Eric.

Bonjour,
oui je sais que mes données ne sont pas "pertinentes" à l'utilisation d'une ANOVA ; mais je l'ai tout de même fait car cela m'était demandé :/ pression de la hiérarchie disons

Bref, pour une utilisation et surtout une compréhension plus personnelle, je voudrais traiter ce jeu de données avec GLM ce qui , je pense, serait plus approprié. Toujours dans l'optique de mettre en évidence les différences significatives de ma variable en fonction des trois facteurs (les stations étant imbriquées dans habitats).

Suis-je sur la bonne voix, si j'utilise GML avec la loi de Poisson puisqu'il s'agit d'un comptage? j'ai aussi transformé mes données en log x+1 pour approcher la loi normale (suffisant pour seulement quelques stations mais pas pour l'analyse globale :/ )

Je suis désolée de vous imposer mon ignorance , ce doit être évident pour vous; mais mes lectures me font douter de tout...

Cordialement
Camille

Bonjour,

Il n'y a pas grand chose d'évident pour qui que ce soit, surtout quand on est dans un domaine pour lequel aucune analyse stat ne convient parfaitement, rassure toi.

mais je l'ai tout de même fait car cela m'était demandé :/ pression de la hiérarchie

Dans ce cas, il est toujours utile de montrer à la "hiérarchie" les incohérences éventuelles (voire résultats faux) d'une méthode inappropriée et si résultat identique alors l'argument est qu'il vaut toujours mieux être plus robuste sur les bases d'analyse que de produire un résultat à l'interprétation juste mais dont les bases d'analyse sont fausse.

Un Glm Poisson semble être un choix adapté vu le descriptif de tes données. Les interactions qui "disparaissent" sont souvent dues à des effets non testables étant donnée le plan d'expérience.

Qu'as tu transformé en log ? La famille d'erreur de type Poisson dans R contient par défaut un lien log qui permet de gérer pas mal de problème dont la non-linéarité de la distribution. Il n'y a pas d'histoire de normalité puisque c'est un comptage. En plus si tu entre une variable log-transformée dans un modèle de Poisson qui attend des entiers en entrée, tu vas générer des problèmes pour rien.

Nik

Merci de la rapidité de réponse,
j'ai testé glm() avec la loi de poisson sur mes données brutes et voila ce que j'obtiens:

Code:

> ANOV_taxa=read.table(file.choose(),sep=";",header=T)
> ANOV_taxa=ANOV_taxa[complete.cases(ANOV_taxa),]
> ANOV_taxa
    recrues stations habitats annees
1        13      fra       FR   2012
2        14      fra       FR   2012
3        28      fra       FR   2012
4        19      fra       FR   2012
(...)
837       4       bz       OS   2014
838       7       bz       OS   2014
839       5       bz       OS   2014
840       3       bz       OS   2014
> ANOV_taxa$recrues<-as.numeric(ANOV_taxa$recrues)
> summary(ANOV_taxa)
    recrues          stations   habitats     annees    
 Min.   :  0.00   ax     : 60   BR:116   Min.   :2012  
 1st Qu.:  4.00   bx     : 60   FR:234   1st Qu.:2012  
 Median : 10.00   by     : 60   MS:115   Median :2013  
 Mean   : 34.35   frc    : 60   OS:356   Mean   :2013  
 3rd Qu.: 24.00   az     : 59            3rd Qu.:2014  
 Max.   :697.00   brb    : 59            Max.   :2014  
                  (Other):463                          
> station=as.factor(ANOV_taxa$stations)
> habitat=as.factor(ANOV_taxa$habitats)
> an=as.factor(ANOV_taxa$annees)
> station
Levels: ax ay az bra brb bx by bz fra frb frc frd msa msb
> habitat
 Levels: BR FR MS OS
> an
Levels: 2012 2013 2014

> res=glm(ANOV_taxa$recrues~habitat*an,family=poisson)
> summary(res)

Call:
glm(formula = ANOV_taxa$recrues ~ habitat * an, family = poisson)

Deviance Residuals:
     Min        1Q    Median        3Q       Max  
-15.6883   -2.3862   -0.9306    1.2489   27.7553  

Coefficients:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)       2.09495    0.05547  37.768  < 2e-16 ***
habitatFR         0.99040    0.06040  16.397  < 2e-16 ***
habitatMS         0.03755    0.07821   0.480    0.631    
habitatOS         0.71746    0.05981  11.996  < 2e-16 ***
an2013            0.44666    0.07103   6.288 3.21e-10 ***
an2014            0.29534    0.07498   3.939 8.18e-05 ***
habitatFR:an2013 -1.35180    0.08395 -16.103  < 2e-16 ***
habitatMS:an2013  1.05302    0.09361  11.249  < 2e-16 ***
habitatOS:an2013 -1.54292    0.08400 -18.368  < 2e-16 ***
habitatFR:an2014  1.89361    0.07913  23.931  < 2e-16 ***
habitatMS:an2014  2.40565    0.09421  25.534  < 2e-16 ***
habitatOS:an2014 -1.22354    0.08614 -14.204  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 65302  on 820  degrees of freedom
Residual deviance: 17165  on 809  degrees of freedom
AIC: 20611

Number of Fisher Scoring iterations: 5



est-ce correct? ou dois je utiliser la fonction : family=quasipoisson?
De plus à part les facteurs seuls, je ne vois pas comment interpréter les résultats, avec:

Code:

habitatFR:an2013 -1.35180    0.08395 -16.103  < 2e-16 ***
habitatMS:an2013  1.05302    0.09361  11.249  < 2e-16 ***
habitatOS:an2013 -1.54292    0.08400 -18.368  < 2e-16 ***
habitatFR:an2014  1.89361    0.07913  23.931  < 2e-16 ***
habitatMS:an2014  2.40565    0.09421  25.534  < 2e-16 ***
habitatOS:an2014 -1.22354    0.08614 -14.204  < 2e-16 ***

et l'habitat BR ainsi que l'année 2012 ne sont pas présents, puis-je tout de même émettre une interprétation?

Merci

En plus de faire un summary() sur l'objet GLM, il serait intéressant de le mouliner avec la fonction anova() et l'argument "test="Chisq":

Code:

anova(res, test="Chisq")

Pour ce qui est de quasipoisson ou non, les estimateurs des pentes resterons les mêmes, ce sont leur SE qui vont changer. L'idée est de ré-ajuster le modèle avec quasipoisson et voir la valeur du "Despersion parameter" dans la sortie de summary() qui devrait idéalement être proche de 1.0.

HTH, Eric.

Bonjour,

ici clairement tes données sont sur-dispersées et il te faut t'orienter vers du quasipoisson ou de la négative binomiale (fonction glm.nb de la librairie MASS). Que l'habitat BR et l'année 2012 n'apparaissent pas c'est normal. Comme ils sont les premiers par ordre alphabétique, ils servent de références. Les autres paramètres sont là pour te dire ce qu'il faut ajouetr ou retrancher si tu es d'une autre année ou d'un autre habitat (par rapport à 2012 et BR).

cdlt

Très bien ,
donc d'après ce que vous me dites j'ai testé :

Code:

> anova(res,test="Chisq")
Analysis of Deviance Table

Model: poisson, link: log

Response: ANOV_taxa$recrues

Terms added sequentially (first to last)


           Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                         820      65302              
habitat     3  20514.3       817      44788 < 2.2e-16 ***
an          2  19433.4       815      25354 < 2.2e-16 ***
habitat:an  6   8189.5       809      17165 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


j'obtiens de la variations spatiale et temporelle de mes individus; ce qui es cohérent


Ensuite je fais le même test mais avec "quasipoisson" pour voir:

Code:

> res1=glm(ANOV_taxa$recrues~habitat*an,family=quasipoisson)
> summary(res1)

Call:
glm(formula = ANOV_taxa$recrues ~ habitat * an, family = quasipoisson)

Deviance Residuals:
     Min        1Q    Median        3Q       Max  
-15.6883   -2.3862   -0.9306    1.2489   27.7553  

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       2.09495    0.26387   7.939 6.78e-15 ***
habitatFR         0.99040    0.28733   3.447 0.000596 ***
habitatMS         0.03755    0.37203   0.101 0.919628    
habitatOS         0.71746    0.28452   2.522 0.011871 *  
an2013            0.44666    0.33789   1.322 0.186576    
an2014            0.29534    0.35666   0.828 0.407863    
habitatFR:an2013 -1.35180    0.39934  -3.385 0.000746 ***
habitatMS:an2013  1.05302    0.44531   2.365 0.018280 *  
habitatOS:an2013 -1.54292    0.39960  -3.861 0.000122 ***
habitatFR:an2014  1.89361    0.37640   5.031 6.02e-07 ***
habitatMS:an2014  2.40565    0.44817   5.368 1.04e-07 ***
habitatOS:an2014 -1.22354    0.40978  -2.986 0.002913 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasipoisson family taken to be 22.62864)

    Null deviance: 65302  on 820  degrees of freedom
Residual deviance: 17165  on 809  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 5

> anova(res1,test="Chisq")
Analysis of Deviance Table

Model: quasipoisson, link: log

Response: ANOV_taxa$recrues

Terms added sequentially (first to last)


           Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                         820      65302              
habitat     3  20514.3       817      44788 < 2.2e-16 ***
an          2  19433.4       815      25354 < 2.2e-16 ***
habitat:an  6   8189.5       809      17165 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Ce qui change les résultats.... mais ayant des données sur-dispersées , je pense que je devrais favoriser le glm() avec quasipoisson ,non?

Enfin pour l'interprétation; voire si j'ai bien compris, pour :

Code:

habitatOS         0.71746    0.28452   2.522 0.011871 *  
an2013            0.44666    0.33789   1.322 0.186576  
(...)  
habitatFR:an2013 -1.35180    0.39934  -3.385 0.000746 ***

-habitat OS est différent significativement par rapport à la référence qu'est l'habitat
BR
- l'année 2013 n'est pas significativement différente de la référence qu'est 2012
- l'habitat FR en 2013 est différent de ??? BR en 2012 ??

je crois que je me trompe :/