R: Curvas de Respuesta Principal (PRC) 3. Significatividad de los ejes

Se nos quedó en el tintero en nuestros análisis PRC un tema algo ecabroso y complicado que es el análisis de la significatividad de los ejes. El test de Monte Carlo para el primer eje no tiene problema en vegan con la función anova.cca:

## Significación primer eje

anova(mod, strata = week, perm.max = 499, first = T)

## Permutation test for rda under reduced model
## Permutations stratified within 'week'
## 
## Model: prc(response = pyrifos, treatment = dose, time = week)
##          Df   Var    F N.Perm Pr(>F)   
## RDA1      1  25.3 15.1    199  0.005 **
## Residual 77 129.0                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

En la última versión de vegan (2.0-5) se introdujo la modificación adecuada en esta función para que se pudiera utilizar el parámetro by con la opción “axis” con un objeto proveniente de de la función prc. Supuestamente con este parámetro la función realiza los test de Monte Carlo para todos los ejes extrayendo los scores de cada eje y utilizándolos como covariables para el siguiente en cada paso, tal como debe hacerse para garantizar la independencia de cada uno de ellos:

anova(mod, strata = week, perm.max = 499, by = "axis")

## Model: rda(formula = pyrifos ~ dose * week + Condition(week))
##          Df   Var     F N.Perm Pr(>F)   
## RDA1      1  25.3 15.10    199  0.005 **
## RDA2      1   8.3  4.95    199  0.005 **
## RDA3      1   6.0  3.61    199  0.005 **
## RDA4      1   4.8  2.85    199  0.005 **
## RDA5      1   4.1  2.48    199  0.005 **
## RDA6      1   3.9  2.30    199  0.005 **
## RDA7      1   3.6  2.14    199  0.010 **
## RDA8      1   3.3  1.99    199  0.005 **
## RDA9      1   3.1  1.84    299  0.017 * 
## RDA10     1   2.6  1.52    499  0.048 * 
## RDA11     1   2.5  1.47    499  0.066 . 
## RDA12     1   2.2  1.32    199  0.115   
## RDA13     1   2.1  1.27     99  0.130   
## RDA14     1   1.9  1.16     99  0.340   
## RDA15     1   1.8  1.07     99  0.430   
## RDA16     1   1.6  0.97     99  0.570   
## RDA17     1   1.6  0.94     99  0.590   
## RDA18     1   1.4  0.86     99  0.730   
## RDA19     1   1.4  0.83     99  0.700   
## RDA20     1   1.3  0.77     99  0.810   
## RDA21     1   1.2  0.72     99  0.870   
## RDA22     1   1.1  0.68     99  0.940   
## RDA23     1   1.0  0.60     99  0.980   
## RDA24     1   0.9  0.55     99  1.000   
## RDA25     1   0.9  0.51     99  1.000   
## RDA26     1   0.8  0.47     99  0.980   
## RDA27     1   0.7  0.45     99  1.000   
## RDA28     1   0.7  0.43     99  1.000   
## RDA29     1   0.7  0.41     99  1.000   
## RDA30     1   0.6  0.36     99  1.000   
## RDA31     1   0.6  0.35     99  1.000   
## RDA32     1   0.5  0.32     99  1.000   
## RDA33     1   0.5  0.31     99  1.000   
## RDA34     1   0.4  0.26     99  1.000   
## RDA35     1   0.4  0.25     99  1.000   
## RDA36     1   0.4  0.24     99  1.000   
## RDA37     1   0.4  0.22     99  1.000   
## RDA38     1   0.3  0.20     99  1.000   
## RDA39     1   0.3  0.20     99  1.000   
## RDA40     1   0.3  0.18     99  1.000   
## RDA41     1   0.3  0.17     99  1.000   
## RDA42     1   0.3  0.16     99  1.000   
## RDA43     1   0.2  0.12     99  1.000   
## RDA44     1   0.2  0.11     99  1.000   
## Residual 77 129.0                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sin embargo, ya adelanto que los resultados de estos análisis no son correctos. Esto lo he comprobado con resultados obtenidos con CANOCO de otra serie de datos propios. El caso es que la función en vegan parece no tomar las puntuaciones adecuadas de cada eje. He dejado esta cuestión en la página de del foro de vegan por si se tratara de un bug.

A la espera de conocer qué está pasando realmente con el parámetro by de la función anova.cca, tenemos una solución alternativa para comprobar la significatividad de los ejes de nuestra PRC subsiguientes al primero. Lo primero es saber que una PRC no es más que una RDA de la forma rda (species ~ Treatment * Year + Condition (Year)). Lo segundo que necesitamos saber es dónde están las puntuaciones de los puntos de muestreo que son combinaciones lineales de las variables ambientales y que son las que tenemos que utilizar como covariables en el análisis del siguiente eje (las SamE de CANOCO). El resultado de la PRC es un objeto cca.object y dichas puntuaciones se encuentran en cca.object$CCA$u. Así que solo tenemos que extraerlas y realizar un RDA añadiéndolas como covariable (argumento Condition). Por último, realizaremos el test de Monte Carlo para el primer eje de este análisis que correspondería a nuestro segundo eje de PRC.

## Significación segundo eje

SamE <- mod$CCA$u
SamE <- SamE[, 1]
secondaxis <- rda(pyrifos ~ dose * week + Condition(week) + Condition(SamE))
anova(secondaxis, strata = week, perm.max = 499, first = T)

## Permutation test for rda under reduced model
## Permutations stratified within 'week'
## 
## Model: rda(formula = pyrifos ~ dose * week + Condition(week) + Condition(SamE))
##          Df   Var    F N.Perm Pr(>F)
## RDA1      1   8.3 4.95     99    0.7
## Residual 77 129.0

He comprobado, ya dije, esto con datos propios porque carezco de los resultados del artículo de Ter Braak y van den Brink en el que están basados datos de pyrifos (solo aparecen en él las p y el porcentaje de varianza recogida por los ejes), y tanto los resultados de las F en vegan y CANOCO son similares y los eigenvalues de la prc y de la posterior rda coinciden.

Solo queda saber qué pasa con by=“axis”, porque sería bastante más cómodo hacerlo así.

Por último, como habréis comprobado el formato de las operaciones en R de este post se ve algo más vistoso. Lo he conseguido utilzando el paquete knitr y la interface gráfica RStudio. knitr le permite a uno elaborar mediante markdown informes en HTML a partir de operaciones con R, como véis, resaltando de forma diferente las órdenes en la línea de comandos y los resultados. Permite además incrustar las gráficas en el mismo documento, pero no lo tengo aún muy dominado. Por el momento, parece que el código HTML que genera es aceptado bastante bien por Blogger