BIOENV (Clarke and Ainsworth 1993)

bio.env.r

bio.env <- function(fix.mat, var.mat, 
fix.dist.method="bray", var.dist.method="euclidean",
scale.fix=FALSE, scale.var=TRUE,
output.best=10,
var.max=ncol(var.mat)
){
	if(dim(fix.mat)[1] != dim(var.mat)[1]){stop("fixed and variable matrices must have the same number of rows")}
	if(var.max > dim(var.mat)[2]){stop("var.max cannot be larger than the number of variables (columns) in var.mat")}

	require(vegan)

	combn.sum <- sum(factorial(ncol(var.mat))/(factorial(1:var.max)*factorial(ncol(var.mat)-1:var.max)))
	
	if(scale.fix){fix.mat<-scale(fix.mat)}else{fix.mat<-fix.mat}
	if(scale.var){var.mat<-scale(var.mat)}else{var.mat<-var.mat}
	fix.dist <- vegdist(fix.mat, method=fix.dist.method)
	RES_TOT <- c()
	best.i.comb <- c()
	iter <- 0
	for(i in 1:var.max){
		var.comb <- combn(1:ncol(var.mat), i, simplify=FALSE)
		RES <- data.frame(var.incl=rep(NA, length(var.comb)), n.var=i, rho=0)
		for(f in 1:length(var.comb)){
			iter <- iter+1
			var.dist <- vegdist(as.matrix(var.mat[,var.comb[[f]]]), method=var.dist.method)
			temp <- suppressWarnings(cor.test(fix.dist, var.dist, method="spearman"))
			RES$var.incl[f] <- paste(var.comb[[f]], collapse=",")
			RES$rho[f] <- temp$estimate
			if(iter %% 100 == 0){print(paste(round(iter/combn.sum*100, 3), "% finished"))}
		}

		order.rho <- order(RES$rho, decreasing=TRUE)
		best.i.comb <- c(best.i.comb, RES$var.incl[order.rho[1]])
		if(length(order.rho) > output.best){
			RES_TOT <- rbind(RES_TOT, RES[order.rho[1:output.best],])
		} else {
			RES_TOT <- rbind(RES_TOT, RES)
		}
	}
	rownames(RES_TOT)<-NULL
		
	if(dim(RES_TOT)[1] > output.best){
		order.by.best <- order(RES_TOT$rho, decreasing=TRUE)[1:output.best]
	} else {
		order.by.best <- order(RES_TOT$rho, decreasing=TRUE)
	}
	OBB <- RES_TOT[order.by.best,]
	rownames(OBB) <- NULL

	order.by.i.comb <- match(best.i.comb, RES_TOT$var.incl)
	OBC <- RES_TOT[order.by.i.comb,]
	rownames(OBC) <- NULL

	out <- list(
		order.by.best=OBB,
		order.by.i.comb=OBC,
		best.model.vars=paste(colnames(var.mat)[as.numeric(unlist(strsplit(OBB$var.incl[1], ",")))], collapse=",") ,
		best.model.rho=OBB$rho[1]
	)
	out
}