R/spseg.R

In spatialkernel: Non-Parametric Estimation of Spatial Segregation in a Multivariate Point Process

#' Integrated Functions for Spatial Segregation Analysis
#' Spatial segregation analysis to be performed by a single function and
#' presentations by associated \code{plot} functions.
#' @param pts an object that contains the points. This could be a two-column
#'     matrix or a ppp object from spatstat.
#' @param marks numeric/character vector of the types of the point in the data.
#' @param h numeric vector of the kernel smoothing bandwidth at which to
#'     calculate the cross-validated log-likelihood function.
#' @param opt integer, 1 to select bandwidth; 2 to calculate type-specific 
#'     probabilities; and 3 to do the Monte Carlo segregation test.
#' @param ntest integer with default 100, number of simulations for the 
#'     Monte Carlo test.
#' @param poly matrix containing the \code{x,y}-coordinates of the polygonal
#'     boundary of the data.
#' @param delta spacing distance of grid points at which to calculate the estimated
#'     type-specific probabilities for \code{\link{image}} plot.
#' @param proc logical with default \code{TRUE} to print the
#'     processing message.
#' @param obj list of the returning value of \code{spseg}.
#' @param types numeric/character types of the marks of data points to plot
#'     the estimated type-specific probabilities, default to plot all types.
#' @param sup logical with default \code{FALSE}, if \code{TRUE} to superimpose data points 
#'     on the estimated type-specific probability surface.
#' @param quan numeric, the pointwise significance levels to add contours to 
#'     \code{\link{image}} plot of the estimated type-specific probability
#'     surface, with default of \code{c(0.05, 0.95)}.
#' @param col list of colors such as that generated by \code{risk.colors}.
#' @param breaks a set of breakpoints for the \code{col}: must give one more 
#'     breakpoint than colour.
#' @param \dots other arguments concerning \code{\link{plot}} and 
#'     \code{\link{points}}
#' @details \code{spseg} implements a complete spatial segregation analysis by
#'   selecting bandwidth, calculating the type-specific probabilities, and then
#'   carrying out the Monte Carlo test of spatial segregation and pointwise
#'   significance. Some \code{plot} functions are also provided here so that
#'   users can easily present the results.
#'   
#'   These functions are provided only for the convenience of users. Users can
#'   instead use individual functions to implement the analysis step by step and
#'   plot the diagrams as they wish.
#'   
#'   Examples of how to use \code{spseg} and present results using \code{plot}
#'   functions are presented in \code{\link{spatialkernel-package}}.
#' @return  A list with components 
#' \describe{
#'   \item{hcv}{bandwidth selected by the cross-validated log-likelihood
#'     function.}
#'   \item{gridx,gridy}{\code{x, y} coordinate vectors at which the grid points 
#'     are generated at which to calculate the type-specific probabilities 
#'     and pointwise segregation test \emph{p}-value.}
#'   \item{p}{estimated type-specific probabilities at grid points generated
#'     by vectors \code{gridx, gridy}.}
#'   \item{pvalue}{\emph{p}-value of the Monte Carlo spatial segregation test.}
#'   \item{stpvalue}{pointwise \emph{p}-value of the Monte Carlo spatial 
#'     segregation test.}
#'   \item{...}{copy of \code{pts, marks, h, opt}.}
#' }
#' @note
#'   Setting \code{h} to a unique value may force \code{spseg} to skip the
#'   selecting bandwidth step, go straight to calculate the type-specific 
#'   probabilities and then test the spatial segregation with this fixed
#'   value of bandwidth.
#' @seealso \code{\link{cvloglk}}, \code{\link{phat}}, \code{\link{mcseg.test}}, 
#'   \code{\link{pinpoly}}, \code{\link{risk.colors}}, and \code{\link{metre}}
#' @keywords spatial nonparametric hplot
#' @rdname spseg
#' @export
spseg.matrix <- function(pts, marks, h, opt=2, ntest=100, poly=NULL, 
  delta=min(apply(apply(pts, 2, range), 2, diff))/100, proc=TRUE, ...)
{
##use data within a polygon 
##select bandwidth, calculate phat, and do MC test 
##opt=1,2,3 to do task one, two or three
  if(!(opt %in% 1:3)) stop("\nargument opt must be one of 1, 2, or 3.\n")
    ans <- list(pts=pts, marks=marks, h=h, opt=opt)
    ## opt==1
    if(length(h) > 1) {
      if(proc) cat("\nCalculating cross-validated likelihood function\n")
      cv <- cvlogl(pts, marks, h)$cv
      hopt <- h[which.max(cv)]
      ans <- c(list(cv=cv, hcv=hopt), ans)
    } else {
      hopt <- h ## for opt=2, phat
    }
    if(opt >= 2) { ##phat and mcseg.test
      if(!is.null(poly)) ans$poly <- poly
      ## opt==2
      if(proc) cat("\nCalculating type-specific probabilities\n")
      mtypes <- unique(marks)
      m <- length(mtypes)
      if(is.null(poly)) {
        xyrng <- apply(pts, 2, range)
      } else xyrng <- apply(poly, 2, range)
      gridxy <- list(gridx=seq(xyrng[1,1]+delta/2, xyrng[2,1]-delta/2, by=delta),
        gridy=seq(xyrng[1,2]+delta/2, xyrng[2,2]-delta/2, by=delta))
      gridpts <- as.matrix(expand.grid(gridxy))
      ngrid <- nrow(gridpts)
      if(is.null(poly)) gridndx <- rep(TRUE, 1:ngrid) else 
        gridndx <- which(pinpoly(poly, gridpts)>0)
      p <- matrix(NA, ncol=m, nrow=ngrid)
      tmp <- phat(gridpts[gridndx,], pts, marks, hopt)$p
      p[gridndx,] <- tmp; colnames(p) <- colnames(tmp)
      ans <- c(gridxy, list(p=p), ans)
      if(opt==3) {
        ## opt==3
        if(proc) cat("\nMonte Carlo testing\n")
        stp <- matrix(NA, ncol=m, nrow=ngrid)
        mc <- mcseg.test(pts, marks, h, stpts=gridpts[gridndx,], 
          ntest=ntest, proc=proc)[c("pvalue", "stpvalue")]
        stp[gridndx,] <- mc$stpvalue; colnames(stp) <- colnames(mc$stpvalue)
        ans <- c(list(pvalue=mc$pvalue, stpvalue=stp, ntest=ntest), ans)
      }
    }
    if(proc) cat("\n")
    ans
}

