Fonction 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:

  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:

> 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:

> 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:

> 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 !