[R-sig-eco] Species Area Curve (recreating EstimateS output)

[stephen sefick]

I would like some advice on the "appropriate" method for constructing
species area curves.  I have written a function based on the vegan
function specaccum, which I will include later in the email.  I would
like to fit the michalis-menton equation to the data to extrapolate to
more area than we have sampled.  When should I fit the nls model?  I
am more or less trying to duplicate functionality of EstimateS, trying
to carry out my advisors wishes, and above all trying to accurately
assess the area needed to sample effectively the benthic diversity in
coastal plain streams of the Southeastern United States.

take 100 permutations of the species accumulation (randomly sample the
stations and build the accumulation from this)
then:

1.take the mean of the permuted accumulation data and fit the model
2.fit the model to all of the permutations separately and then take
the mean of the predicted
3.fit the model to all of the permuted data together

Thanks very much for any and all help

Stephen

function for 2 above:
#feed regular vegan community matrix into this...

perm_species <- function(comm, sites=6, permutations=1000,
predict=100, surber_num=3, surber_size=0.092){
    surber_station <- surber_num*surber_size
    x <- comm
    x <- as.matrix(x)
    x <- x[, colSums(x) > 0, drop=FALSE]
    n <- nrow(x)
    p <- ncol(x)
    if (p == 1) {
        x <- t(x)
        n <- nrow(x)
        p <- ncol(x)
    }
    accumulator <- function(x, ind) {
        rowSums(apply(x[ind, ], 2, cumsum) > 0)
    }
        perm <- array(dim = c(n, permutations))
        for (i in 1:permutations) {
            perm[, i] <- accumulator(x, sample(n))
        }
	
	
	nls.fit <- array(dim = c(predict, permutations))
	predict2 <- seq(from=surber_station, by=surber_station, length.out=predict)	
	for(i in 1:permutations){	
	sites2 <- 1:sites	
	square_meters <- seq(from=surber_station, by=surber_station ,length.out=sites)	
	s <- data.frame(x=square_meters, y=perm[,i])	
	#o <- nls(y~a*(x^b), start=list(a=10, b=1), data=s)	
	o <- nls(y~((a*x)/(b+x)), start=list(a=10, b=1), data=s)	
	nls.fit[,i] <- (coef(o)["a"]*predict2)/(coef(o)["b"]+predict2)	
	}
	
	sites <- 1:n
        specaccum <- apply(nls.fit, 1, mean)
        sdaccum <- apply(nls.fit, 1, sd)		
	error <- qnorm(0.975)*sdaccum/sqrt(permutations)	
	left <- specaccum-error
	right <- specaccum+error	
	out <- data.frame(specaccum, sdaccum, left, right)
	out
}

-- 
Stephen Sefick

Let's not spend our time and resources thinking about things that are
so little or so large that all they really do for us is puff us up and
make us feel like gods.  We are mammals, and have not exhausted the
annoying little problems of being mammals.

								-K. Mullis