R help - May 2004 - Rhelp: Need help interpreting plots in spatstat

Hello everybody--

I have been playing with my data in spatstat, and what I'd like to 
present is a basic exploratory spatial analysis. I have used the 
following code, using a ppp.object called tsdspoints.  The code 
develops the simulations and the envelopes I want, but  I don't 
understand my first plot here, the [tsds.ghat.short$r, 
tsds.ghat.short$raw]...I cobbled together this code from some class 
examples used by a prof of mine, and I don't really know how to 
interpret what is going on. The subsequent lines, I'm pretty sure, 
represent the empirical Ghat and the simulation envelope...but that 
first plot stumps me. 

You can see a pdf of the resulting plot if you download it from my 
website : www.bol.ucla.edu/~lschwei/sample.pdf WARNING: this is 
immediate download. 

THE CODE: ######

ghat.env <- function(n, s, r, win=owin(myOwin){
hold <- matrix(0, s, length(r))
for(i in 1:s){
hold[i,] <- Gest(runifpoint(n, win=myOwin), r=r)$raw
}
mn <- apply(hold, 2, mean)
Up <- apply(hold, 2, max)
Down <- apply(hold, 2, min)
return(data.frame(mn, Up, Down))
}
tsds.ghat <- Gest(tsdspoints)
tsds.ghat.short <- tsds.ghat[tsds.ghat$r <= .15, ]
tsds.cdf <- 1-exp(-40*pi*(tsds.ghat.short$r^2))
plot(tsds.ghat.short$r, tsds.ghat.short$raw, xlab="Distance
(degrees)",
ylab="Ghat", main=" Transport, Storage, and Disposal sites
(TSDs)")
lines(tsds.ghat.short$r, tsds.cdf)
tsds.env <- ghat.env(n=40, s=100, r=tsds.ghat$r)
tsds.env.short <- tsds.env[tsds.ghat$r <= .15, ]
lines(tsds.ghat.short$r, tsds.env.short$Up, lty=5)
lines(tsds.ghat.short$r, tsds.env.short$Down, lty=5)
###############

Many Thanks, 

Lisa
who was working on her stinking dissertation while everybody else was 
learning cool things in Vienna.

Dear Lisa

'spatstat' has extensive help files which would help you
to resolve this.

help(Gest) says that the output value of Gest() includes a 
vector Gest()$raw which is the raw (i.e. without edge correction)
estimate of the G function. It is not a very good estimate of G
because it can be severely biased. But it can be used for
hypothesis testing.

 > ghat.env <- function(n, s, r, win=owin(myOwin){
 > hold <- matrix(0, s, length(r))
 > for(i in 1:s){
 > hold[i,] <- Gest(runifpoint(n, win=myOwin), r=r)$raw
 > }
 > mn <- apply(hold, 2, mean)
 > Up <- apply(hold, 2, max)
 > Down <- apply(hold, 2, min)
 > return(data.frame(mn, Up, Down))
 > }

    This defines a function 'ghat.env' which generates 's'
independent
    simulated realisations of a pattern of 'n' points in the window
'myOwin',
    computes the raw estimate of the G function for each realisation,
    then takes the pointwise sample mean, maximum, and minimum of these
    estimates. If the estimates are denoted G1(r), ... Gn(r)
    then mn[i] is the sample mean of G1(r[i]), ...., Gn(r[i]).
    
 > tsds.ghat <- Gest(tsdspoints)
 > tsds.ghat.short <- tsds.ghat[tsds.ghat$r <= .15, ]

    These compute the G function estimates for your dataset 'tddspoints'
    and discard the estimates for r > 0.15.
    The result is a data frame containing several different estimates of G,
    including $raw.

 > tsds.cdf <- 1-exp(-40*pi*(tsds.ghat.short$r^2))

    This computes the theoretical G function for a Poisson process
    with intensity 40. 

 > plot(tsds.ghat.short$r, tsds.ghat.short$raw, xlab="Distance
(degrees)",
 > ylab="Ghat", main=" Transport, Storage, and Disposal sites
(TSDs)")

    This plots the RAW estimate of G() from your dataset 'tddspoints'
    against r. 

 > lines(tsds.ghat.short$r, tsds.cdf)

    This superimposes a plot of the theoretical G function for a Poisson
    process with intensity 40.

    ****NOTE**** THESE FUNCTIONS ARE NOT COMPARABLE!!!!
    The raw estimate is a severely biased estimate of G.
    
 > tsds.env <- ghat.env(n=40, s=100, r=tsds.ghat$r)
   
    This invokes the function ghat.env as described above,
    with n=40. Thus it simulates realisations of a random pattern
    of a FIXED number n=40 of points. 

    The results of these simulations are not strictly
    comparable with the Poisson process of intensity 40.
    A Poisson process has a random number of points.

    The mean number of points in the window equals 40 times the area
    of the window - I hope your window has area 1 or these simulations
    are completely incompatible with the 'tsds.cdf' curve.

 > tsds.env.short <- tsds.env[tsds.ghat$r <= .15, ]
   
    Discards estimates of G(r) for r > 0.15

 > lines(tsds.ghat.short$r, tsds.env.short$Up, lty=5)
 > lines(tsds.ghat.short$r, tsds.env.short$Down, lty=5)

    Plots the pointwise maxima and minima of the simulations
    in a classical Monte Carlo test plot.

--------------------------------------------------
I suggest you use the following code instead.

I've changed $raw to $trans which gives you the 'translation
correction',
yielding an unbiased estimator of G(r). If you really want the raw estimate,
edit $trans to $raw everywhere. 

regards
Adrian Baddeley


-------------------------------

ghat.env <- function(n, s, r, win=owin(myOwin){
   hold <- matrix(0, s, length(r))
   for(i in 1:s){
      hold[i,] <- Gest(runifpoint(n, win=win), r=r)$trans
   }
   mn <- apply(hold, 2, mean)
   Up <- apply(hold, 2, max)
   Down <- apply(hold, 2, min)
   return(data.frame(mn, Up, Down))
}

tsds.ghat <- Gest(tsdspoints)
tsds.ghat.short <- tsds.ghat[tsds.ghat$r <= .15, ]

# empirical G function 
plot(tsds.ghat.short$r, tsds.ghat.short$trans, 
     xlab="Distance (degrees)", ylab="Ghat", 
     main=" Transport, Storage, and Disposal sites (TSDs)",
     type="l", lwd=2)

# 100 simulations of n points 

forty <- tsdspoints$n
tsds.env <- ghat.env(n=forty, s=100, r=tsds.ghat$r)
tsds.env.short <- tsds.env[tsds.ghat$r <= .15, ]

# upper and lower envelopes
lines(tsds.ghat.short$r, tsds.env.short$Up, lty=5)
lines(tsds.ghat.short$r, tsds.env.short$Down, lty=5)

# mean G function from simulations of n points
lines(tsds.ghat.short$r, tsds.env.short$mn,lty=2)