R/Aindex.R

In GmAMisc: 'Gianmarco Alberti' Miscellaneous

#'R function for calculating the Hodder-Okell's A index of spatial association
#'
#'The function allows to calculate the Hodder-Okell's A index of spatial association between the
#'features of two point patterns.\cr
#'
#'The functions takes as input two point patterns (SpatialPointDataframe class) and calculate the A
#'index. Details about the latter are provided by:\cr Orton C. 1980, "Mathematics in Archeology",
#'Glasgow: William Collins Sons & Co Ltd, pp. 154-155\cr Blankholm P. 1990, "Intrasite spatial
#'Analysis in Theory and Practice", Aarhus: Aarhus University Press, pp. 130-135.\cr
#'
#'The A index is about equal to 1 when the two patterns are randomly mingled; it is smaller than 1
#'when the two patterns are segregated; it is larger than 1 when the features of the two point
#'patterns tend to occur together. The computational details are provided by Blankholm's book cited
#'above (page 132).\cr
#'
#'The significance of the A index is calculated via the randomized approach devised by:\cr Kintigh K
#'W. 1990, “Intrasite Spatial Analysis: A Commentary of Major Methids”. In Voorrips A, “Mathematics
#'and Information Science in Archaeology: A Flexible Framework”, Studies in Modern Archaeology 3:
#'165-200\cr
#'
#'Given two patterns A and B being analysed, the procedure keeps the points location unchanged and
#'randomly assigns the points to either pattern. The random re-assignment is performed B times (199
#'by default) and each time the A index is calculated. One-tailed and two-tailed p values are
#'calculated following the procedure described by Baddeley et al., "Spatial Point Patterns.
#'Methodology and Applications with R", CRC Press 2016, p. 387.\cr
#'
#'@param x Point pattern (SpatialPointDataframe class).
#'@param y Point pattern (SpatialPointDataframe class).
#'@param B Number of permutations (199 by default).
#'@param studyplot Feature (of polygon type; SpatialPolygonsDataFrame class) representing the study
#'  area; if not provided, the study area is internally worked out as the bounding polygon based on
#'  the union the convex hulls of the x and y patterns. This is only used for visualization purpose,
#'  should the user want to plot the two point patterns within the actual study area.
#'@param addmap FALSE (default) or TRUE if the user does not want or wants a map of the study area
#'  and of the two patterns to be displayed.
#'
#'@return The function produces:\cr -an histogram showing the frequency distribution of the
#'  randomized A index, with vertical reference lines representing the 0.025th and 0.975th quantile
#'  of the distribution. A black dot represents the observed A index. At the bottom of the chart the
#'  randomized p values are reported;\cr -optionally (setting the 'addmap' parameter to TRUE), a map
#'  showing the point patterns (and the study area, if supplied).
#'
#'@keywords Aindex
#'
#'@export
#'
#'@importFrom graphics hist rug abline points par layout text dotchart axis
#'@importFrom rgeos gConvexHull gDistance
#'@importFrom sp coordinates
#'@importFrom raster union pointDistance plot text
#'@importFrom utils txtProgressBar setTxtProgressBar
#'@importFrom stats quantile
#'@importFrom stats AIC as.formula binomial coef coefficients confint cooks.distance dbinom dist fitted formula glm hatvalues kruskal.test median na.omit pbinom predict.glm qnorm quantile r2dtable rstandard sd t.test update
#'@importFrom methods as
#'
#' @examples
#' # calculate the Hodder-Okell's A index for the two patterns, and plot the map
#' Aindex(springs, points, addmap=TRUE)
#'
#'@seealso \code{\link{distRandSign}} , \code{\link{distCovarModel}} ,
#'  \code{\link{pointsCovarModel}}
#'
Aindex <- function(x, y, studyplot=NULL, B=199, addmap=FALSE){

  if(is.null(studyplot)==TRUE){
    # build the convex hull of the union of the convex hulls of the two features
    ch <- gConvexHull(union(gConvexHull(x), gConvexHull(y)))
    region <- ch
  } else {
    region <- studyplot
  }

  #size of patter x
  nx <- length(x)

  #interpoint distance matrix
  distxx <- gDistance(x, x, byid=TRUE)

  #sum the interpoint distance matrix
  sum.distxx <- sum(distxx)

  #calculation of the "average" interpoint distance for pattern x
  rxx <- sum.distxx / (nx * (nx-1))

  #size of patter y
  ny <- length(y)

  #interpoint distance matrix
  distyy <- gDistance(y, y, byid=TRUE)

  #sum the interpoint distance matrix
  sum.distyy <- sum(distyy)

  #calculation of the "average" interpoint distance for pattern y
  ryy <- sum.distyy / (ny * (ny-1))

  #matrix of the distance between points in pattern x to the points in pattern y
  distxy <- gDistance(x, y, byid=TRUE)

  #sum of the interpoint distance matrix
  sum.distxy <- sum(distxy)

  #calculation of the "average" interpoint
  rxy <- sum.distxy / (nx * ny)

  #Aindex
  Aind <- (rxx * ryy) / (rxy^2)

  #pool the coordinates of points in pattern x and y to be used later within the pointDistance function
  pooledData <- as.matrix(rbind(coordinates(x), coordinates(y)))

  #set the progress bar to used inside the loop
  pb <- txtProgressBar(min = 0, max = B, style = 3)

  #create an empty sloth for the permuted A index
  Aind.perm <- numeric(length=B)

  for(i in 1:B){
    #create a random set of number to be used as index to extract the values from the pooledData matrix
    #the set of number has size equal to the number of points in pattern x
    index <- sample(1: (nx + ny), size=nx, replace=F)
    #extract from the pooledData matrix the rows corresponding to the generated random index
    x.perm <- pooledData[index,]
    #extract from the pooledData matrix all the rows that do not correspond to the random index
    #the procedure eventually creates two sets of points randomly shuffling the entries of the pooledData matrix
    y.perm <- pooledData[-index,]

    #repeat all the steps needed to calculate the A index using the random samples
    distxx.perm <- pointDistance(x.perm, lonlat=FALSE, allpairs=TRUE)
    sum.distxx.perm <- sum(distxx.perm)
    rxx.perm <- sum.distxx.perm / (nx * (nx-1))

    distyy.perm <- pointDistance(y.perm, lonlat=FALSE, allpairs=TRUE)
    sum.distyy.perm <- sum(distyy.perm)
    ryy.perm <- sum.distyy.perm / (ny * (ny-1))

    distxy.perm <- pointDistance(x.perm, y.perm, lonlat=FALSE, allpairs=TRUE)
    sum.distxy.perm <- sum(distxy.perm)
    rxy.perm <- sum.distxy / (nx * ny)
    Aind.perm[i] <- (rxx.perm * ryy.perm) / (rxy.perm^2)

    setTxtProgressBar(pb, i)
  }

  #calculate the single tailed and the two-tailed p values
  p.lower <- (1 + sum (Aind.perm < Aind)) / (1 + B)
  p.upper <- (1 + sum (Aind.perm > Aind)) / (1 + B)
  two.tailed.p <- 2 * (min(p.lower, p.upper))

  if(addmap==TRUE){
    par(mfrow=c(1,2))
    plot(region,
         main="Map of point patterns plus study area",
         cex.main=0.9, col=NA,
         border="red",
         lty=2,
         axes=TRUE)
    plot(x,
         add=TRUE,
         pch=20)
    plot(y,
         add=TRUE,
         pch=20,
         col="red")
  }

  #plot the histogram of the permuted distirbution of the A index
  graphics::hist(Aind.perm, main=paste0("Freq. distribution of Hodder-Okell's A index\n across ", B, " permutations"),
       sub=paste0("A-index: ", round(Aind, 3),"\n p-value segregated: ", p.lower, "; p-value associated: ", p.upper, "\np-value different from randomly mingled: ", two.tailed.p, "\nA<1 (segregated); A=1 (randomly minlged); A>1 (associated)"),
       xlab="",
       cex.main=0.90,
       cex.sub=0.70,
       cex.axis=0.85)
  rug(Aind.perm, col = "#0000FF")
  abline(v=stats::quantile(Aind.perm, 0.025), lty=2, col="blue")
  abline(v=stats::quantile(Aind.perm, 0.975), lty=2, col="blue")
  points(x=Aind, y=0, pch=20, col = "black")

  # restore the original graphical device's settings
  par(mfrow = c(1,1))
}