Análise RLQ + 4th Corner

library(ade4)
library(kableExtra)
library(tidyverse)
data(aravo)

Contexto dos dados:

O dataset aravo representa comunidades vegetais em um gradiente de derretimento da neve nos Alpes Franceses. Ele contém:

- matriz L: abundância de espécies (sites x espécies)
- matriz R: variáveis ambientais medidas em cada site
- matriz Q: atributos funcionais medidos em cada espécie

Pergunta ecológica central: “Em que medida a duração do derretimento da neve afeta as comunidades de gramíneas alpinas, a partir da filtragem de atributos funcionais?”

Exploração inicial:

Matriz L: 75 sites x 82 espécies

dim(aravo$spe)
## [1] 75 82
kable(head(aravo$spe[1:10], 10)) %>%
  kable_styling(full_width = F, bootstrap_options = c("striped", "hover", "condensed", "responsive"), font_size = 10)
Agro.rupe Alop.alpi Anth.nipp Heli.sede Aven.vers Care.rosa Care.foet Care.parv Care.rupe Care.semp
AR07 0 0 0 0 0 1 0 0 3 0
AR71 0 0 0 0 0 2 0 1 0 0
AR26 3 0 1 0 1 2 0 0 0 1
AR54 0 0 0 2 0 2 0 0 0 0
AR60 0 0 0 0 0 0 0 0 0 0
AR70 0 0 0 0 0 3 0 1 1 0
AR22 2 0 0 0 2 2 0 1 0 0
AR25 0 0 0 0 0 2 0 1 3 1
AR27 0 0 0 0 1 3 0 1 0 0
AR28 3 0 0 0 2 3 0 0 0 2

Matriz Q: 82 espécies x 8 atributos funcionais

dim(aravo$traits)  
## [1] 82  8
kable(head(aravo$traits, 10)) %>%
  kable_styling(full_width = F, bootstrap_options = c("striped", "hover", "condensed", "responsive"), font_size = 10)
Height Spread Angle Area Thick SLA N_mass Seed
Agro.rupe 6 10 80 60.0 0.12 8.1 218.70 0.08
Alop.alpi 5 20 20 190.9 0.20 15.1 203.85 0.21
Anth.nipp 15 5 50 280.0 0.08 18.0 219.60 0.54
Heli.sede 0 30 80 600.0 0.20 10.6 233.20 1.72
Aven.vers 12 30 60 420.0 0.14 12.5 156.25 1.17
Care.rosa 30 20 80 180.0 0.40 6.5 208.65 1.68
Care.foet 12 40 45 280.0 0.14 19.0 397.10 0.23
Care.parv 25 10 80 520.0 0.20 13.7 237.01 0.61
Care.rupe 9 5 70 75.0 0.20 8.8 181.28 0.64
Care.semp 20 20 50 600.0 0.20 11.5 182.85 0.83

Matriz R: 75 sites x 6 variáveis ambientais (4 contínuas + 2 qualitativas)

dim(aravo$env)
## [1] 75  6
kable(head(aravo$env, 10)) %>%
  kable_styling(full_width = F, bootstrap_options = c("striped", "hover", "condensed", "responsive"), font_size = 10)
Aspect Slope Form PhysD ZoogD Snow
AR07 7 2 1 50 no 140
AR71 1 35 3 40 no 140
AR26 5 0 3 20 no 140
AR54 9 30 3 80 no 140
AR60 9 5 1 80 no 140
AR70 1 30 3 40 no 140
AR22 5 10 3 70 no 160
AR25 6 12 2 15 no 140
AR27 9 2 1 65 no 140
AR28 1 5 3 60 no 150

Etapa 1: Análises preliminares em cada matriz

(1) Matriz L

Usamos uma Correspondence Analysis (CA) A CA é apropriada para dados de abundância ou presença/ausência, preservando a estrutura de composição relativa das comunidades. USa distância qui-quadrado Diferente da PCA, que é linear, a CA lida bem com dados onde espécies têm respostas unimodais ao gradiente

afcL.aravo <- dudi.coa(aravo$spe, scannf = FALSE)
summary(afcL.aravo)
## Class: coa dudi
## Call: dudi.coa(df = aravo$spe, scannf = FALSE)
## 
## Total inertia: 4.214
## 
## Eigenvalues:
##     Ax1     Ax2     Ax3     Ax4     Ax5 
##  0.6607  0.4204  0.3463  0.3037  0.2444 
## 
## Projected inertia (%):
##     Ax1     Ax2     Ax3     Ax4     Ax5 
##  15.677   9.976   8.218   7.205   5.800 
## 
## Cumulative projected inertia (%):
##     Ax1   Ax1:2   Ax1:3   Ax1:4   Ax1:5 
##   15.68   25.65   33.87   41.08   46.88 
## 
## (Only 5 dimensions (out of 74) are shown)

