CODEWITHSUNDEEP
rspatstatmodel-fitting
Camille Gontier
Camille Gontier

Reputation: 155

Constraints on the parameters of a fitted model in spatstat

In the R package spatstat, I would like to impose some constraints on the parameters of a model fitted on my data.

My data is a multitype ppp object which types are called BRP and GluR. As can be seen on the plot below, points are (roughly) organized into concentric circles. BRP objects are localized on a ring with radius 128.4, while GluR are organized on an outer ring with radius 261.8.enter image description here

Using the ppm function, I would thus like to fit a model where the intensity function for type BRP scales with ~abs(sqrt(x^2+y^2)-128.4), and where the intensity function for type GluR scales with ~abs(sqrt(x^2+y^2)-261.8). However, when using the following trend formula:

fit_NMJ <- ppm(NMJ,~marks*(abs(sqrt(x^2+y^2)-128.4)+abs(sqrt(x^2+y^2)-261.8)))

I obtain cross-interactions between GluR and ~abs(sqrt(x^2+y^2)-radius_BRP ), and between BRP and ~abs(sqrt(x^2+y^2)-radius_GluR) (which I want to avoid):

> radius_GluR <- 261.8
> radius_BRP <- 128.4
> fit_NMJ <- ppm(NMJ,~marks*(abs(sqrt(x^2+y^2)-radius_BRP )+abs(sqrt(x^2+y^2)-radius_GluR )))
> fit_NMJ
Nonstationary multitype Poisson process

Possible marks: ‘BRP’ and ‘GluR’

Log intensity:  ~marks * (abs(sqrt(x^2 + y^2) - radius_BRP) + abs(sqrt(x^2 + y^2) - radius_GluR))

Fitted trend coefficients:
                                 (Intercept)                                    marksGluR 
                                -9.850558359                                  0.231315028 
           abs(sqrt(x^2 + y^2) - radius_BRP)           abs(sqrt(x^2 + y^2) - radius_GluR) 
                                -0.044039428                                  0.009055454 
 marksGluR:abs(sqrt(x^2 + y^2) - radius_BRP) marksGluR:abs(sqrt(x^2 + y^2) - radius_GluR) 
                                 0.039652754                                 -0.018198305 

                                                 Estimate       S.E.       CI95.lo      CI95.hi Ztest       Zval
(Intercept)                                  -9.850558359 1.92806275 -13.629491901 -6.071624817   *** -5.1090445
marksGluR                                     0.231315028 2.00684749  -3.702033781  4.164663836        0.1152629
abs(sqrt(x^2 + y^2) - radius_BRP)            -0.044039428 0.01749010  -0.078319402 -0.009759454     * -2.5179626
abs(sqrt(x^2 + y^2) - radius_GluR)            0.009055454 0.01392686  -0.018240694  0.036351602        0.6502150
marksGluR:abs(sqrt(x^2 + y^2) - radius_BRP)   0.039652754 0.01784817   0.004670991  0.074634517     *  2.2216710
marksGluR:abs(sqrt(x^2 + y^2) - radius_GluR) -0.018198305 0.01478904  -0.047184295  0.010787685       -1.2305263

i.e. I would like to constraint the parameters marksGluR:abs(sqrt(x^2 + y^2) - radius_BRP) and marksBRP:abs(sqrt(x^2 + y^2) - radius_GluR) to be zero. Is it possible to specify this in the trend formula ? Does it also work for multiple independent datasets (when using the mppm function) ?

I also have the option to separately fit these two types, but when fitting different non-Poisson models with point interactions I might miss some inter-type interactions (as, for instance, when fitting a MultiStrauss model).

Upvotes: 0

Views: 69

Answers (1)

Adrian Baddeley
Adrian Baddeley

Reputation: 2973

You just need to construct canonical covariates directly.

First let's abbreviate

g <- function(x, y, marks, R, type) {
       I(marks == type) * abs(sqrt(x^2+y^2) - R)
     }
gGluR <- function(x, y, marks) {
    g(x, y, marks, radius_GluR, "GluR")
}
gBRP <- function(x, y, marks) {
    g(x, y, marks, radius_BRP, "BRP")
}

Then

ppm(NMJ ~ marks + gGluR + gBRP)

Upvotes: 2

Related Questions