I'm trying to use multivariate techniques to compare two datasets (same structure) that were collected using different sampling techniques. I'd like to compute a PCoA for the first dataset (D1), and then see how the data from dataset 2 (D2) compares, relative to the distance matrix computed for D1. Essentially, I want to calculate PCoA scores for D2, based on the distance matrix of D1.

I'm using the "vegan" package in R.

Thanks in advance.

#calculate distance matrix and PCoA for "D2"
serb.d2 <- vegdist(serbbiomB2, "bray") #distance matrix for D2
serb.pcoa.2 <- cmdscale(serb.d2, k=4) #PCoA scores for D2
#how do I calculate PCoA scores for "D1", based on "serb.d2"?
  • $\begingroup$ What is your question? $\endgroup$ Commented Nov 7, 2014 at 20:40
  • $\begingroup$ I think I understand conceptually what you want to do, at a higher level, but you can't really do anything to a dataset with the dissimilarities from the first data set. I explain a directly analogue of what you want to do using a constrained PCoA in my answer and suggest an alternative which avoids even doing the embedding in principal coordinates. $\endgroup$ Commented Nov 7, 2014 at 21:42

1 Answer 1


If you are going to go down that route, I would suggest you look at capscale() in the vegan package that you are using.

Put all your x and y data in the same object, say XY. Then have another data frame (DF) with the indicator variable you need, say technique, which has levels awesome and cheap (say) --- essentially, this is just a single column with values awesome or cheap indicating which technique was used for each row of XY.

Now you can do

fit <- capscale(XY ~ technique, data = DF, distance = "bray")

This will fit a constrained PCoA (constrained analysis of principal coordinates) and test whether the amount of variation in the dissimilarities of your data that can be explained by which technique was used to generate the data is larger than you'd expect to see if there were no difference between techniques.

Here's an example from vegan using the Dutch dune data set

data(dune, dune.env)
fit <- capscale(dune ~ Management, data = dune.env, distance = "bray")

This gives us

> fit
Call: capscale(formula = dune ~ Management, data = dune.env, distance =

              Inertia Proportion Rank
Total          4.2990                
Real Total     4.5940     1.0000     
Constrained    1.5000     0.3265    3
Unconstrained  3.0940     0.6735   14
Imaginary     -0.2950               5
Inertia is squared Bray distance 

Eigenvalues for constrained axes:
  CAP1   CAP2   CAP3 
0.8998 0.4533 0.1471 

Eigenvalues for unconstrained axes:
  MDS1   MDS2   MDS3   MDS4   MDS5   MDS6   MDS7   MDS8   MDS9  MDS10  MDS11 
1.2730 0.4874 0.3784 0.3100 0.2100 0.1507 0.0855 0.0729 0.0601 0.0322 0.0172 
 MDS12  MDS13  MDS14 
0.0098 0.0044 0.0024


> set.seed(23)
> anova(fit)
Permutation test for capscale under reduced model
Permutation: free
Number of permutations: 999

Model: capscale(formula = dune ~ Management, data = dune.env, distance = "bray")
         Df Variance     F Pr(>F)   
Model     3   1.5001 2.586  0.003 **
Residual 16   3.0938                
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The anova() is using a permutation test to compare the ratio of Model to Residual variance (as indicated by the pseudo $F$ entry in the table).

If you do this with the current CRAN version of vegan you'll see a slightly different display for anova() --- I've got the development version (2.2-0 to be) running here.

You could also do this analysis on the raw dissimilarity matrices without embedding them in an Euclidean space (what PCoA does), using the adonis() function in vegan. See its help page ?adonis for details.


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.