anova et critères de tri

[E.H. [compte supprimé]]

Bonjour à tous,

Je débute sous R, mais j'ai pas le choix du logiciel :-) J'ai lu les divers manuels de base, les cours de univ.lyon1, etc... mais je trouve pas réponse à mes questions.
Je suis donc confronté à un probleme.
Je dois faire des anova sur un ensemble de données importées sous R.

Ma base se présente ainsi :

Code : Tout sélectionner

> dim(focale)
[1] 507  20
> names(focale)
 [1] "Id"    "An"    "Mo"    "Ma"    "Li"    "He"    "Es"    "Ta"    "Co"    "Ty"    "Ag"    "Po"    "asrF1" "asrF2" "asrF3" "asrF4" "asrF5" "asrF6" "asrF7" "asrF8"

8 variables sont de type caractère à plusieurs niveaux (variables qualitatives). Par exemple :

Code : Tout sélectionner

> levels(focale$Ma)
[1] "B" "M"
> levels(focale$Es)
[1] "C" "M"
> levels(focale$He)
[1] "T1" "T2" "T3" "T4"

etc...

Les autres variables sont de type numérique (quantitatif).

Mes variables dépendantes seront successivement asrF1, asrF2, ...
Mes variables indépendantes sont toutes les autres.

---

Alors voilà les problèmes auquel je me confronte :

1)---

Je souhaite réaliser une anova de type :

Code : Tout sélectionner

> ana1 <- aov(asrF1~Ty, data=focale)

, mais en rajoutant des critères ! et c'est là que je suis perdu ... :'(

En clair, je souhaiterais par exemple faire une anova avec asrF1 dépendant, Ty indépendant, mais en choisissant comme critère l'un des niveaux de l'une des autres varibles (par explemple, faire cette anova sur les valeur où la variable focale$Es = "C", puis la même anova, mais cette fois avec la variable focale$Es = "M").
Comment faire cette anova avec critères ajoutés ?

Si ce n'est pas possible de rajouter ce type de critères directement dans l'anova, comment créer un nouvelle objet (focale_tri_Es_C) ne comprenant que les valeurs triées où focale$Es = "C", puis un autre (focale_tri_Es_M) ne comprenant que les valeurs triées où focale$Es = "M" ?

2)---

Là, c'est pour une aide stats tout court, jsuis un peu à la rue dans l'explication des résultats.

Sur cet exemple :

Code : Tout sélectionner

> ana1 <- aov(asrF1~Ty, data=focale); summary(ana1)

             Df Sum Sq Mean Sq F value   Pr(>F)   
Ty            4  3.098   0.775  3.9671 0.003521 **
Residuals   502 98.013   0.195                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Comment interpréter cela ?
Est ce : asrF1 ne diffère pas significativement selon le niveau de Ty pour un risque 0.35% ? ou alors : le risque de rejetter H0 alors qu'elle est vrai est de 0,35%; donc on rejette H0 et on accetpte l'hypo alternative (et les **, comment je l'interprète ?)

et si au contraire, j'avais eu Pr(>F) > 0.05, comment l'interpréter ?
Est ce : asrF1 differe selon le niveau de Ty avec risque erreur 5% ? Si c'est bon, comment savoir quel niveau de Ty est supérieur aux autres ?

---

Merci d'avance pour vos réponses ;-)

[Logez Maxime]

Bonjour,

J’avoue ne pas très bien comprendre ta première question mais je peux répondre à la deuxième question.
Ici asrF1 diffère significativement entre les niveaux de Ty.
Lorsque que tu pratiques un test statistique la statistique que tu calcules (ici la valeur de F) suit une certaine loi statistique sous l’hypothèse nulle H0 (qui est qu’il n’y a pas de vraie différence). Tu compares donc la valeur de la statistique obtenue avec une valeur théorique au seuil que tu fixes (5%, 1% ...). Dans le cas d’une loi de F de paramètre 4, 502 avec un seuil alpha = 5% la valeur théorique est de qf(0.95,4,502) = 2.389696. La valeur observée est supérieure à la valeur théorique : tu rejettes H0. La probabilité qui t’est donnée est la probabilité d’observer une valeur supérieure à celle que tu as obtenue selon une loi F de paramètre 4, 502 : donc P(F(4,502)>Fobs) = 0.003521. La probabilité d’avoir une valeur supérieure à la valeur seuil 2.389696 est de 5% (le seuil que tu te fixes). Ici la probabilité que tu obtiennes une valeur supérieure à 3.9671 est de 0.003521, ce qui est inférieure a 5% et même à 1% (c’est ce que signifie les deux étoiles trois étoiles pour 0.1%, ...), donc tu rejettes H0. Si au contraire Pr(>F) >5% alors la probabilité d’observer des valeurs qui suivent une loi F de paramètres … supérieures à ta valeur observée est supérieur au seuil alpha donc la valeur de F que tu as obtenues dans ce cas fait partie des valeurs possibles sous l’hypothèse H0, donc tu acceptes H0.

En espéant t'avoir aidé.

Maxime

[E.H. [compte supprimé]]

merci pour ce début de réponse :)

