Estimating multiStrauss interaction radii
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