R/Linda.R

In Linda: Line Transect-Based Nearest Neighbor Distance Analysis

##############################################################################
#distributional aggregation measure for line transect-derived biodiversity data
#Line transect-based nearest neighboring distance analysis
##############################################################################


##############################################################################
### The Cluster-Random-Regular Continuum
#  Clark, P. J., & Evans, F. C. (1954). Distance to nearest neighbor as a measure of spatial relationships in populations. Ecology, 35, 445-453.
#  R is the Clark and Evans competition index
#  the input xy contains x and y coordinates of recorded organisms in a single line transect
#
CECI <- function(xy, area=NULL, method = "NP")
{
mat=xy

  if(is.data.frame(mat) | is.matrix(mat))
  {
  n = dim(mat)[1] # numbers of points
  }else
  {
  n=length(mat)
  }
  if(is.null(area))
  {
  xmin=min(xy[,1])
  xmax=max(xy[,1])
  ymin=min(xy[,2])
  ymax=max(xy[,2])
  area=(xmax-xmin)*(ymax-ymin)
  }else
  {
  area = area   # study area (square meter)
  }
  Rho = n/area  # density of points per square meter
  # minimal distance between a target point to its neighbours
  if(method == "NP"){
    dd <- as.matrix(dist(mat, p=2))
    diag(dd) <- NA
    dmin = apply(dd, 1, min, na.rm = TRUE)
  } else if (method == "P"){
    
    t1 <- Sys.time()
    #library(doParallel)
    #detectCores()
    x=0
    cl <- makeCluster(as.integer(detectCores()/3))
    registerDoParallel(cl)
    Res <- foreach(x = 1:n,.combine = c) %dopar% {
      d <- sqrt((mat[x,1] - mat[,1])^2 + (mat[x,2] - mat[,2])^2)
      dmin = sort(d)[2]
      dmin
    }
    t2 <- Sys.time()-t1
    t2
    stopCluster(cl)
    dmin = Res
  }
  # R is the Clark and Evans competition index
  R = (1/n) * sum(dmin) * 2 * sqrt(Rho)
  names(R) = "Clark & Evans Competition Index"
  sigma=.26136/sqrt(n*Rho)
  t1 = (1/n) * sum(dmin)
  t2 = 1/(2*sqrt(Rho))
  C=(t1-t2)/sigma
  
  return(list(R = R, t1=t1,t2=t2,c=C,p = 2*pnorm(-abs(C))))
}


###################################################################################################
#find nearest neighbor distance in a single line transect and conduct hypothesis testing
#xy is two column matrix, the first is x, and the second is y
#collected from sequential sampling, so xy are ordered in a time-forward way (backward is fine)
#in other words, rows of xy matrix should be sorted according to sampling sequential orders.
#L is the length of the line transect, can be given or calculated from the original data (if L=NULL)
#L=NULL is recommended, because the given line transect length might over-estimate the aggregation pattern!
NND <- function(xy,L=NULL){
  n = dim(xy)[1]
  li=vector()
  dd=0
  if(n > 2)
  {
	li[1]=sqrt((xy[1,1]-xy[2,1])^2+(xy[1,2]-xy[2,2])^2)
	dd=dd+li[1]
	for(i in 2:(n-1))
	{
	d1=sqrt((xy[i,1]-xy[i+1,1])^2+(xy[i,2]-xy[i+1,2])^2)
	d2=sqrt((xy[i,1]-xy[i-1,1])^2+(xy[i,2]-xy[i-1,2])^2)
	li=c(li,min(d1,d2))
	dd=dd+d1
	}#i
	li=c(li,sqrt((xy[n-1,1]-xy[n,1])^2+(xy[n-1,2]-xy[n,2])^2))
	#dd=dd+li[n]
	###############
	if(is.null(L))
	{
	L=dd #excluding the first and last distance (which are used repeatedly)
	}
  }else if(n==2)
  {
    li=sqrt((xy[1,1]-xy[2,1])^2+(xy[1,2]-xy[2,2])^2)
	if(is.null(L))
	{
	L=li
	}
  }
  ###############
  if(n>1)
  {
  ra=mean(li)
  lam <- (n-1)/L
  re <- 1/(2*lam)
  sig <- 1/(2*lam*sqrt(n)) 
  c <- (ra-re)/sig
  return(list(R = ra/re, ra = ra, re = re, sig=sig,c=c,p = 2*pnorm(-abs(c)),df=n,L=L))
  }else
  {
  return(0)
  }
}##
##



#######################################################################################
#multiple line transects for conduct hypothesis testing
#lxy is three-column matrix, the first one is the line-transect ID
#the second and third ones are x and y respectively
LNND <- function(lxy){
lid=unique(lxy[,1])
N=dim(lxy)[1]
K=length(lid)
t1_lst=t2_lst=sigre=vector()
for(i in 1:K)
{
ids=which(lxy[,1]==lid[i])
res=NND(lxy[ids,-1])
t1_lst[i]=res$ra*(res$df)
t2_lst[i]=res$re*(res$df)
sigre[i]=res$sig^2*(res$df)^2 #variance of each line transect
}#i
sigj=sqrt(sum(sigre))/N #standard error
#sigj=sqrt(sum(sigre))/sqrt(N) #pool-level standard deviation
  #sig <- 1/(2*lam*sqrt(n)) 
  #sig2 <- n*(1/(2*lam))^2 #sigma_lE

c=((sum(t1_lst)/(N))-(sum(t2_lst)/(N)))/sigj

#####
return(list(LR = sum(t1_lst)/sum(t2_lst), t1 = sum(t1_lst)/N, t2 = sum(t2_lst)/N, 
sig=sigj,c =c , p = 2*pnorm(-abs(c)),df = N))
}##
##


################################################################################
## the estimator of pi by Solow (2000)
## z is a vector of sequential species' label in the line-transect sample
#see Chen et al. (2019) MEE
v.est = function(z){
  m = length(z)
  TT = sum(z[-m]==z[-1])
  v = TT/(m-1)
  est = v
  est
}


################################################################################
#### pai estimator
## z is a vector of sequential species' label in the line-transect sample
##see Chen et al. (2023) Frontiers in Plant Science
pi.est = function(z){
  m = length(z)
  Xi=table(z) ## Xi is the species abundance data
  TT = sum(z[-m]==z[-1])
  v = TT/(m-1)
  b = m*(1-v)
  Y=m^2-sum(Xi^2)
  if((m+1)^2*b^2>8*b*Y) {
    est = 1-((m+1)*b+sqrt((m+1)^2*b^2-8*b*Y))/2/Y
  } else {
    est = v
  }
  if(v==1 || est > v){
    est = v
  }
  if(est<0) est = 0
  est
}