[R] envelope in spatstat

[Adrian Baddeley]

Tom Richardson wrote:

> This query is regarding the use of the 'envelope' function in Spatstat.

For queries about a contributed package, please contact the package author.

> My data can be represented as a point process with CONTINUOUS marks:
> 
> points <- ppp(x=x,y=y, marks=m, window= wind)
> 
> However the marks are alignments (lines), and so have to be treated
> differently to normal scalar marks. Hence to create a mcf object with the
> appropriate test function for alignment marks, I input 'func' (below)
> suggested by Stoyan & Penttinen (1989):
> 
> func <- function(m1,m2) { sin(abs(m1-m2))^2}
> mcf <- markcorr(points, func, normalise = TRUE, method="density")
> 
> So far, so good. 

I presume the "alignment" values 'm' are numerical values that represent 
angles (in radians). Otherwise the code above would be incorrect, 
because sin() works on angles in radians.

> Using 'envelope' and 'rlabel' I would like to
> check if the pattern in the data is  lost when randomly relabeling the mark
> for each point. If the test function, 'func' were the usual G(m1,m2)=m1*m2,
> then the following would work:
> 
> E <- envelope(points,  markcorr, nsim=20, simulate=expression(rlabel(points)))
> 
> However, in the above 'markcorr' calculates  G(m1,m2)=m1*m2  by default. So
> need to specify the test function to be G(m1,m2)= sin(abs(m1-m2))^2.
> 
> [ ....]
 >
> So how do I tell 'envelope' that I want to specify the mark correlation
> test function ??

First look at help(markcorr). It says that the test function for the 
mark correlation is specified by the second argument to markcorr(), 
which has the formal name 'f'. So you want 'f' to be 'func'.

Now look at help(envelope). It says that the arguments "..." will be 
passed through to the function that is evaluated for each point pattern. 
So call envelope(points, markcorr, f=func) and the argument f=func will 
be passed to markcorr.

Thus:
        E <- envelope(points,  markcorr, nsim=20, f=func,
              simulate=expression(rlabel(points)))


Two other comments:
    (1) for efficiency you can replace sin(abs(m1-m2))^2
         by  sin(m1-m2)^2.

    (2) if you take nsim=20 then the significance level is
        1/21 = 4.8% (one-sided) or 2/21 = 9.5% (two-sided).
        If you want a significance level of 5%, you need either
        nsim=19 (one-sided) or nsim=39 (two-sided).

Adrian Baddeley