A saída mostra quantos eixos são necessários para resumir a variância (inércia). Aqui, os dois primeiros eixos descrevem gradientes de composição florística. Inertia total: análoga à variância total de uma PCA. OS cinco primeiros eixos explicam ~47% da variância total

(2) Matriz R

Usamos Hill-Smith (dudi.hillsmith), que é uma “PCA para dados mistos” (quantitativos + qualitativos).

É necessário usar os pesos dos sites calculados na CA de L, garantindo compatibilidade entre análises. Cada site em R é ponderado pelo total de abundâncias no sítio em L. A grande diferença entre a RLQ e rodar simplesmente uma PCA com dados ambientais (R) e atributos (Q) está justamente no uso dos pesos da matriz de abundância (L) Se não usar os pesos, a ordenação em R e de Q reflete apenas a variação interna de cada tabela, sem levar em conta o quanto cada site contribui (em espécies) ou o quanto cada espécie realmente aparece (em abundância).

acpR.aravo <- dudi.hillsmith(aravo$env, row.w = afcL.aravo$lw, scannf = FALSE)
summary(acpR.aravo)
## Class: mix dudi
## Call: dudi.hillsmith(df = aravo$env, row.w = afcL.aravo$lw, scannf = FALSE)
## 
## Total inertia: 10
## 
## Eigenvalues:
##     Ax1     Ax2     Ax3     Ax4     Ax5 
##  2.2658  1.8474  1.3064  1.0994  0.9308 
## 
## Projected inertia (%):
##     Ax1     Ax2     Ax3     Ax4     Ax5 
##  22.658  18.474  13.064  10.994   9.308 
## 
## Cumulative projected inertia (%):
##     Ax1   Ax1:2   Ax1:3   Ax1:4   Ax1:5 
##   22.66   41.13   54.20   65.19   74.50 
## 
## (Only 5 dimensions (out of 10) are shown)

Os eixos representam combinações de variáveis (ex: duração do derretimento, altitude, etc.).

(3) Matriz Q Usamos PCA (dados todos contínuos). Se fossem dados fuzzy, poderia utilizar uma dudi.fca. Pesos das espécies vêm da CA de L (col.w).

acpQ.aravo <- dudi.pca(aravo$traits, row.w = afcL.aravo$cw, scannf = FALSE)
summary(acpQ.aravo)
## Class: pca dudi
## Call: dudi.pca(df = aravo$traits, row.w = afcL.aravo$cw, scannf = FALSE)
## 
## Total inertia: 8
## 
## Eigenvalues:
##     Ax1     Ax2     Ax3     Ax4     Ax5 
##  2.4090  1.4983  1.2264  1.0194  0.7469 
## 
## Projected inertia (%):
##     Ax1     Ax2     Ax3     Ax4     Ax5 
##  30.113  18.729  15.330  12.743   9.336 
## 
## Cumulative projected inertia (%):
##     Ax1   Ax1:2   Ax1:3   Ax1:4   Ax1:5 
##   30.11   48.84   64.17   76.91   86.25 
## 
## (Only 5 dimensions (out of 8) are shown)

Mostra quais atributos explicam maior parte da variabilidade entre espécies. Eixo 1 normalmente reflete trade-offs funcionais clássicos (ex: altura vs. SLA).

Etapa 2: RLQ Analysis

A RLQ procura maximizar a covariância entre variáveis ambientais e atributos, ponderada pela presença/abundância das espécies.

rlq.aravo <- rlq(acpR.aravo, afcL.aravo, acpQ.aravo, scannf = FALSE)
plot(rlq.aravo)

Eixo RLQ1: - lado negativo: espécies como Poa supina, Alchemilla pentaphyllea associadas a SLA alto, baixo porte e sementes pequenas a habitats com derretimento tardio. - lado positivo: espécies de ambientes com derretimento precoce, geralmente mais altas e com sementes maiores.

