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

About this list Date view Thread view Subject view Author view Attachment view

From: Adrian Baddeley (adrian@maths.uwa.edu.au)
Date: Mon 31 May 2004 - 12:10:46 EST


Message-id: <16570.37926.351650.674832@maths.uwa.edu.au>

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)

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


About this list Date view Thread view Subject view Author view Attachment view

This archive was generated by hypermail 2.1.3 : Mon 31 May 2004 - 23:05:14 EST