[R-sig-eco] Fwd: change dissimilarity index in custom function

cuca smith cuca2305 at gmail.com
Fri Jan 9 13:12:19 CET 2015


Dear all,

I want to fit a GLM for a typical 'distance decay' pattern. I'm going to
use the following custom function to obtain the model parameters through
bootstrapping (Millar et al. 2011. Ecology), but I need to change the
"binary" distance measure to other dissimilarity indices (Simpson and
Sorensen, using e.g. betapart package). In my case, counts is a
presence-absence matrix and gradient is an environmental variable.

Bstrap=function(gradient,counts,nboots=1000,measure="binary",dist="dist",
                like.pairs=T) {
  dist=match.fun(dist)
  if(length(gradient)!=nrow(as.matrix(counts)))
    stop("\nDifferent number of sites detected in input data")
  nsites=length(gradient)
  Results=matrix(NA,nrow=nboots,ncol=2)
  nPairs=rep(NA,nboots)
  for(i in 1:nboots) {
    indices=sample(1:nsites,replace=T)
    d=as.vector(dist(gradient[indices],method="euclidean"))
    s=as.vector(1-dist(counts[indices,],method=measure))
    df=data.frame(d,s)
    n=nrow(df)
    #remove like site pairs (zero distance) if requested
    if(!like.pairs) df=subset(df,subset=c(d>0))
    n.used=nrow(df)
    if(n.used<2) stop("Less than 2 site-pairs. Can not fit glm")
    #log-binomial fit
    g=glm(s~d,family=binomial(link=log),data=df)
    Results[i,]=coef(g); nPairs[i]=n.used }
  #Get parameters of interest
  Results=cbind(Results[,1],-Results[,2],
              exp(Results[,1]),-log(2)/Results[,2])
  colnames(Results)=c("a","beta","s0","halfd")
  if(!like.pairs) {
   cat("\n",100*(n-mean(nPairs))/n,"% removed due to like pairs")
   cat("\n Bstrap s.e.s have been corrected for the smaller no. of
pairs\n\n") }
  Variances=apply(Results,2,wt.var,w=nPairs/n)*mean(nPairs)/n
  Summary=rbind(apply(Results,2,weighted.mean,w=nPairs/n),sqrt(Variances))
  rownames(Summary)=c("Bstrap mean","Bstrap s.e.")
  return(list(Fits=Results,Summary=Summary,
               CtrlList=list(nboots=nboots,like.pairs=like.pairs)))
}

Any suggestion will be helpful!!

Thanks in advance,
Silvia

	[[alternative HTML version deleted]]



More information about the R-sig-ecology mailing list