# R/kstest2d.R In sppmix: Modeling Spatial Poisson and Related Point Processes

#### Documented in kstest2dkstest2dsurf

#' Nonparametric Goodness-of-fit test between two point patterns
#'
#' @description
#' This function performs a two-dimensional Kolmogorov-Smirnov goodness-of-fit
#' test on two point patterns.
#'
#' For examples see
#'
#' \url{http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html
#' #kstest2d}
#'
#' @param x1,x2 Objects of class \code{\link[spatstat]{ppp}}.
#' @param showinfo Logical variable. Requests to display the test conclusion based on the value of the p-value. Default is TRUE.
#' @return A list with class "htest" containing the following components:
#' \item{statistic}{Value of the KS statistic}
#' \item{p.value}{The p-value of the test}
#' \item{alternative}{A character string describing the alternative hypothesis}
#'
#' @references Peacock, J.A. (1983). Two-dimensional goodness-of-fit testing in
#' astronomy. Monthly Notices Royal Astronomy Society, 202, 615-627.
#'
#' Adapted from Matlab code by Dylan Muir.
#' @author Jiaxun Chen, Sakis Micheas
#' @examples
#' \donttest{
#' # generate two point patterns
#' mix1 <- rnormmix(3, sig0 = .01, df = 5, xlim=c(0, 5), ylim=c(0, 5))
#' intsurf1=to_int_surf(mix1,lambda = 40, win =spatstat::owin( c(0, 5),c(0, 5)))
#' mix2 <- rnormmix(8, sig0 = .01, df = 10, xlim=c(0, 5),ylim=c(0, 5))
#' intsurf2=to_int_surf(mix2,lambda = 50, win =spatstat::owin( c(0, 5),c(0, 5)))
#' #generate patterns from the two different models
#' pp1 <- rsppmix(intsurf1)
#' pp2 <- rsppmix(intsurf2)
#' pp3 <- rsppmix(intsurf2)#pp3 is from the same model as pp2
#' # Test for goodness of fit, p-value should be small
#' kstest2d(pp1, pp2)
#' # Test for goodness of fit, p-value should be large
#' kstest2d(pp2, pp3)}
#'
#' @export
kstest2d <- function(x1, x2,showinfo=TRUE) {
ecdf2d <- function(x,edge){
count=cbind(as.numeric(x$x>=edge[1] & x$y>=edge[2]),
as.numeric(x$x<=edge[1] & x$y>=edge[2]),
as.numeric(x$x<=edge[1] & x$y<=edge[2]),
as.numeric(x$x>=edge[1] & x$y<=edge[2]))
return(count)
}

if((!spatstat::is.ppp(x1)) | (!spatstat::is.ppp(x2))){
stop("Point pattern must be of class ppp.")
}
n1 <- length(x1$x) n2 <- length(x2$x)
pp1 <- cbind(x1$x,x1$y)
pp2 <- cbind(x2$x,x2$y)
KSstat=-Inf
for (i in 1:(n1+n2)){
if(i <= n1) {
edge=pp1[i,]
} else {
edge=pp2[i-n1,]
}
vfcdf1=apply(ecdf2d(x1,edge), 2, sum)/n1
vfcdf2=apply(ecdf2d(x2,edge), 2, sum)/n2
vfdiff=abs(vfcdf2-vfcdf1)
fksts=max(vfdiff)
if (fksts>KSstat) KSstat=fksts
}

#Peacock Z calculation and P estimation
n <- n1*n2/(n1+n2)
Zn <- sqrt(n) * KSstat
Zinf <- Zn/(1-0.53 * n^(-0.9))
pValue=2*exp(-2 * (Zinf - 0.5)^2)
# Clip invalid values for P
if (pValue > 1.0)
pValue = 1.0;

RVAL <- list(statistic=c(KSstat = KSstat),p.value = pValue,method =
"Two sample 2D Kolmogorov-Smirnov Test",data.name =
paste(sep=", ",deparse(substitute(x1)),deparse(substitute(x2))
),alternative="No goodness of fit")
class(RVAL) <- "htest"
if(showinfo)
{
if (RVAL$p.value > .1) { cat("\nLarge p-value, the patterns arise from the SAME model.\n") } else { cat("\nSmall p-value, the patterns arise from DIFFERENT models.\n") } } return(RVAL) } #' Nonparametric Goodness-of-fit test for a point pattern against a surface #' #' @description #' This function performs a two-dimensional Kolmogorov-Smirnov goodness-of-fit #' test for a point pattern against a given intensity surface. #' #' For examples see #' #' \url{http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html #' #kstest2dsurf} #' #' @param pp Object of class \code{\link[spatstat]{ppp}}. #' @param intsurf Object of class \code{intensity_surface}. #' @param truncate Requests to truncate the generated point patterns #' to be within the window of the intensity object intsurf. Default is FALSE. #' @param iters Number of point patterns to generate and #' compare against \code{pp}. The larger this #' value is, the more test performed, and thus the #' more reliable the result. #' #' @author Sakis Micheas #' @seealso \code{\link{rmixsurf}}, #' \code{\link{kstest2d}}, #' \code{\link{rsppmix}}, #' \code{\link{plotmix_2d}} #' @examples #' \donttest{ #' # generate two intensity surfaces; assume the same window [-3,3]x[-3,3] #' mixsurf1 <- rmixsurf(m = 3, lambda=100,xlim=c(-3,3),ylim=c(-3,3)) #' plot(mixsurf1) #' mixsurf2 <- rmixsurf(m = 5, lambda=200,xlim=c(-3,3),ylim=c(-3,3)) #' plot(mixsurf2) #' #generate point patterns from the two different models #' pp1 <- rsppmix(mixsurf1, truncate=FALSE) #' plotmix_2d(mixsurf1,pp1,colors=TRUE) #' pp2 <- rsppmix(mixsurf2, truncate=FALSE) #' plotmix_2d(mixsurf2,pp2,colors=TRUE) #' # Test for goodness of fit, p-value should be small #' kstest2d(pp1, pp2) #' # Test each pattern for gof against both Poisson models #' kstest2dsurf(pp1, mixsurf1)#correct model for pp1 #' kstest2dsurf(pp1, mixsurf2)#wrong model for pp1 #' kstest2dsurf(pp2, mixsurf2)#correct model for pp2 #' kstest2dsurf(pp2, mixsurf1)#wrong model for pp2} #' #' @export kstest2dsurf <- function(pp, intsurf,truncate = FALSE,iters=500) { if((!spatstat::is.ppp(pp)) || (!is.intensity_surface(intsurf))) { stop("Bad point pattern object or intensity surface object.") } gof=0 for(i in 1:iters) { pp1 <- suppressWarnings(rsppmix(intsurf,truncate = truncate)) res=kstest2d(pp, pp1,showinfo=FALSE) if(res$p.value>.1)
gof=gof+1
}
cat("K-S suggests gof",100*gof/iters,"% of the time.\n",iters,"K-S tests performed.")
}


## Try the sppmix package in your browser

Any scripts or data that you put into this service are public.

sppmix documentation built on Jan. 13, 2021, 10:04 p.m.