[Logez Maxime]

Bonjour,


Juste un petit complément d'informations, la p-value renvoie la probabilité de rejetter H0 alors que H0 est vraie (ce qui correspond à l'erreur de première espèce : "alpha"). Il convient donc de comparer cette probabilité avec le seuil qu'on se fixe. Si on se fixe un seuil de 1% et que le p-value est inférieure à ce seuil alors logiquement on rejette H0, si par contre cette p-value est supérieur au seuil alors on ne peut pas rejetter H0. J'espère que cela clarifera un peu les choses.

J'avoue ne pas comprendre tes autres questions, quels critères veux-tu rajouter dans tes anova, ne veux tu pas en fait changer la variable de référence dans l'anova ? Désolé mais ce n'est pas clair pour moi.

Maxime

[Nicolas Péru]

Pour répondre à la première question...

Lorsqu'il s'agit de rajouter des critères de sélection il faut utiliser les crochets "[...]"

Dans ton cas cela donne :

Code : Tout sélectionner

ana1 <- anova (lm(asrF1[focale$Es=="C"]~Ty[focale$Es=="C"], data=focale))


Note bien que le signe "=" se traduit dans un R par "==".

Ensuite tu peux combiner 2 critères comme ceci

Code : Tout sélectionner

ana1 <- anova (lm(asrF1[focale$Es=="C" & focale$He=="T1"]~Ty[focale$Es=="C" & focale$He=="T1"], data=focale))


Attention à bien spécifier le critère de sélection pour toutes les variables que tu utilises.

N.B : pour ma part je préfère utiliser la fonction "anova" plutôt que "aov" car je n'ai jamais bien compris ce que fait cette dernière. En plus en utilisant anova tu passes forcément par un objet de classe 'lm' qui te donnes des outils de validation de ton modèle complémentaire de la simple p-value.

[Logez Maxime]

Bonjour,


Je n'avais pas vi que focale$Es était une autre variable catégorielle que celles impliquées dans le modèle. Plutôt que faire des tests séparés sur focale$Es=="C" et ensuite sur focale$Es=="M", je te conseille d'intégrer cette variable directement dans l'anova et de pratiquer une anova a deux facteurs plutôt que faire deux anovas a la suite les unes des autres:

Code : Tout sélectionner

