[R-sig-Geo] Estimating multiStrauss interaction radii

[Rolf Turner]

On 17/08/12 14:11, "José M. Blanco Moreno" wrote:

     <SNIP>
>> (c) Gamma parameters greater than 1 are OK for interactions between 
>> points of
>> different type.  For points of the same type they are not just 
>> undesirable, they
>> are verboten.  The model is undefined if any gamma_{ii} > 1.
> So if I understand it properly, I can discard any model fitting any 
> gamma_{ii} > 1. However, that seems to be happening to me/my data 
> quite often. If it is not asking too much, should I take it as an 
> indication of what? Recalcitrant data, improperly chosen model, any of 
> the previous?
If you get a gamma_{ii} estimate bigger than 1, then this is basically 
saying that
gamma_{ii} is *equal* to 1, i.e. there is no interaction amongst points 
of type i.

If there is indeed no interaction then the estimate of gamma_{ii} will 
be bigger
than 1 as often as not --- or even oftener!  E.g.:

set.seed(42)
X <- rpoispp(200)
marks(X) <- factor(sample(c("a","b"),npoints(X),TRUE))
M <- matrix(0.05,2,2)
FIT <- ppm(X,inter=MultiStrauss(radii=M))

You get gamma_{11} = 1.2386.  Which doesn't really make sense.  (You get
gamma_{22} = 0.8116, which is OK.

In such cases you should "project" the fitted model onto the valid 
parameter space:

     FIT_proj <- ppm(X,inter=MultiStrauss(radii=M),project=TRUE)

The log pseudolikelihood for FIT_proj is a wee bit smaller than the 
(essentially meaningless)
log pseudolikelihood for FIT.  (This will always happen.)

I guess that in your profiling you should set project=TRUE in your calls 
to ppm().

But I'm actually now a bit puzzled by the results that I get from using 
project=TRUE
in the forgoing artificial example.  I need to do some more checking 
into this.  I'll
get back to you on this issue ASAP.
>
>
>>
>> (d) I have no idea how to handle the problem of getting two local maxima
>> which are widely separated (and of similar heights, presumably).
> Yes, I did not make the point clear enough. Since they are about the 
> same height, changing the grid a bit can indicate one or the other, 
> depending on "alignment". But I wasn't aware (or I did not remember) 
> that they have to accour at interpoint distances (likely that I 
> haven't read enough), so I will check again.

The crucial quantity here is "s_{ij}(X)" = the number of r-close pairs 
in the pattern X
whose members are of type i and type j respectively.  As you change r, 
the values of
s_{ij}(X) will tend to change --- as r increases there will be more 
r-close pairs.  But
the number of r-close pairs can increase only when r "crosses" a 
boundary equal to
an interpoint distance.  And if none of the s_{ij}(X) change, the model 
fit doesn't
change.

     <SNIP>
> Just after sending the e-mail I experimented with optim, but not too 
> much. Anyway, the "preliminary results" were simply discouraging, but 
> yes, I would like to know how you did, because my force brute approach 
> using optim did not work, plainly, it sucks.

     I didn't try using the optim() approach on MultiStrauss models; I 
*think* I tried
     it on Geyer models, optimising over the interaction radius and the 
saturation
     constant.  .... just went away and searched for what I did, and I 
see in my notes
     the comment "Doesn't work worth a tinker's dam!"

     Now thinking on it a bit more, I have *vague* recollections that I 
tried using
     method="L-BFGS-B" and had a *bit* more success, but it wasn't 
really very
     good.  I can't find any trace of my efforts in this regard so I 
think I must have
     given it up as a bad job and scrubbed any code that I wrote.

     So basically you're right.  This idea sucks, and I don't think it's 
worth pursuing
     it any further.

         cheers,

             Rolf