Fonction [vegan] : ordiR2step
Bonjour à tous,
je voudrais faire une cca mais je n'arrive pas à utiliser la fonction ordiR2step pour trouver le meilleur modèle.
Les données
Un tableau avec les données d'abondance de 48 espèces de diatomées dans 30 sites
Code:
1
2
3
4
5
6
7
| AC001A AC013A AC013E AM011A AM012A AS001A AU002A AU003B CC001A CC002A CC9997 CM004A CO001A CY002A CY003A CY009A CY011A
4 0.00 0.55 0 0.74 0.92 1.66 4.60 0 0.00 0.00 0.00 0 0.18 1.11 0.00 0 0.00
7 0.36 3.40 0 1.07 8.05 0.36 0.00 0 2.15 3.40 0.00 0 3.94 1.97 3.04 0 0.72
31 0.90 1.08 0 0.90 5.39 0.00 0.00 0 0.00 0.18 0.18 0 0.72 0.00 0.00 0 0.00
34 0.17 0.52 0 0.69 0.35 0.00 0.00 0 9.69 7.96 15.23 0 0.52 3.46 2.77 0 9.17
37 0.00 6.84 0 2.54 2.34 0.19 0.00 0 0.00 0.00 0.00 0 2.15 0.59 0.19 0 0.00
42 0.18 0.91 0 0.00 0.73 0.36 14.03 0 0.00 0.00 0.00 0 0.91 0.00 0.00 0 0.00 |
Remarques : beaucoup de double zéros et d'espèces rares
Un tableau avec les données environnementales de 13 paramètres physico-chimiques dans les mêmes 30 sites
Code:
1
2
3
4
5
6
7
8
| > head(diatomEN)
Depth pH Cond Na K Mg Ca Cl SO4 Alk TP NO3 Chla
4 3.0 8.013077 481.1915 386.00 79.25 462.00 3565.75 618.50 1635.50 3457.692 77.93539 926.1538 33.04
7 1.2 8.001538 569.1262 1244.50 243.25 603.25 5353.00 1641.75 1946.75 3513.846 167.03847 1569.2308 198.08
120 4.0 8.590769 527.7184 1028.00 103.25 192.50 4478.00 957.75 940.25 4198.461 475.53308 870.7692 64.08
34 1.0 7.829231 442.4185 1239.25 145.00 303.25 2821.75 1129.00 979.75 2236.923 460.62692 1416.1538 100.96
31 1.2 8.156154 502.6900 619.00 241.75 387.75 4472.00 626.25 1211.50 3995.385 635.25153 933.0769 272.00
42 1.7 8.245455 275.2609 387.25 44.75 140.00 2159.00 487.00 361.75 1908.182 67.96636 934.5455 153.44 |
Remarques : pas la même amplitude et pas les mêmes unités
Ces données sont issues d'une étude de Bennion (1994),*DOI:*10.1007/BF00026729
La CCA et la sélection
Voila mon script et mes sorties :
Code:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
| > cca.all <- cca (downweight(diatomAB) ~ ., data=diatomEN)
> cca.all
Call: cca(formula = downweight(diatomAB) ~ Depth + pH + Cond + Na + K + Mg + Ca + Cl + SO4 + Alk + TP + NO3 + Chla, data = diatomEN)
Inertia Proportion Rank
Total 4.9958 1.0000
Constrained 2.2217 0.4447 13
Unconstrained 2.7740 0.5553 16
Inertia is scaled Chi-square
Eigenvalues for constrained axes:
CCA1 CCA2 CCA3 CCA4 CCA5 CCA6 CCA7 CCA8 CCA9 CCA10 CCA11 CCA12 CCA13
0.3687 0.2978 0.2694 0.2659 0.2061 0.1832 0.1554 0.1326 0.0936 0.0818 0.0752 0.0603 0.0318
Eigenvalues for unconstrained axes:
CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8 CA9 CA10 CA11 CA12 CA13 CA14 CA15 CA16
0.4642 0.3707 0.3403 0.2514 0.2339 0.1933 0.1794 0.1677 0.1514 0.1034 0.0813 0.0646 0.0606 0.0569 0.0323 0.0226
> anova (cca.all, permutations = 99999)
Permutation test for cca under reduced model
Permutation: free
Number of permutations: 99999
Model: cca(formula = downweight(diatomAB) ~ Depth + pH + Cond + Na + K + Mg + Ca + Cl + SO4 + Alk + TP + NO3 + Chla, data = diatomEN)
Df ChiSquare F Pr(>F)
Model 13 2.2217 0.9857 0.5354
Residual 16 2.7740
> adjRsq.cca <- RsquareAdj (cca.all)$adj.r.square # adjR2 as stopping criteria used in ordiR2step
> cca.0 <- cca (diatomAB ~ 1, data = diatomEN) # empty model only with intercept
> sel.osR2.cca <- ordiR2step (cca.0, scope = formula (cca.all), direction = 'forward', R2scope = adjRsq.cca, permutations = 99, trace = F)
> sel.osR2.cca
Call: cca(formula = downweight(diatomAB) ~ 1, data = diatomEN)
Inertia Rank
Total 5.812
Unconstrained 5.812 29
Inertia is scaled Chi-square
Eigenvalues for unconstrained axes:
CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
0.6780 0.5733 0.4857 0.4019 0.3967 0.3862 0.3629 0.3124
(Showing 8 of 29 unconstrained eigenvalues)> sel.osR2.cca_adj <- sel.osR2.cca
> sel.osR2.cca_adj$anova$`Pr(>F)` <- p.adjust (sel.osR2.cca$anova$`Pr(>F)`, method = 'holm', n = ncol (diatomEN))
> sel.osR2.cca_adj$anova
$`Pr(>F)`
[1] NA |
Les questions
Pourquoi il me sélectionne le modèle avec les intercept alors qu'il existe un meilleur modèle ?
Code:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
| > cca1 <- cca (downweight(diatomAB) ~ K + Cl + Chla , data=diatomEN)
> cca1
Call: cca(formula = downweight(diatomAB) ~ K + Cl + Chla, data = diatomEN)
Inertia Proportion Rank
Total 4.9958 1.0000
Constrained 0.7398 0.1481 3
Unconstrained 4.2560 0.8519 26
Inertia is scaled Chi-square
Eigenvalues for constrained axes:
CCA1 CCA2 CCA3
0.28750 0.23998 0.21234
Eigenvalues for unconstrained axes:
CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
0.5789 0.4712 0.4147 0.3274 0.2976 0.2764 0.2333 0.2151
(Showing 8 of 26 unconstrained eigenvalues)
> anova (cca1, permutations = 99999)
Permutation test for cca under reduced model
Permutation: free
Number of permutations: 99999
Model: cca(formula = downweight(diatomAB) ~ K + Cl + Chla, data = diatomEN)
Df ChiSquare F Pr(>F)
Model 3 0.7398 1.5065 0.03005 *
Residual 26 4.2560
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 |
Autre question, j'ai entendu parlé d'une fonction qui utilise le R2 (non ajusté il me semble) ET le critère AIC pour sélectionner le meilleur modèle. Connaissez-vous cette fonction ?
Merci d'avance !
