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

Date: Mon 31 May 2004 - 12:10:46 EST

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

```

Dear Lisa

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

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

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