NNid <- function(xy, group=NULL) {
 
# find id of nearest neighbor, given xy positions as a matrix
# uses either of two functions,
# depending on whether splancs or spatialstats is available

# if group is non null, then finds the neighbor that is withing
#   the same group as the point.  Typical use is restricting
#   neighbors to same side of the street.
#  group only works with spatial stats

# if splancs, uncomment the following line
  nndistG(xy)$neighs

# if spatial stats, uncomment the following lines
#  if (is.null(group)) {
#    temp <- find.neighbor(xy,k=2);
#    as.vector(temp[temp[,1] != temp[,2],2]);		
# remove self neighbors, return id of neighbor
#  N.B. assumes find.neighbor returns list in order of points 1, 2,
#    }
#   else {
#     temp <- find.neighbor(xy,k=10);
#	# find more neighbors than needed
#     temp <- temp[temp[,1] != temp[,2],];
	# remove self 
#     temp <- temp[group[temp[,1]] == group[temp[,2]],];
	# keep only pairs in same group
#     nndist <- tapply(temp[,3], temp[,1], min);
	# minimum distance for each point
#     nnpos <- tapply(temp[,3],temp[,1]);
	# gives cell number for each set of points
     
#     as.vector(temp[temp[,3]==nndist[nnpos],2]);
	# return the point # for the neighbor that has the smallest 
	# nndist among nieghbors in same group
#     }
  } 



NNinfo <- function(xy,label, group = NULL) {

  #  Info for nearest neighbor test of segregation, Naa only
  #  data is xy coordinates in xy, labels in label
  #    group (if want to restrict neighbors to only consider
  #           points in same group as the point )

  # output is a list with Naa = obs. count, 
  #   ENaa = expected count
  #   VarNaa = variance of count
  #  Nbb = obs count,
  #   ENbb = expected count
  #   VarNbb = variance of count
  #   CovNaaNbb = covariance

  nnid <- NNid(xy, group);		# get the nearest neighbor id

  n <- table(label)
  N <- sum(n);

  l <- sort(unique(label));

  if (length(l) > 2) {stop("Currently, only 2 labels supported")};
  if (dim(xy)[1] != N) {stop("Number of labels not same as number
of points")};
  
  Naa <- sum((label==l[1])*(label[nnid] == l[1]));
  ENaa <- n[1]*(n[1]-1)/(N-1);

  Nbb <- sum((label==l[2])*(label[nnid] == l[2]));
  ENbb <- n[2]*(n[2]-1)/(N-1);

  R <- sum((1:N) == nnid[nnid]);	# number of refl NN
  Q1 <- table(c(1:6,table(nnid)))
  Q <- 2*Q1[2] + 6*Q1[3] + 12*Q1[4] + 20*Q1[5] + 30*Q1[6] - 70;
    # get # shared neighbors,  extra bit (1:6) is to ensure
     # that outer table() always generates a length 6 vector
    #  correct for the extra by subtracting (2+6+12+20+30)=70
  P11 <- ENaa / N;
  P111 <- P11 * (n[1] - 2)/(N-2);
  P1111 <- P111 * (n[1] - 3) / (N - 3);

  VarNaa <- (N+R)*P11 + (2*N-2*R + Q)*P111 + 
		(N*(N-3)-Q+R)*P1111 - ENaa^2;

  P22 <- ENbb / N;
  P222 <- P22 * (n[2] - 2)/(N-2);
  P2222 <- P222 * (n[2] - 3) / (N - 3);
  P1122 <- P11*n[2]*(n[2]-1) / ((N-2)*(N-3));

  VarNbb <- (N+R)*P22 + (2*N-2*R + Q)*P222 + 
		(N*(N-3)-Q+R)*P2222 - ENbb^2;

  CovNaaNbb <- (N*N - 3*N - Q + R)*P1122 - ENaa*ENbb;

  list(Naa = Naa,ENaa = ENaa, VarNaa = VarNaa, 
       Nbb=Nbb, ENbb = ENbb, VarNbb = VarNbb, CovNaaNbb = CovNaaNbb,
       labels = l )
  }

NNtest <- function(l) {
# compute test statistics from list output from NNinfo

  M <- matrix(NA,nrow=3,ncol=3,
	dimnames=list(c(paste("Z ",l$labels[1],l$labels[1],sep=''),
			paste("Z ",l$labels[2],l$labels[2],sep=''),"2 d.f."),
	   c("Test. stat.","2 sided P value", "1 sided P value")))
  M[1,1] <- (l$Naa - l$ENaa) / sqrt(l$VarNaa);
  M[2,1] <- (l$Nbb - l$ENbb) / sqrt(l$VarNbb);
  r <- l$CovNaaNbb / sqrt(l$VarNaa * l$VarNbb);
  M[3,1] <- (M[1,1]^2 + M[2,1]^2 - 2*r*M[1,1]*M[2,1])/(1-r*r);
  M[1,2] <- 2*(1-pnorm(abs(M[1,1])));
  M[1,3] <- 1-pnorm(M[1,1]);
  M[2,2] <- 2*(1-pnorm(abs(M[2,1])));
  M[2,3] <- 1-pnorm(M[2,1]);
  M[3,2] <- 1-pchisq(M[3,1],2);
  M
  }