title: summing 2 distance matrices for getting a third 'overall' distance matrix (ecological context)
body: <p>I am ecologist, using mainly the vegan R package.</p>
  
  <p>I have 2 matrices (sample x abundances) (See data below):</p>
  
  <blockquote>
    <blockquote>
      <p>matrix 1/ nrow= 6replicates*24sites, ncol=15 species abundances (fish)
      matrix 2/ nrow= 3replicates*24sites, ncol=10 species abundances (invertebrates)</p>
    </blockquote>
  </blockquote>
  
  <p>The sites are the same in both matrices. I want to get the overall bray-curtis dissimilarity (considering both matrices) among pairs of sites. I see 2 options:</p>
  
  <p>option 1, averaging over replicates (at the site scale) fishes and macro-invertebrates abundances, cbind the two mean abundances matrix (nrow=24sites, ncol=15+10 mean abundances) and calculating bray-curtis.</p>
  
  <p>option 2, for each assemblage, computing bray-curtis dissimilarity among pairs of sites, computing distances among sites centroids. Then summing up the 2 distance matrix.</p>
  
  <p>In case I am not clear, I did these 2 operations in the R codes below.</p>
  
  <p>Please, could you tell me if the option 2 is correct and more appropriate than option 1.</p>
  
  <p>thank you in advance.</p>
  
  <p>Pierre</p>
  
  <h1>here is below the R code exemples</h1>
  
  <h1>generating data</h1>
  
  <pre><code>library(plyr);library(vegan)
  
  #assemblage 1: 15 fish species, 6 replicates per site
  a1.env=data.frame(
    Habitat=paste("H",gl(2,12*6),sep=""),
    Site=paste("S",gl(24,6),sep=""),
    Replicate=rep(paste("R",1:6,sep=""),24))
  
  summary(a1.env)
  
  a1.bio=as.data.frame(replicate(15,rpois(144,sample(1:10,1))))
  
  names(a1.bio)=paste("F",1:15,sep="")
  
  a1.bio[1:72,]=2*a1.bio[1:72,]
  
  #assemblage 2: 10 taxa of macro-invertebrates, 3 replicates per site
  
  a2.env=a1.env[a1.env$Replicate%in%c("R1","R2","R3"),]
  
  summary(a2.env)
  
  a2.bio=as.data.frame(replicate(10,rpois(72,sample(10:100,1))))
  
  names(a2.bio)=paste("I",1:10,sep="")
  
  a2.bio[1:36,]=0.5*a2.bio[1:36,]
  
  
  #environmental data at the sit scale
  
  env=unique(a1.env[,c("Habitat","Site")])
  
  env=env[order(env$Site),]
  </code></pre>
  
  <h1>OPTION 1, averaging abundances and cbind</h1>
  
  <pre><code>a1.bio.mean=ddply(cbind(a1.bio,a1.env),.(Habitat,Site),numcolwise(mean))
  
  a1.bio.mean=a1.bio.mean[order(a1.bio.mean$Site),]
  
  a2.bio.mean=ddply(cbind(a2.bio,a2.env),.(Habitat,Site),numcolwise(mean))
  
  a2.bio.mean=a2.bio.mean[order(a2.bio.mean$Site),]
  
  bio.mean=cbind(a1.bio.mean[,-c(1:2)],a2.bio.mean[,-c(1:2)])
  
  dist.mean=vegdist(sqrt(bio.mean),"bray")
  </code></pre>
  
  <h1>OPTION 2, computing for each assemblage distance among centroids and summing the 2 distances matrix</h1>
  
  <pre><code>a1.dist=vegdist(sqrt(a1.bio),"bray")
  
  a1.coord.centroid=betadisper(a1.dist,a1.env$Site)$centroids
  
  a1.dist.centroid=vegdist(a1.coord.centroid,"eucl")
  
  a2.dist=vegdist(sqrt(a2.bio),"bray")
  
  a2.coord.centroid=betadisper(a2.dist,a2.env$Site)$centroids
  
  a2.dist.centroid=vegdist(a2.coord.centroid,"eucl")
  </code></pre>
  
  <blockquote>
    <blockquote>
      <blockquote>
        <p>summing up the two distance matrices using Gavin Simpson 's fuse()    </p>
      </blockquote>
    </blockquote>
  </blockquote>
  
  <pre><code>dist.centroid=fuse(a1.dist.centroid,a2.dist.centroid,weights=c(15/25,10/25))
  </code></pre>
  
  <blockquote>
    <blockquote>
      <blockquote>
        <p>summing up the two euclidean distance matrices (thanks to Jari Oksanen correction)</p>
      </blockquote>
    </blockquote>
  </blockquote>
  
  <pre><code>dist.centroid=sqrt(a1.dist.centroid^2 + a2.dist.centroid^2)
  </code></pre>
  
  <h1>and the 'coord.centroid' below for further distance-based analysis (is it correct ?)</h1>
  
  <pre><code>coord.centroid=cmdscale(dist.centroid,k=23,add=TRUE)
  </code></pre>
  
  <h1>COMPARING OPTION 1 AND 2</h1>
  
  <pre><code>pco.mean=cmdscale(vegdist(sqrt(bio.mean),"bray"))
  
  pco.centroid=cmdscale(dist.centroid)
  
  comparison=procrustes(pco.centroid,pco.mean)
  
  protest(pco.centroid,pco.mean)
  </code></pre>