[R] Re: Rhelp: Need help interpreting plots in spatstat

[Adrian Baddeley]

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)