Why does PROTEST report high significance when PROCRUSTES residuals are high?

[NJOram]

I am trying to test the concurrence of a plant and a microbial community (2 matrices, d1, d2) with Procrustes analysis and protest in vegan. Although the Procrustes residuals are quite high, Protest tells me that the communities are highly significant (with a sum of squares near zero). I am confused why the residuals are high (signaling that the concurrence is low) while the Protest function says otherwise. Can anyone shed light on this issue?

An example:
example_data1.csv
example_data2.csv

# Read Data
# 1-36 sites (rows), 6 species (columns)
d1<-read.csv("example_data1.csv")%>%
  dplyr::select(-X)%>%
  column_to_rownames(var = "ID")%>%
  as.matrix()

# 1-36 sites (rows), 13988 species (columns)
d2<-read.csv("example_data2.csv")%>%
  dplyr::select(-X)%>%
  column_to_rownames(var = "ID")%>%
  as.matrix()

# Transform and run a PCA on each community (each matrix)
# Extract scores 
# PCA 1
mat1<- decostand (d1, 'hell')
pca1<-rda(mat1)
summary(pca1)
plot(pca1)
scores_pca1<-scores(pca1, display = "sites", choices = c(1:6))%>%
  as.data.frame() 

# PCA2 
mat2 <- decostand (d2, 'hell')
pca2<-rda(mat2)
summary(pca2)
plot(pca2)
scores_pca2<-scores(pca2, display = "sites", choices = c(1:6))%>%
  as.data.frame()

###### Procrusts
set.seed(12345)
pro<-procrustes(scores_pca1,scores_pca2)
summary(pro)
plot(pro, kind = 1)
plot(pro, kind = 2)

# Protest to test significance
prot<-protest(x = scores_pca1, y = scores_pca2,permutations = 9999)
prot

Call:
protest(permutations = 9999, x = scores_pca1, y = scores_pca2) 

Procrustes Sum of Squares (m12 squared):        2.22e-16 
Correlation in a symmetric Procrustes rotation:     1 
Significance:  1e-04 

Permutation: free
Number of permutations: 9999

[jarioksa]

For a protets result prot (like in your example), use:

densityplot(permustats(prot))

This will display you the density distribution of randomized values and the line for the observed value. If the observed value is at the tail of distribution then p-value will be low.

You claim that residuals are large, but that is contrary to the results you posted. The last output says that the Procrustes sum of squares is 2.2 × 10^-16 and correlation is 1. NB, protests is based on symmetric Procrustes analysis unlike your prior examples of procrustes().

[NJOram]

Thanks Jari!
If I run
densityplot(permustats(prot))
I get this distribution.. which I think would indicate that something is wrong - why is everything clustered around zero?

I understand that the residuals resulting from protest() are very small (prot in my example). However, the residuals for 'procrustes' are not, right? My goal is to test if significant variation in mat2 (microbial community) is explained by mat1 (the plant community). Is there a better approach to this? Apologies if this is an extremely basic question...

[gavinsimpson]

I see the problem here. The arguments to protest() are not x and y, but are X and Y. When you call protest() as:

prot <- protest(x = scores_pca1, y = scores_pca2, permutations = 9999)

what seems to happen is that scores_pca1 gets used for X and Y and of course that will end up with a procrustes correlation of 1 and sums of squares effectively 0.

I don't know why you are extracting the scores as both procrustes() and protest() will happily work with ordination objects. If you call this as intended you get the right answer:

r$> set.seed(12345)                                                             
r$> prot <- protest(X = scores_pca1, Y = scores_pca2, permutations = 9999)       
r$> prot                                                                        

Call:
protest(X = scores_pca1, Y = scores_pca2, permutations = 9999) 

Procrustes Sum of Squares (m12 squared):        0.8139 
Correlation in a symmetric Procrustes rotation: 0.4314 
Significance:  0.025 

Permutation: free
Number of permutations: 9999

which is what you get if you call this with the ordination objects and pass choices = 1:6 through to scores()`:

r$> set.seed(12345)                                                             
r$> prot <- protest(X = pca1, Y = pca2, permutations = 9999, choices = 1:6)      
r$> prot                                                                        

Call:
protest(X = pca1, Y = pca2, permutations = 9999, choices = 1:6) 

Procrustes Sum of Squares (m12 squared):        0.8139 
Correlation in a symmetric Procrustes rotation: 0.4314 
Significance:  0.025 

Permutation: free
Number of permutations: 9999

@jarioksa perhaps we should add something to protest() to throw a warning/message if X and Y are the same after the scores() calls?

[NJOram]

Thanks so much Gavin - I completely missed the capitals here. The scores were extracted because the PCAs are based on matrices with different numbers of columns. However, the choices = 1:6 is smoother so I will go with this.

Thanks again!