[S] Multitype version of Ripley's K Function

Adrian Baddeley (adrian@maths.uwa.edu.au)
Fri, 16 Oct 1998 09:52:33 +0800


Steve Friedman <friedman@gis.umn.edu> asked:

> Has anyone modified the K-function in the Spatial Stats module to
> account for more than one species?

Yes. Here is a simple version which I hand out in class.

It works only for rectangular windows, and uses the `border correction'
rather than Ripley's edge correction.

I'll post a more general version to statlib soon.

----
Adrian Baddeley, Mathematics, University of Western Australia
<http://maths.uwa.edu.au/~adrian/>

=====================================================================
#
# Kcross.S Estimation of cross K function
# for bivariate point pattern
#
# $Revision: 1.2 $ $Date: 1998/07/21 06:15:43 $
#
#
# -------- functions ----------------------------------------
#
# Kcross.b1() compute estimate of K by border correction,
# denominator = count of points
#
# Kcross.b2() compute estimate of K by border correction,
# denominator = area of eroded window
#
# crossdist() compute matrix of distances between
# each pair of data points
#
# bdist() compute vector of distances from
# each data point to boundary of window
#
#
# -------- standard arguments ------------------------------
# x1,y1 coordinates of input data points (type 1)
# x2,y2 coordinates of input data points (type 2)
# xrange,yrange dimensions of observation window
# r distance values at which to compute K
#
# -------- example -----------------------------------------
#
# x1 <- runif(20)
# y1 <- runif(20)
# x2 <- runif(30)
# y2 <- runif(30)
# plot(c(0,1),c(0,1),type="n")
# points(x1,y1,pch=2)
# points(x2,y2,pch=3)
# r <- seq(0,0.25,length=50)
# k <- Kest.b2(x1,y1,x2,y2,r=r)
# plot(r, k, type="l")
# lines(r, pi * r^2, lty=2)
#
# ------------------------------------------------------------------------

"bdist"<-
function(x, y, xrange=c(0,1), yrange=c(0,1))
{
xd <- pmin(x - min(xrange), max(xrange) - x)
yd <- pmin(y - min(yrange), max(yrange) - y)
d <- pmin(xd, yd)
return(d)
}

"crossdist"<-
function(x1, y1, x2, y2)
{
if(length(x1) != length(y1))
stop("lengths of x1 and y1 do not match")
if(length(x2) != length(y2))
stop("lengths of x2 and y2 do not match")
n1 <- length(x1)
n2 <- length(x2)
X1 <- matrix(rep(x1, n2), ncol = n2)
Y1 <- matrix(rep(y1, n2), ncol = n2)
X2 <- matrix(rep(x2, n1), ncol = n1)
Y2 <- matrix(rep(y2, n1), ncol = n1)
d <- sqrt((X1 - t(X2))^2 + (Y1 - t(Y2))^2)
return(d)
}

"Kcross.b1"<-
function(x1, y1, x2, y2, xrange=c(0,1), yrange=c(0,1), r)
{
if(length(x1) != length(y1))
stop("lengths of x1 and y1 do not match")
if(length(x2) != length(y2))
stop("lengths of x2 and y2 do not match")
n1 <- length(x1)
n2 <- length(x2)
width <- abs(diff(xrange))
height <- abs(diff(yrange))
area <- width * height
lambda1 <- n1/area
lambda2 <- n2/area
d <- crossdist(x1, y1, x2, y2)
b <- bdist(x1, y1, xrange, yrange)
k <- rep(0, length(r))
for(i in 1:length(r)) {
close <- (d <= r[i])
nclose <- apply(close, 1, sum)
bok <- (b > r[i])
numerator <- sum(nclose[bok])
denominator <- sum(bok)
k[i] <- numerator/denominator
}
k <- k/lambda2
return(k)
}

"Kcross.b2"<-
function(x1, y1, x2, y2, xrange=c(0,1), yrange=c(0,1), r)
{
if(length(x1) != length(y1))
stop("lengths of x1 and y1 do not match")
if(length(x2) != length(y2))
stop("lengths of x2 and y2 do not match")
n1 <- length(x1)
n2 <- length(x2)
width <- abs(diff(xrange))
height <- abs(diff(yrange))
area <- width * height
lambda1 <- n1/area
lambda2 <- n2/area
d <- crossdist(x1, y1, x2, y2)
b <- bdist(x1, y1, xrange, yrange)
numerator <- numeric(length(r))
denominator <- (width - 2 * r) * (height - 2 * r)
for(i in 1:length(r)) {
close <- (d <= r[i])
nclose <- apply(close, 1, sum)
bok <- (b > r[i])
numerator[i] <- sum(nclose[bok])
}
k <- numerator/(denominator * n1 * n2)
return(k)
}

-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu. To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message: unsubscribe s-news