Search code examples
rpointconfidence-intervalspatstatcox

Problem in estimating the confidence intervals of a log-gaussian cox process with kppm on spatstat on R


I am fitting a log gaussian cox point process with the kppm function on R. spatstat for point patterns. It is supposed to give me the intensity coefficients prediction and the confidence intervals both for the intensity of the latent non homogeneous poisson, then the coefficients of interaction of the cox process and the new confidence intervals (for the same coefficients) which take into account the interaction.

the commands i typed are:

quad_fes <- quadscheme(data = ppp_fes, dummy = test_dummy_points)
kppm_all_fes<-kppm(
 X=quad_fes,
  trend=~log(daylight_P50sp_im) + log(daylight_Hits_im ) + log(d_fes_im) + SpeedLimit_im + strada_tip_im + log(P1_im) + log(E1_im),
  clusters="LGCP", method="mincon")
kppm_all_fes
confint(kppm_all_fes)

The covariates are pixeled images (not defined on every point of the region, but it wasn't a problem for the estimate??) they are obtained from a rasterization with as.im. ppp_fes is a point pattern. type ppp

The R output doesn't give me the second set of confidence intervals (with interaction) what could it be? it just gives columns of NA

The output is

    Inhomogeneous Cox point process model
Fitted to point pattern dataset ‘quad_fes’
Fitted by minimum contrast
    Summary statistic: inhomogeneous K-function

Log intensity:  ~log(daylight_P50sp_im) + log(daylight_Hits_im) + log(d_fes_im) + SpeedLimit_im + 
strada_tip_im + log(P1_im) + log(E1_im)

Fitted trend coefficients:
           (Intercept) log(daylight_P50sp_im)  log(daylight_Hits_im)          log(d_fes_im) 
         -8.4336803232          -0.7195136965           0.0597531881          -0.1869696161 
         SpeedLimit_im          strada_tip_im             log(P1_im)             log(E1_im) 
         -0.0100578297           0.0005152514           0.3538175166          -0.0706547930 

Cox model: log-Gaussian Cox process
    Covariance model: exponential
Fitted covariance parameters:
     var    scale 
9.505158 6.238008 
Fitted mean of log of random intensity: [pixel image]

And for confint, it is

                          2.5 % 97.5 %
(Intercept)               NA     NA
log(daylight_P50sp_im)    NA     NA
log(daylight_Hits_im)     NA     NA
log(d_fes_im)             NA     NA
SpeedLimit_im             NA     NA
strada_tip_im             NA     NA
log(P1_im)                NA     NA
log(E1_im)                NA     NA

I can't find any info on Baddeley (spatial point patterns, 2015) on how to solve this.. I know it's a long shot but pleaseeeeeeee help I really want to find a solution.

I have tried using just one covariate at the time and checking how my covariates could be different from the examples on Baddeley's book but I cannot find a solution.

EDIT the warnings I get while using kppm are

Warning: Values of the covariates ‘daylight_P50sp_im’, ‘daylight_Hits_im’, ‘SpeedLimit_im’, ‘strada_tip_im’, ‘P1_im’, ‘E1_im’ were NA or undefined at 0.57% (7 out of 1236) of the quadrature points. Occurred while executing: ppm.quad(Q = X, trend = trend, covariates = covariates, forcefit = TRUE, Warning: Values for 2 query points lying outside the pixel image domain were estimated by projection to the nearest pixel

var(kppm_all_fes)

returns Error in var(kppm_all_fes) : is.atomic(x) is not TRUE

cvar(kppm_all_fes)

returns

(Intercept) log(daylight_P50sp_im) log(daylight_Hits_im) log(d_fes_im) SpeedLimit_im
(Intercept)                     NA                     NA                    NA            NA            NA
log(daylight_P50sp_im)          NA                     NA                    NA            NA            NA
log(daylight_Hits_im)           NA                     NA                    NA            NA            NA
log(d_fes_im)                   NA                     NA                    NA            NA            NA
SpeedLimit_im                   NA                     NA                    NA            NA            NA
strada_tip_im                   NA                     NA                    NA            NA            NA
log(P1_im)                      NA                     NA                    NA            NA            NA
log(E1_im)                      NA                     NA                    NA            NA            NA
                       strada_tip_im log(P1_im) log(E1_im)
(Intercept)                       NA         NA         NA
log(daylight_P50sp_im)            NA         NA         NA
log(daylight_Hits_im)             NA         NA         NA
log(d_fer_im)                     NA         NA         NA
SpeedLimit_im                     NA         NA         NA
strada_tip_im                     NA         NA         NA
log(P1_im)                        NA         NA         NA
log(E1_im)                        NA         NA         NA

Solution

  • Confidence intervals are based on (1) the parameter estimates and (2) the estimated standard errors (square roots of the variances) of these parameter estimates.

    In your example, the parameter estimates are okay; they are printed in the output as the "fitted trend coefficients" and could be extracted as a vector using coef(kppm_all_fes).

    It appears that the calculated standard errors are NA in your example. This could be because the calculated variances are NA, or are negative (because the square root of a negative number is returned as NA). You can check this by typing

    vcov(kppm_all_fes)
    

    which should return a square matrix (8 x 8 in your example) of non-negative finite numbers for the variances.

    There should have been some kind of warning message when the calculation was executed. This would also give information about what is going wrong.

    The variance estimates can be infinite or NA in some special cases. One of the special cases is a power law model, as explained in section 9.3.8 of Baddeley Rubak and Turner (2015). A power-law relationship is represented by the formula ~ log(Z) where Z is the original covariate. Your example model involves power-law relationships in some of the variables. So this could be the explanation.

    To solve your problem completely, I would need to have access to your data, or a Minimal Reproducible Example of the problem.

    EDIT: in this case, after getting access to the data, we were able to determine that the problem was a bug in vcov.kppm in its handling of NA values in the covariates. The bug has been fixed in spatstat.model 3.2-8.004. I will leave the original answer above, because it identifies the "legitimate" reasons for NA values in the standard errors.