Why is constrained RDA symmetrical? I am having a very weird behaviour when performing constrained RDA. The hypothesis I am testing is: is there a difference between treated and control samples? The experimental design comprise 10 sites, where a control and a treatment were sampled.
As far as I understand, this is the classical example of random effects caused by the sampling site so my model would be:
data_agglomerated_at_taxa_level ~ Treatment + Condition(Site)
However, something very interesting happens: as you can see from the picture:

what I am observing are samples that are symmetrical (I and III quadrants, II and IV quadrants). Each couple of symmetrical data point is the couple $(control_{s}, treatment_{s})$, $s = site$.
If I remove the constraint on the Site, then the data points are not symmetrical anymore so, as far as I can tell, there is something happening with the Site. So said, I can't really understand what that might be. Any ideas?
PS: i am performing the analysis with Phyloseq ordinate() function and the issue seems to be related to the type of data in the sample_data(phylo_object).
UPDATE:
i am adding a reproducible example, and another issue that i found:
library("vegan")

# create fake otu
otu_fake <- matrix(0, nrow=10, ncol=100)

# zero inflated random fake otu data
for (i in 1:100) {
    zero_infl <- ifelse(rbinom(n = 10, size = 1, prob = 0.90)==1, 0, 1)
    ones <- which(zero_infl==1)
    otu_fake[ones, i] <- sample.int(10, length(ones))
}
colnames(otu_fake) <- paste0("otu_", 1:100)
rownames(otu_fake) <- paste0("sample_", 1:10)

# create fake meta
fake_meta <- data.frame(sample=paste0("sample_", 1:10), treatment=c(rep("CTR", 5), rep("TRT", 5)), site=c(1:5, 1:5))

# perform RDA
rda_veg <- rda(formula = otu_fake ~ treatment+Condition(site), data = fake_meta)

# plot it
dev.new()
plot(rda_veg, main="Site as integer")

# cast site to factors
fake_meta$site <- as.factor(fake_meta$site)

# perform RDA
rda_veg <- rda(formula = otu_fake ~ treatment+Condition(site), data = fake_meta)

# plot it
dev.new()
plot(rda_veg, main="Site as factors")

This is an example of the experiment i am working on. The funny thing is that the symmetry appears only after casting the site variable to factor (or character).
So, updated question is: is this due to the design/mathematical structure of the data, is there something else, or am I missing something?
UPDATE on the update
I provide here more code where the symmetry issue can be tested, while the factor/integer problem remains.
To see the symmetry, samples needs to equal to 10. For higher values, the symmetry is lost.
library("vegan")

# choose a number of taxa and samples
taxas <- 100
samples <- 30

# with samples = 10 the setting is:
# 1 site comprise 1 control and 1 treated sample

# with samples = 20, 30, 40... the setting becomes:
# 1 site comprise 2 controls and 2 treated samples

# create fake otu
otu_fake <- matrix(0, nrow=samples, ncol=taxas)

# zero inflated random fake otu data
for (i in 1:taxas) {
    zero_infl <- ifelse(rbinom(n = samples, size = 1, prob = 0.90)==1, 0, 1)
    ones <- which(zero_infl==1)
    otu_fake[ones, i] <- sample.int(samples, length(ones))
}
colnames(otu_fake) <- paste0("otu_", 1:taxas)
rownames(otu_fake) <- paste0("sample_", 1:samples)

# create fake meta
fake_meta <- data.frame(sample=paste0("sample_", 1:samples), treatment=c(rep("CTR", samples/2), rep("TRT", samples/2)), site=1:5)

# perform RDA
rda_veg <- rda(formula = otu_fake ~ treatment+Condition(site), data = fake_meta)

# plot it
dev.new()
plot(rda_veg, main=paste("Site as integers with ", samples, " samples"))

# cast site to factors
fake_meta$site <- as.factor(fake_meta$site)

# perform RDA
rda_veg <- rda(formula = otu_fake ~ treatment+Condition(site), data = fake_meta)

# plot it
dev.new()
plot(rda_veg, main=paste("Site as factors with ", samples, " samples"))

 A: This really will happen always with this specific design – or if I understood your design correctly: you did not have a reproducible example. I understood that for each "Site" you have only two units, one Control, another Treatment. You have only one constrained axis and the second axis is unconstrained PC. So the first axis will separate Control and Treatment. The effect of Condition(Site) is to make the mean of sites equal which in this multivariate setting means that their mean will be at the origin of the constrained axis (RDA1). Second axis, PC1, is unconstrained, and it shows how the observed values of your SUs differ from that mean, and Control and Treatment are equally far away from their mean or the origin. If you draw a line from Control and Treatment point of a Site, the line must go through their mean – or the origin - and distance from the origin is equal for Control and Treatment because this is how mean (average) is defined for two observations.
Was so long that I try to summarize: by the definition of partial RDA, the mean of each Site is at the RDA1 axis value 0, and by the definition of mean, the distance to origin is equal to the mean for both points of a Site.
I think that it may make sense to show only the RDA1 axis.
EDIT: To clarify, the reason for symmetry is that you Condition on mean of pairs of points. That mean will be at the origin (0,0), and the two points will be at equal distances and opposite directions from the origin. This concerns in particular unconstrained axes, or residual PC's.
You can produce this with any random seed using:
x <- matrix(rnorm(10*15), nrow=10) # random data
cl <- gl(5,2) # class for pairs of rows
library(vegan)
pcpair <- rda(x ~ Condition(cl)) # center for pairs of points
plot(pcpair, dis="si")
ordisegments(pcpair, cl) # connect two points of the class 

A: I don't know much about RDA and CCA, but this might help you into the right direction.
The function seems to apply some transformation before performing PCA. This transformation seems to center the data in a way such that the average PC score per site is zero.
### rda uses Ybar, some centered and re-scaled values per column
### the object contains 'rda_veg[6]$pCCA$Fit' this is a vector
### that resembles the center of the pca scores per site
### when you subtract it then you get the symmetry

### perform pca with this
otu_transform = rda_veg[3]$Ybar-rda_veg[6]$pCCA$Fit


pc = prcomp(otu_transform, center = FALSE)$rotation

### plot pc1 and pc2
### still not exactly the same but similar
sites1 = otu_transform[,] %*% pc[,1]
sites2 = otu_transform[,] %*% pc[,2] 
plot(sites1,sites2, col = fake_meta$site)


The upper graph pc1 vs pc2 is made with your code but adding a choices parameter plot(rda_veg, choices = c(2,3), main="Site as factors") 
