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
Remarques : beaucoup de double zéros et d'espèces rares
Code : Sélectionner tout - Visualiser dans une fenêtre à part
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
Un tableau avec les données environnementales de 13 paramètres physico-chimiques dans les mêmes 30 sites
Remarques : pas la même amplitude et pas les mêmes unités
Code : Sélectionner tout - Visualiser dans une fenêtre à part
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
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 :
Les questions
Code : Sélectionner tout - Visualiser dans une fenêtre à part
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
Pourquoi il me sélectionne le modèle avec les intercept alors qu'il existe un meilleur modèle ?
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 ?
Code : Sélectionner tout - Visualiser dans une fenêtre à part
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
Merci d'avance !
Partager