ana1 <- anova(lm(asrF1~Ty+ES) # pour un modèle additif
ana1 <- anova(lm(asrF1~Ty*ES) # pour un modèle avec interaction

Cela me semble plus juste.

Maxime

[Logez Maxime]

Re,


Pour ce qui est des comparaisons entre les moyennes des modalités d'un facteur, par exemple entre les modalités de Ty pour asfr1, tu as plusieurs possibilités de test, qui sont assez bien résumé ici : http://pbil.univ-lyon1.fr/R/liens/AvnerDea_bio.pdf

Une des premières choses à faire serait de faire un boxplot(asfr1~Ty) pour avoir une idéen de quelle(s) moyenne(s) sont susceptible d'être différentes.

Maxime

[Renaud Lancelot]

pour répondre à la première question...

Lorsqu'il s'agit de rajouter des critères de sélection il faut utiliser les crochets "[...]"

Dans ton cas cela donne :

Code : Tout sélectionner

ana1 <- anova (lm(asrF1[focale$Es=="C"]~Ty[focale$Es=="C"], data=focale))


Note bien que le signe "=" se traduit dans un R par "==".

Ensuite tu peux combiner 2 critères comme ceci

Code : Tout sélectionner

ana1 <- anova (lm(asrF1[focale$Es=="C" & focale$He=="T1"]~Ty[focale$Es=="C" & focale$He=="T1"], data=focale))
i


Aaargh ! Cette syntaxe est hyper-lourde et source d'erreur. il vaut mieux utiliser l'argument subset, qui existe dans la plupart des fcts de modélisation et évite d'autres désagréments (avec les niveaux de facteurs non représentés dans le subset, par exemple):

Code : Tout sélectionner

m <- lm(asrF1 ~ Ty, data = focale, subset = Es == "C" & He == "T1")
anova(m, test = "F")


Renaud

[E.H. [compte supprimé]]

je teste tout ca

Code : Tout sélectionner

> m <- lm(asrF1 ~ Ty, data = focale, subset = Es == "C")
> anova(m, test = "F")

Code : Tout sélectionner

Analysis of Variance Table

Response: asrF1
           Df Sum Sq Mean Sq F value    Pr(>F)   
Ty          4  4.403   1.101  6.8755 2.295e-05 ***
Residuals 414 66.281   0.160                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


on en conclue donc que Ty a une influence significative sur les variations de asrF1. dans ce cas, je suis intéressé pour évaluer l'influence des modalités de Ty (comparaison multiple)... Suis je obligé de réaliser les comparaisons 2 à 2( student), ou est il possible de le faire automatiquement pour l'ensemble des valeurs sous R ?

En revanche, y'a un truc que je saisis pas. Je mêne en parallèle les même analyses sous XLSTAT, et mes résultats sont différents en ce qui concerne les paramètresde modèle.
Voyez par vous même :

Sous R:

Code : Tout sélectionner

> summary(m)

Call:
lm(formula = asrF1 ~ Ty, data = focale, subset = Es == "C")

Residuals:
    Min      1Q  Median      3Q     Max
-0.3802 -0.3222 -0.1078  0.1297  1.2486

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  0.38024    0.08343   4.558 6.82e-06 ***
TyF         -0.27246    0.09012  -3.023  0.00265 **
TyFS        -0.18534    0.11089  -1.671  0.09540 . 
TyM         -0.05802    0.08794  -0.660  0.50981   
TyNI        -0.20148    0.12077  -1.668  0.09601 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4001 on 414 degrees of freedom
Multiple R-Squared: 0.06229,    Adjusted R-squared: 0.05323
F-statistic: 6.876 on 4 and 414 DF,  p-value: 2.295e-05


sous XLSTAT :

Code : Tout sélectionner

Paramètres du modèle :                  
                  
Source   Valeur   Ecart-type   t   Pr > |t|   
Constante      0,322     0,028   11,586   < 0,0001   
Ty-M              0,000       0,000            
Ty-F              -0,214   0,044   -4,877   < 0,0001   
Ty-FS             -0,127   0,078   -1,629   0,104   
Ty-E              0,058   0,088   0,660   0,510   
Ty-NI             -0,143   0,092   -1,566   0,118   

question : pourquoi mes estimations ne sont pas les mêmes sout XLSTAT et R, d'autant que R a exclu la valeur Ty-E ?

---

méthodes :

en ce qui concerne les anova, la littérature précise qu'elles ne sont valable que si la normalité et l'invariance des variances sont vérifiées.
dans mon cas, ce n'ets pas le cas. meme après transformation arcsin(racine(frequence)) , la normalité n'est pas vérifiée, et l'invariance des variance n'est pas vraie pour l'ensemble de mes données (mais elle l'est pour asrF1).

de ce fait, ai je vraiment le droit de réaliser une anova sur mes données ? mes résultats sont ils exploitables ?

merci a vous ;-)

[Nicolas Péru]

bonjour,

si tu n'as pas homoscedasticité des variances alors tu peux t'arreter là avec les anovas. Autant la normalité n'est pas toujours un critère fort surtout si ton échantillon est grand (peu de gens la teste) autant l'homoscedasticité est incontournable.

