First, I did DB-RDA for a community dataset with two factors and one covariate. The results showed the interaction effect was significant. Then, I need to do the pairwise comparison for all the levels. I posted my R code here and my question is how to modify the pairwise.adonis function below to let it allow my covariate (claysilt). Thanks much!

all_id <- obj2_1[, c("cnty", "site_id", "rep", "System", "Depth", "Year","vxd","vxy","dxy","vxdxy", "claysilt")]
all_fame_m <- obj2_1[, c("m027","m030","m035","m036","m037","m038",
famem <- cbind.data.frame(all_id, all_fame_m)

obj2_1_dbRDA <- capscale(all_fame_m~System*Year+Condition(claysilt), distance = "bray", data=obj2_1)
anova(obj2_1_dbRDA, by="term", perm=4999)

pairwise.adonis <- function(x,factors, sim.method, p.adjust.m)
  co = as.matrix(combn(unique(factors),2))
  pairs = c()
  F.Model =c()
  R2 = c()
  p.value = c()

  for(elem in 1:ncol(co)){
    ad = adonis(x[factors %in% c(as.character(co[1,elem]),as.character(co[2,elem])),] ~ 
                  factors[factors %in% c(as.character(co[1,elem]),as.character(co[2,elem]))] , 
                permutations = 9999, method =sim.method);
    pairs = c(pairs,paste(co[1,elem],'vs',co[2,elem]));
    F.Model =c(F.Model,ad$aov.tab[1,4]);
    R2 = c(R2,ad$aov.tab[1,5]);
    p.value = c(p.value,ad$aov.tab[1,6])
  p.adjusted = p.adjust(p.value,method=p.adjust.m)
  pairw.res = data.frame(pairs,F.Model,R2,p.value,p.adjusted)
PW.Adonis1=pairwise.adonis(all_fame_m,famem$vxy,sim.method="bray",p.adjust.m = "bonferroni")
### vxy is a variable I combined System * Year, so it has all the levels.
  • $\begingroup$ Just to clarify, is what you're wanting a way to test for the effect due to covariate? $\endgroup$ – Todd Burus Feb 21 '20 at 4:37
  • $\begingroup$ I have factor A (two levels) and B (4 levels) and a covariate C (continuous variable). Actually I'm not very interested in the covariate. Including it in the db-rda was because it influences my community composition, which was from the experiment design we cannot control. LIke the covariate in the ANCOVA. Here my data was community species composition. what I want to do is to compare all levels of my two-way interactions (A*B). $\endgroup$ – sunshine Feb 21 '20 at 4:41
  • $\begingroup$ You should be able to test what you want by fitting a full model with covariate in R and looking at the Type III ANOVA. You'll need to have a list of contrasts in there as well. $\endgroup$ – Todd Burus Feb 21 '20 at 5:11
  • $\begingroup$ I know how to do it if I have one single dependent variable but now I have a multivariate species data. As I mentioned in the question above, there is an R function "multiconstrained" can do the multiple comparisons for all level but this function doesn't allow me to add the covariate. My statistical knowledge and R coding skill only allow me to use the R package available online. I don't know how to write my own code to do so. Do you have more specific codes? Thanks a lot! $\endgroup$ – sunshine Feb 21 '20 at 5:22
  • $\begingroup$ Posted a general code for this below. If you need to go further, let me know. $\endgroup$ – Todd Burus Feb 21 '20 at 6:20

In general, if your data (called dat) has two factors A and B and a covariate C, then you can test using an ANCOVA as follows.

Anova(aov(y~C + A*B, data=dat, contrasts=list(A=contr.sum, B=contr.sum)), type="III")

This lets you test the significance of the interaction terms (i.e. a null that they are all equal to 0) by looking at the test statistic and P-value for A:B.


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