[R-sig-Geo] directional mantel correlogram

[Tyler Smith]

Hi,

First off, I'm just learning spatial statistics, so I apologise if I'm
not using the correct terminology.

I'm working on an ecological problem, and I'm trying to distinguish
between two models of community assembly. The first is a completely
intrinsic process, with no relationship between species and the
environment. As such, I expect this process should produce an
isotropic pattern. The alternative process links species abundances to
the underlying environmental gradient, which is spatially structured.
In this case I think the community composition should be anisotropic,
since it is based on an underlying spatial gradient.

It seems to me that I could evaluate this problem with a directional
mantel correlogram, using a measure of community similarity contrasted
against varying distance and direction classes. I found the mgram()
function in package ecodist that computes omnidirectional mantel
correlograms. I've extended this function by adding a second loop to
'mask' the distance matrix for each of the selected angles. I've
pasted my work below, along with a simple test case (which runs in a
minute or two on my machine).

Questions for r-geo gurus:

Is this an appropriate approach to the problem?
Is my implementation correct?

Thanks for your time and patience,

Tyler

dmgram <-
  function (species.d, space.xy, breaks, nclass, stepsize, 
  	   nperm = 1000, mrank = FALSE, nboot = 500, pboot = 0.9, 
	   cboot = 0.95, alternative = "two.sided", trace = FALSE, 
	   alpha = c(0, 45, 90, 135))
{
  space.d <- dist(space.xy)
  space.ang <- pw.angle(space.xy)
  species.d <- as.vector(species.d)
  space <- as.vector(space.d)

  ## calculate break points for specified angles, given that breaks
  ## are at the midway between successive angles
  
  ang.breaks <- numeric(length(alpha))

  for (i in 1:(length(alpha) - 1))
    ang.breaks[i] <- mean(alpha[i:(i+1 )])

  ang.breaks[length(ang.breaks)] <- mean(c(180 + alpha[1], 
  				    	   alpha[length(alpha)]))

  ## any angles between the last break point and 180 are flipped so
  ## that the first angle class is contiguous
  space.ang[space.ang > ang.breaks[length(ang.breaks)]] <- 
            space.ang[space.ang > ang.breaks[length(ang.breaks)]] - 180
  
  if (missing(breaks)) {
    if (missing(nclass)) {
      if (missing(stepsize)) {
	## reduce the calculated number of classes to account for
    	## additional angle categories
        nclass <- round(1 + 3.3 * log10(length(space.d))/length(alpha))
        stepsize <- (max(space.d) - min(space.d))/nclass
      }
      else {
        nclass <- round((max(space.d) - min(space.d))/stepsize)
      }
    }
    else {
      if (missing(stepsize)) {
        stepsize <- (max(space.d) - min(space.d))/nclass
      }
    }
    breaks <- seq(0, stepsize * nclass, stepsize)
  }
  else {
    nclass <- length(breaks) - 1
  }

  answer.m <- matrix(0, ncol = 7, nrow = nclass * length(alpha))
  dimnames(answer.m) <- list(NULL, c("lag", "ngroup", "mantelr", 
  		     		   "pval", "llim", "ulim","dir"))
  answer.m[, 4] <- rep(1, nrow(answer.m))

  for (j in 1:length(ang.breaks)) {
    if (j == 1) {
      amin <- -180  ## set the minimum value for the first angle class
    } else
      amin <- ang.breaks[j - 1]
    amax <- ang.breaks[j]
    space.aclass <- matrix(1, nrow = nrow(space.ang), 
    		    	      ncol = ncol(space.ang))

    space.aclass[space.ang <= amin] <- 0
    space.aclass[space.ang > amax] <- 0
    diag(space.aclass) <- 0
    space.aclass <- as.dist(space.aclass)
    ## space.aclass serves as a mask to zero out all distances that
    ## are not in the current angle class    

    for (i in 1:nclass) {
      row.num <- (j - 1) * nclass + i
      dmin <- breaks[i]
      dmax <- breaks[i + 1]
      answer.m[row.num, 1] <- (dmin + dmax)/2
      space.dclass <- rep(0, length(space.d))
      space.dclass[(space.d * space.aclass) <= dmin] <- 1
      space.dclass[(space.d * space.aclass) > dmax] <- 1
      ngroup <- length(space.dclass) - sum(space.dclass)
      answer.m[row.num, 2] <- ngroup
      if (ngroup > 0) {
        mant <- mantel(species.d ~ space.dclass, nperm = nperm, 
	     	mrank = mrank, nboot = nboot, pboot = pboot, 
		cboot = cboot)
        answer.m[row.num, 3] <- mant[1]
        if (alternative == "two.sided") 
          answer.m[row.num, 4] <- mant[4]
        else answer.m[row.num, 4] <- mant[2]
        answer.m[row.num, 5] <- mant[5]
        answer.m[row.num, 6] <- mant[6]
      }
      answer.m[row.num,7] <- alpha[j]
      if (trace) 
        cat(i, "\t", answer.m[i, 2], "\t", answer.m[i, 3], 
            "\n")
    }
  }
  
  results <- list(mgram = answer.m, resids = NA)
  class(results) <- "dmgram"
  results
}

pw.angle <- function(x) {
  res <- matrix(NA, nrow = nrow(x), ncol = nrow(x))
  for (i in 1:(nrow(x)-1))
    for (j in (i+1):nrow(x))
      res[i,j] <- atan((x[i,2] - x[j,2])/(x[i,1] - x[j,1])) * 180 / pi
  res[lower.tri(res)] <- t(res)[lower.tri(t(res))]
  diag(res) <- 0
  res[res < 0] = 180 + res[res < 0]
  return(res)
}

plot.dmgram <- function (x, pval = 0.05, xlab = "Distance", 
	       		ylab = "Mantel r", ...) {
  x <- x$mgram
  pval.v <- x[, 4]
  ang.v <- x[,7]
  ang.lev <- unique(ang.v)
  ang.num <- length(ang.lev)
  
  plot(x[,1], x[,3], type = 'n', xlab = xlab, ylab = ylab, ...)
  
  for (i in 1:ang.num) {
    points(x[x[,7] == ang.lev[i], 1], x[x[,7] == ang.lev[i], 3], 
    	     col = i, type = 'l')
    points(x[pval.v <= pval & ang.v == ang.lev[i], 1], 
             x[pval.v <= pval & ang.v == ang.lev[i], 3], 
	     pch = 16, col = i)
    points(x[pval.v > pval & ang.v == ang.lev[i], 1], 
    	     x[pval.v > pval & ang.v == ang.lev[i], 3], 
	     pch = 1, col = i)
  }
  invisible()
}

## test case for an artificial isotropic pattern:

coords <- cbind(rep(1:20, each = 20), rep(1:20, 20)) 
coords2 <- cbind(jitter(coords[,1], factor = 10), jitter(coords[,2], 
	   	   factor = 10))
coords.dist <- dist(coords)
coords2.dist <- dist(coords2)
dmangram.coords <- dmgram(coords2.dist, coords)

plot(dmangram.coords)