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 : 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.
Adrian Baddeley
2004-May-31 02:10 UTC
[R] Re: Rhelp: Need help interpreting plots in spatstat
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)