pour te répondre quant à la différence entre R et XLSTAT (bien que je ne connaisse pas ce dernier) c'est certainement parce que par défaut XLSTAT donne une valeur des paramètres par rapport au cas où il n'y aurait que l'intercept dans ton modèle. R de son côté se base toujours sur un "témoin" qui est par défaut la première modalité de ton facteur quand tu travailles avec des variables catégorielles. c'est la raison pour laquelle il te fait "disparaitre" Ty-E des results car elle sert par défaut de témoin (c'est la première dans l'ordre alphabétique)
Donc pour tes tests il va falloir t'orienter vers du non-paramétrique ou alors te contenter de faire des boxplots et là tu peux utiliser l'option "notch = T" pour faire apparaitre l'intervalle de confiance de la médiane. Ensuite si les 2 intervalles de confiance se chevauchent c'est que les 2 valeurs ne diffèrent pas significativement l'une de l'autre (généralement au seuil de 5%). Dans le cas contraire les valeurs sont différentes. Une simple approche graphique suffirait d'après cette publication :
Chambers, J. M., Cleveland, W. S., Kleiner, B. and Tukey, P. A.
(1983) _Graphical Methods for Data Analysis._ Wadsworth &
Brooks/Cole.
Donc à toi de voir ensuite ce don tu as besoin...mais il n'est pas toujours nécessaire de tester ce qui est évident : toujours faire des graphiques de ses données :)

...dans ce cas, je suis intéressé pour évaluer l'influence des modalités de Ty (comparaison multiple)...

Là je ne te suis pas trop car j'ai l'impression que c'est ce que tu fais plus bas dans ton post. Donc réexplique différemment ce que tu veux stp.

voilà.

[Logez Maxime]

Bonjour,


Pour répondre a ta question concernant les différences des paramètres, il ne s'agit pas d'un problème de prise en compte de l'intercept mais tout simplement de la modalité qui a servi de référence, si tu regardes les résultats de XLSTAT tu vois Ty-M égale à 0 alors que si tu regardes les résultats de R tu ne vois pas Ty-E ! Ce qui signifie que la variable qui a servi de référence dans les deux cas ne sont pas les mêmes. Il est donc évident que les valeurs des paramètres ne seront pas les mêmes si tu veux comparer les coefficients il te faut modifier les constrats de la variable ou alors modifier l'ordre des levels de Ty :

Code : Tout sélectionner

asFr1 <- rnorm(100)
Ty <- as.factor(sample(c("E","F","FS","M","NI"),100,rep=T))
Ty
[1] E  NI FS M  M  FS NI M  FS M  E  FS F  NI NI NI F  M  F  NI FS M  FS FS F  E  M  F  E  FS F  NI M  M  FS M  F  F  FS FS FS E  FS M  NI FS F
 [48] FS M  F  E  E  E  FS E  NI F  F  E  NI NI NI E  NI M  FS NI M  E  M  M  NI F  M  F  M  E  NI M  M  F  M  E  NI F  E  F  FS FS NI M  NI NI F
 [95] FS NI E  FS F  M
Levels: E F FS M NI

summary(lm(asFr~Ty)) # la référence de Ty sera E (premier level de Ty) donc n'apparait pas dans le summary
Call:
lm(formula = asFr1 ~ Ty)

Residuals:
     Min       1Q   Median       3Q      Max
-1.95058 -0.78142 -0.01678  0.62207  1.92067

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.08257    0.23970   0.344    0.731
TyF         -0.25691    0.32534  -0.790    0.432
TyFS         0.19113    0.31817   0.601    0.549
TyM         -0.42801    0.31214  -1.371    0.174
TyNI         0.12602    0.31817   0.396    0.693

Residual standard error: 0.9588 on 95 degrees of freedom
Multiple R-Squared: 0.06346,    Adjusted R-squared: 0.02402
F-statistic: 1.609 on 4 and 95 DF,  p-value: 0.1783

Ty2 <- relevel(Ty,4) # permet de modifier l'odre des levels en mettant comme premier level le quatrième donc M
summary(lm(asFr1~Ty2))
Call:
lm(formula = asFr1 ~ Ty2)

Residuals:
     Min       1Q   Median       3Q      Max
-1.95058 -0.78142 -0.01678  0.62207  1.92067

Coefficients:
            Estimate Std. Error t value Pr(>|t|) 
(Intercept)  -0.3454     0.1999  -1.728   0.0873 .
Ty2E          0.4280     0.3121   1.371   0.1735 
Ty2F          0.1711     0.2972   0.576   0.5662 
Ty2FS         0.6191     0.2894   2.139   0.0350 *
Ty2NI         0.5540     0.2894   1.914   0.0586 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9588 on 95 degrees of freedom
Multiple R-Squared: 0.06346,    Adjusted R-squared: 0.02402
F-statistic: 1.609 on 4 and 95 DF,  p-value: 0.1783

 # une autre possibilité sans toucher aux facteurs : modifier les contrasts
