Multivariate gradient analysis for matrix nestedness

The following R code implements a permutation-based multivariate approach to gradient analysis for inferring the relative importance of putative drivers of matrix nestedness.  (Previous univariate methods of gradient analysis assume that explanatory variables are uncorrelated.)  In the paper we used a structural equation model (path analysis) to disentangle the direct and indirect influences of the row- and column-specific variables (since many variables were causally dependent), but multiple regression can also be used if no indirect effects are suspected.

Source: Novak, Moore & Leidy (2011). Nestedness patterns and the dual nature of community reassembly in California streams: a multivariate permutation-based approach. Global Change Biology, 17, 3714-3723.

Download: PermuationGradientAnalysis.R

library(vegan) # required for NODF calculation

########### Define a few functions ###########
# Spearman correlation of observed order to sorted order
spear<-function(x){ cor(seq(1,length(x)),order(x,decreasing=T),method='spearman') }

# NODF function from vegan with friendly output
NODF<-function(mat){out<-t(as.matrix(nestednodf(mat,order=FALSE)$statistic[c(3,2,1)])); colnames(out)<-c('NODF','NODF.row','NODF.col');out}

# Permutation function that spits out NODF and associated spearman correlations
rshuffle.spear<-function(M,R,C,order=FALSE,printall=TRUE,reps=3){
  Out<-dim(0)
  r<-nrow(M)
  c<-ncol(M)
  for (p in 1:reps){
    if(order==FALSE){
      temp.row.order<-sample(r)
      temp.col.order<-sample(c)
    }
    if(order==TRUE){
      temp.row.order<-order(apply(M,1,sum),decreasing=TRUE)
      temp.col.order<-order(apply(M,2,sum),decreasing=TRUE)
    }
  temp.M<-M[temp.row.order,temp.col.order]
  temp.R<-R[temp.row.order,]
  temp.C<-C[,temp.col.order]
  NODF.out<-NODF(temp.M)
  Out<-rbind(Out,cbind(NODF.out, t(as.matrix(apply(temp.R,2,spear))), t(as.matrix(apply(temp.C,1,spear)))))
  Out<-data.frame(Out)
  colnames(Out)[-c(1:3)]<-c(colnames(R),rownames(C))
  temp.M<-cbind(temp.M,apply(temp.M,1,sum))
  temp.M<-rbind(temp.M,apply(temp.M,2,sum))
  if(printall==TRUE){print(temp.M); print(temp.R); print(temp.C)}
 }
 return(Out)
}
########### # Make up some data & analyze it ###########
m<-n<-8
fill<-0.46
M<-matrix(rbinom(n*m,1,fill),nrow=n,byrow=T)
M<-M[order(apply(M,1,sum),decreasing=TRUE),order(apply(M,2,sum),decreasing=TRUE)]
R<-matrix(c(sample(seq(0.1,n/10,0.1)),sample(seq(1,n,1)),sample(seq(10,n*10,10))),nrow=n,byrow=F)
C<-matrix(c(sample(seq(0.1,m/10,0.1)),sample(seq(1,m,1)),sample(seq(10,m*10,10))),ncol=n,byrow=T)
colnames(R)<-LETTERS[1:ncol(R)]
rownames(R)<-paste('Sp',seq(1:n),sep='')
rownames(C)<-letters[1:nrow(C)]
colnames(C)<-paste('Site',seq(1:m),sep='')
rshuffle.spear(M,R,C,order=TRUE)
rshuffle.spear(M,R,C,order=FALSE,printall=FALSE)