summary(rlq.aravo)
## RLQ analysis
## 
## Class: rlq dudi
## Call: rlq(dudiR = acpR.aravo, dudiL = afcL.aravo, dudiQ = acpQ.aravo, 
##     scannf = FALSE)
## 
## Total inertia: 1.578
## 
## Eigenvalues:
##      Ax1      Ax2      Ax3      Ax4      Ax5 
## 1.367618 0.154783 0.045189 0.006948 0.001993 
## 
## Projected inertia (%):
##     Ax1     Ax2     Ax3     Ax4     Ax5 
## 86.6481  9.8066  2.8630  0.4402  0.1263 
## 
## Cumulative projected inertia (%):
##     Ax1   Ax1:2   Ax1:3   Ax1:4   Ax1:5 
##   86.65   96.45   99.32   99.76   99.88 
## 
## (Only 5 dimensions (out of 8) are shown)
## 
## 
## Eigenvalues decomposition:
##        eig     covar      sdR      sdQ      corr
## 1 1.367618 1.1694520 1.463896 1.530028 0.5221231
## 2 0.154783 0.3934246 1.240738 1.151650 0.2753347
## 
## Inertia & coinertia R (acpR.aravo):
##    inertia      max     ratio
## 1  2.14299 2.265776 0.9458085
## 12 3.68242 4.113141 0.8952819
## 
## Inertia & coinertia Q (acpQ.aravo):
##     inertia      max     ratio
## 1  2.340986 2.409033 0.9717534
## 12 3.667285 3.907346 0.9385616
## 
## Correlation L (afcL.aravo):
##        corr       max     ratio
## 1 0.5221231 0.8128245 0.6423565
## 2 0.2753347 0.6484013 0.4246362

Mostra inércia explicada e correlações entre R e Q.
Valores altos indicam forte correspondência entre atributos e ambiente.

# Visualizações adicionais
par(mfrow = c(1, 3))
s.arrow(rlq.aravo$l1)   # variáveis ambientais projetadas
s.arrow(rlq.aravo$c1)   # atributos funcionais projetados
s.label(rlq.aravo$lQ)   # espécies no espaço de atributos

Teste de significância da RLQ

randtest(rlq.aravo, nrepet = 999, modeltype = 6)
## class: krandtest lightkrandtest 
## Monte-Carlo tests
## Call: randtest.rlq(xtest = rlq.aravo, nrepet = 999, modeltype = 6)
## 
## Number of tests:   2 
## 
## Adjustment method for multiple comparisons:   none 
## Permutation number:   999 
##      Test     Obs  Std.Obs   Alter Pvalue
## 1 Model 2 1.57836 26.26567 greater  0.001
## 2 Model 4 1.57836 12.96874 greater  0.001

Modelos de permutação: - Modelo 2: permuta linhas de L → testa se composição de spp independe de R. - Modelo 4: permuta colunas de L → testa se atributos independem de R. Rejeitar H0 (p < 0.05) → indica que existe associação significativa entre ambiente e atributos funcionais via composição de espécies. - Modelo 6: permuta linha e colunas

Etapa 3: Fourth-corner analysis

Complementa RLQ testando associações bivariadas atributo-ambiente.

four.comb.aravo <- fourthcorner(aravo$env, aravo$spe, aravo$traits,
                                modeltype = 6, nrepet = 999,
                                p.adjust.method.G = "none",
                                p.adjust.method.D = "none")
par(mfrow = c(1, 1))
plot(four.comb.aravo, alpha = 0.05, stat = "D2")

Quadrados coloridos indicam correlações atributo-ambiente: vermelho = correlação positiva, azul = negativa.

Ajuste para múltiplos testes (FDR)

four.comb.aravo.adj <- p.adjust.4thcorner(four.comb.aravo,
                                          p.adjust.method.G = "fdr",
                                          p.adjust.method.D = "fdr")
plot(four.comb.aravo.adj, alpha = 0.05, stat = "D2")

Versão integrada com RLQ (mais intuitiva)

plot(four.comb.aravo.adj, x.rlq = rlq.aravo, alpha = 0.05,
     stat = "D2", type = "biplot")

Etapa 4: RLQ + Fourth-corner

Testa se a associação R-Q está concentrada em eixos específicos

testQaxes.comb.aravo <- fourthcorner.rlq(rlq.aravo, modeltype = 6,
                                         typetest = "Q.axes", nrepet = 999,
                                         p.adjust.method.G = "fdr",
                                         p.adjust.method.D = "fdr")

testRaxes.comb.aravo <- fourthcorner.rlq(rlq.aravo, modeltype = 6,
                                         typetest = "R.axes", nrepet = 999,
                                         p.adjust.method.G = "fdr",
                                         p.adjust.method.D = "fdr")