contrasts(Ty) # montre que E sert de référence (des 0 de partout)
   F FS M NI
E  0  0 0  0
F  1  0 0  0
FS 0  1 0  0
M  0  0 1  0
NI 0  0 0  1
contr.treatment(levels(Ty),4) # M devient la référence
   E F FS NI
E  1 0  0  0
F  0 1  0  0
FS 0 0  1  0
M  0 0  0  0
NI 0 0  0  1
contrasts(Ty) <- contr.treatment(levels(Ty),4)#  la référence est bien M maintenant
summary(lm(asFr1~Ty))# même chose qu'avec Ty2

Call:
lm(formula = asFr1 ~ Ty)

Residuals:
     Min       1Q   Median       3Q      Max
-1.95058 -0.78142 -0.01678  0.62207  1.92067

Coefficients:
            Estimate Std. Error t value Pr(>|t|) 
(Intercept)  -0.3454     0.1999  -1.728   0.0873 .
TyE           0.4280     0.3121   1.371   0.1735 
TyF           0.1711     0.2972   0.576   0.5662 
TyFS          0.6191     0.2894   2.139   0.0350 *
TyNI          0.5540     0.2894   1.914   0.0586 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9588 on 95 degrees of freedom
Multiple R-Squared: 0.06346,    Adjusted R-squared: 0.02402
F-statistic: 1.609 on 4 and 95 DF,  p-value: 0.1783


Avec quoi as-tu testé l'égalité des variances ? En effet si cette condition n'est pas remplie alors oublie les résultats de ce test et passe en non paramétrique avec un test de Kruskal-Wallis.

Maxime

[E.H. [compte supprimé]]

L'égalité des variance a été testée avecBartlett, et la p-value est <0,0001 :(

un probleme vient déjà de la transformation angulaire de mes données, je viens de lire que

la transformation angulaire n’est pas valable lorsque dans les données, p prend les valeurs 0 ou 1. On l’améliorera en remplaçant, avant de prendre des valeurs angulaires, 0 par (1/4n) et 1 par [1-(1/4n)], où n est le nombre d’observations sur la base desquelles est estimé p pour chaque groupe

Or j'ai effectivement un grand nombre de 0.

Il faut donc que je revoie cela mais je suis bloqué sur cette transfo :

Ma population initiale complete (ensemble des données recueillies) comprend N=507 individus.
Dois je réaliser la transformation à ma base complete avec N=507 et ne plus y toucher par la suite ou alors faire au cas par cas, le cas par cas étant que mes analyses se portent successivement selon différents critères avec un nombre d'individus différents?

Par exemple, si j'analyse uniquement les données provenant du Massif "B" (focale$Ma == "B"), N(B) = 325 , dois je appliquer la valeur extrème 1/(4n) avec N=507ou N=325 ?
Ou encore, si j'analyse uniquement les données concernant Espèce "C" (focale$Es == "C"), N(C) =419 , j'applique quelle valeur de N ? 325 ou 507 ?
Sachant que si je choisit l'espèce C, les modalités focale$Ty sont F, FS, NI, E, M, et ont toutes des effectifs différents. Je choisit quoi pour N ?

:( el_boulet :(

[Logez Maxime]

Re,


Que cherches-tu as tester au juste ?


Maxime

[E.H. [compte supprimé]]

je saisi pas ta question je crois.

je cherche a tester l'influence de différentes variables (Ma, Ty, Ag, Li, He, ... en clair, massif, typologie, lieu, tranche horaire, ...) sur le budget temps (des fréquences) accordé à différents comportements (asrF1, F2, F3...) chez des individus

est ca que tu voulais savoir ?

[Logez Maxime]

Bonjour,

je ne suis pas très au fait de cette transformation mais d'après :

la transformation angulaire n’est pas valable lorsque dans les données, p prend les valeurs 0 ou 1. On l’améliorera en remplaçant, avant de prendre des valeurs angulaires, 0 par (1/4n) et 1 par [1-(1/4n)], où n est le nombre d’observations sur la base desquelles est estimé p pour chaque groupe

n doit s'agir de la valeur qui a servi a calculer tes fréquences. Sinon une alternative est de passer sur des tests no paramétriques, qui sont moins puissant mais qui ne sont pas assujettis a des hypothèses de fonctionnement et en plus ça évite de transformer tes données.

Maxime