#' @rdname spseg
#' @export
plotcv <- function(obj, ...) plot(obj$h, obj$cv, type="l", ...)

#' @rdname spseg
#' @export
plotphat <- function(obj, types=unique(obj$marks), sup=TRUE, col=risk.colors(10), 
    breaks=seq(0,1,length=length(col)+1), ...) 
{
    if(obj$opt<=1) stop("\nRun phatmctest() with argument opt=2 or 3.\n")
    sapply(types, function(x) match.arg(x, unique(obj$marks)))
    m <- length(types)
    for(j in 1:m) {
        if(is.null(obj$poly)) {
            plot(obj$pts[1,], xlab="", ylab="", xlim=range(obj$pts[,1]), ylim=range(obj$pts[,2]),
                 asp=1, main=types[j], type="n")
        } else {
            plot(obj$pts[1,], xlab="", ylab="", xlim=range(obj$poly[,1]), 
                 ylim=range(obj$poly[,2]), asp=1, main=types[j], type="n")
        }
        image(obj$gridx, obj$gridy, matrix(obj$p[,types[j]], ncol=length(obj$gridy)), 
              zlim=0:1, add=TRUE, col=col, breaks=breaks) 
        if(!is.null(obj$poly)) {
            polygon(obj$poly)
            if(sup) {
                ndx <- which(obj$marks == types[j]) ##not gridpts, pts!!!
                points(obj$pts[ndx,], ...)
            }
        }
        if(m>1 && j<m && interactive()) readline("\nPress Enter to plot next one ...")
    }
}

#' @rdname spseg
#' @export
plotmc <- function(obj, types=unique(obj$marks), quan=c(0.05, 0.95), sup=FALSE, 
    col=risk.colors(10), breaks=seq(0,1,length=length(col)+1), ...) 
{
    if(obj$opt<=2) stop("\nRun phatmctest() with argument opt=3 first.\n")
    sapply(types, function(x) match.arg(x, unique(obj$marks)))
    m <- length(types)
    for(j in 1:m) {
        if(is.null(obj$poly)) {
            plot(obj$pts[1,], xlab="", ylab="", xlim=range(obj$pts[,1]), ylim=range(obj$pts[,2]),
                 asp=1, main=types[j], type="n")
        } else {
            plot(obj$pts[1,], xlab="", ylab="", xlim=range(obj$poly[,1]), 
                 ylim=range(obj$poly[,2]), asp=1, main=types[j], type="n")
        }
        image(obj$gridx, obj$gridy, matrix(obj$p[,types[j]], ncol=length(obj$gridy)), 
              zlim=0:1, add=TRUE, col=col, breaks=breaks)
        if(!is.null(quan)) {
            contour(obj$gridx, obj$gridy,
                    matrix(obj$stpvalue[,types[j]], ncol=length(obj$gridy)),
                    levels=quan, add=TRUE, ...)
        } 
        if(!is.null(obj$poly)) {
            polygon(obj$poly)
            if(sup) {
                ndx <- which(obj$marks == types[j]) ##not gridpts, pts!!!
                points(obj$pts[ndx,], ...)
            }
        }
        if(m>1 && j<m && interactive()) readline("\nPress Enter to plot next one ...")
    }
}

##' Run the spatial segregation analysis
##'
##' This is the details of the S3 generic method
##' @title spseg
##' @return spseg results
##' @author Barry Rowlingson
##' @rdname spseg
##' @export
spseg <- function(pts, ...){
    UseMethod("spseg")
}

##' spseg for spatstat objects
##'
##' Does spseg for marked ppp objects
##' @title spatial segregation for ppp objects
##' @return an spseg object
##' @author Barry Rowlingson
##' @rdname spseg
##' @export
spseg.ppp = function(pts, h, opt, ...){
    pw = as.data.frame(spatstat::as.polygonal(pts$win))
    if(ncol(pw)>2){
        stop("spseg can only handle simple ring windows")
    }
    
    xypoly = as.matrix(pw)
    xypts = cbind(pts$x, pts$y)
    m = spatstat::marks(pts)
    # test if only one mark
    m = as.character(m)
    spseg.matrix(xypts, m, h=h, opt=opt, poly=xypoly, ...)
}