2
$\begingroup$

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:

enter image description here

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"))
$\endgroup$
5
  • $\begingroup$ What is on the z-axis the RDA component? Why is symmetry weird? $\endgroup$ Nov 28, 2022 at 16:07
  • $\begingroup$ this I did not verified...also I found it more unexpected than weird. I have never seen anything like that (but I am not an expert). $\endgroup$
    – gabt
    Nov 28, 2022 at 16:26
  • $\begingroup$ The effect also occurs when you plot pc1 versus pc2. The rda function seems to do some sort of centering and scaling that makes all the treatments opposite of each other for each site. $\endgroup$ Nov 29, 2022 at 9:48
  • $\begingroup$ @SextusEmpiricus may I ask you how do you plot pc1 versus pc2? $\endgroup$
    – gabt
    Nov 29, 2022 at 9:50
  • $\begingroup$ plot(rda_veg, choices = c(2,3), main="Site as factors") use the choices parameter $\endgroup$ Nov 29, 2022 at 10:30

2 Answers 2

3
$\begingroup$

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 
$\endgroup$
4
  • $\begingroup$ I see the point: the problem is due to the experimental design which comprise only one control and one treatment for one site. in contrast, what I am observing would not happen for a design where I have, for example, 3 controls and 3 treated for each site. is this correct? (i will also try to provide a working example, but not with the original data which are too big). $\endgroup$
    – gabt
    Nov 28, 2022 at 16:26
  • $\begingroup$ When you parallel out "Site", points of each site will be centred to the origin. You see that easily for two-point Sites, but the situation is similar when you have more points in a Site (class). This centring will appear on all unconstrained PC-axes that are all centred to class-means (Site means). CCA has the same feature, but not as visible because points are weighted, and the weights influence their distances from the origin. $\endgroup$ Nov 28, 2022 at 21:15
  • $\begingroup$ i am working on my understanding of what you said. however I found another issue which is related to the data type which encodes the variable "site" and I am quite puzzled, now. $\endgroup$
    – gabt
    Nov 29, 2022 at 8:34
  • $\begingroup$ I think I got the point. Also, as far as i can see from the help of cca, the variable should be a factor. however if I cast "Site" from integer to factor, the ordination goes sideways. Is there a reason for that? Same applies to the example you provided. If I define cl <- rep(1:5, each=2) then run the RDA, the plot will be very different. A priori, how do I know which one is correct? Maybe this should be an entirely different question, though. $\endgroup$
    – gabt
    Nov 29, 2022 at 13:44
1
$\begingroup$

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)

PCA comparison

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")

$\endgroup$
1
  • $\begingroup$ thank you for the example, i am working on the mathematics! $\endgroup$
    – gabt
    Dec 2, 2022 at 7:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.