print(testQaxes.comb.aravo, stat = "D")
## Fourth-corner Statistics
## ------------------------
## Permutation method  Comb. 2 and 4  ( 999  permutations)
## 
## Adjustment method for multiple comparisons:   fdr 
## call:  fourthcorner.rlq(xtest = rlq.aravo, nrepet = 999, modeltype = 6,      typetest = "Q.axes", p.adjust.method.G = "fdr", p.adjust.method.D = "fdr") 
## 
## ---
## 
##              Test Stat           Obs    Std.Obs     Alter Pvalue
## 1  AxcR1 / Height    r  2.326199e-01  2.3113906 two-sided  0.022
## 2  AxcR2 / Height    r  3.058482e-05 -0.0411559 two-sided  0.978
## 3  AxcR1 / Spread    r -8.340492e-02 -0.9252242 two-sided  0.373
## 4  AxcR2 / Spread    r  2.220989e-02  0.6499801 two-sided  0.514
## 5   AxcR1 / Angle    r  2.925318e-01  2.8948155 two-sided  0.003
## 6   AxcR2 / Angle    r  1.618969e-01  4.1673943 two-sided  0.001
## 7    AxcR1 / Area    r -6.873643e-02 -0.7083137 two-sided  0.472
## 8    AxcR2 / Area    r -2.246150e-01 -5.7424004 two-sided  0.001
## 9   AxcR1 / Thick    r  1.846911e-01  1.8407105 two-sided  0.067
## 10  AxcR2 / Thick    r  7.630935e-02  2.0008325 two-sided  0.043
## 11    AxcR1 / SLA    r -4.983705e-01 -5.1427124 two-sided  0.001
## 12    AxcR2 / SLA    r  3.518139e-02  0.5218505 two-sided  0.594
## 13 AxcR1 / N_mass    r -4.244131e-01 -4.4875164 two-sided  0.001
## 14 AxcR2 / N_mass    r  1.153893e-01  2.0243439 two-sided   0.05
## 15   AxcR1 / Seed    r  1.555753e-01  1.5145685 two-sided  0.132
## 16   AxcR2 / Seed    r -5.490263e-02 -1.3793518 two-sided  0.161
##            Pvalue.adj   
## 1  0.0502857142857143  .
## 2               0.978   
## 3   0.459076923076923   
## 4   0.548266666666667   
## 5               0.008 **
## 6               0.004 **
## 7   0.539428571428571   
## 8               0.004 **
## 9   0.119111111111111   
## 10              0.086  .
## 11              0.004 **
## 12             0.6336   
## 13              0.004 **
## 14 0.0666666666666667  .
## 15             0.2112   
## 16  0.234181818181818   
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(testRaxes.comb.aravo, stat = "D")
## Fourth-corner Statistics
## ------------------------
## Permutation method  Comb. 2 and 4  ( 999  permutations)
## 
## Adjustment method for multiple comparisons:   fdr 
## call:  fourthcorner.rlq(xtest = rlq.aravo, nrepet = 999, modeltype = 6,      typetest = "R.axes", p.adjust.method.G = "fdr", p.adjust.method.D = "fdr") 
## 
## ---
## 
##                  Test   Stat         Obs      Std.Obs     Alter Pvalue
## 1      Aspect / AxcQ1      r -0.01202699 -0.233137230 two-sided  0.798
## 2       Slope / AxcQ1      r  0.23055049  4.503019388 two-sided  0.001
## 3      Form.1 / AxcQ1 Homog.  0.23632349 -0.162744183      less   0.46
## 4      Form.2 / AxcQ1 Homog.  0.08451815 -0.420350341      less  0.357
## 5      Form.3 / AxcQ1 Homog.  0.26620328 -1.236514903      less  0.122
## 6      Form.4 / AxcQ1 Homog.  0.10425439 -2.883182283      less  0.003
## 7      Form.5 / AxcQ1 Homog.  0.13081695 -0.646414069      less  0.296
## 8       PhysD / AxcQ1      r  0.30328900  4.772491319 two-sided  0.001
## 9    ZoogD.no / AxcQ1 Homog.  0.58253649  5.382923466      less      1
## 10 ZoogD.some / AxcQ1 Homog.  0.29769156 -1.844780101      less   0.03
## 11 ZoogD.high / AxcQ1 Homog.  0.11941619 -1.767750017      less  0.029
## 12       Snow / AxcQ1      r -0.51057226 -4.879101841 two-sided  0.001
## 13     Aspect / AxcQ2      r -0.09633526 -2.240878848 two-sided  0.028
## 14      Slope / AxcQ2      r -0.03663561 -0.733050206 two-sided  0.482
## 15     Form.1 / AxcQ2 Homog.  0.17514702 -1.506139244      less   0.05
## 16     Form.2 / AxcQ2 Homog.  0.09144593 -0.007704137      less  0.578
## 17     Form.3 / AxcQ2 Homog.  0.28926185 -0.073269925      less  0.501
## 18     Form.4 / AxcQ2 Homog.  0.18092965  0.906902190      less  0.812
## 19     Form.5 / AxcQ2 Homog.  0.24091276  1.822646546      less   0.95
## 20      PhysD / AxcQ2      r  0.10215843  1.598826091 two-sided   0.11
## 21   ZoogD.no / AxcQ2 Homog.  0.30385849 -2.823094038      less  0.002
## 22 ZoogD.some / AxcQ2 Homog.  0.38464080  0.356147625      less  0.646
## 23 ZoogD.high / AxcQ2 Homog.  0.24392686  2.142584093      less  0.972
## 24       Snow / AxcQ2      r  0.07594269  0.702883457 two-sided   0.49
##           Pvalue.adj   
## 1  0.885818181818182   
## 2              0.008 **
## 3  0.653333333333333   
## 4             0.5535   
## 5  0.225230769230769   
## 6             0.0144  *
## 7             0.4736   
## 8              0.008 **
## 9                  1   
## 10              0.08  .
## 11              0.08  .
## 12             0.008 **
## 13            0.0672  .
## 14 0.653333333333333   
## 15 0.109090909090909   
## 16 0.730105263157895   
## 17 0.707294117647059   
## 18 0.885818181818182   
## 19             0.996   
## 20              0.22   
## 21             0.012  *
## 22             0.816   
## 23                 1   
## 24 0.653333333333333   
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Visualizações em tabela

par(mfrow = c(1, 2))
plot(testQaxes.comb.aravo, alpha = 0.05, type = "table", stat = "D2")
plot(testRaxes.comb.aravo, alpha = 0.05, type = "table", stat = "D2")

Visualizações em biplot

par(mfrow = c(1, 2))
plot(testQaxes.comb.aravo, alpha = 0.05, type = "biplot", stat = "D2",
     col = c("black", "blue", "orange", "green"))
plot(testRaxes.comb.aravo, alpha = 0.05, type = "biplot", stat = "D2",
     col = c("black", "blue", "orange", "green"))

Conclusão geral

A RLQ mostrou que há forte relação entre variáveis ambientais, (principalmente duração do derretimento da neve) e atributos funcionais das espécies alpinas.
A Fourth-corner identificou quais atributos específicos (ex: SLA alto, baixa altura, sementes pequenas) estão mais associados a ambientes de derretimento tardio.
Isso sugere um processo de filtragem ambiental claro: a duração do manto de neve atua como filtro, selecionando espécies com determinado conjunto de atributos funcionais.

Informações importantes

RLQ isolado

Vai mostrar eixos ordinais com gradientes ambientais (ex.: Temp/Flow) e atributos de espécies associados.
É exploratório → não diz se as correlações são estatisticamente significativas.
Exploração inicial, para ver gradientes gerais

4th-corner isolado

Vai imprimir uma tabela com correlações e p-valores entre cada variável ambiental e cada atributo.
Pode dar muitos testes (aqui: 3 × 3 = 9), e alguns podem ser falsos positivos.
Quando há hipóteses específicas de pares ambiente-atributo.

RLQ + 4th-corner

Testa as associações ao longo dos eixos principais do RLQ, não par a par.
Dá resultados mais robustos e uma síntese visual mais clara dos pares que realmente importam.
Integrar exploração + inferência. Reduz ruído, evita inflar testes múltiplos, e dá uma síntese ecológica mais clara.

Referências

Dolédec, S., Chessel, D., ter Braak, C.J.F. & Champely, S. (1996). Matching species traits to environmental variables: a new three-table ordination method. Ecology, 77(4), 1227–1236.

Dray, S., & Dufour, A.B. (2007). The ade4 package: implementing the duality diagram for ecologists. Journal of Statistical Software, 22(4), 1–20.

Dray, S., Choler, P., Dolédec, S., et al. (2014). Combining the fourth‐corner and the RLQ methods for assessing trait responses to environmental variation. Ecology, 95(1), 14–21.