1

I am using spatstat v.1.59-0 to build some point process models, but I am having problems with some of the validation tools, specifically effectfun and and the residual plot parres

The components of the model are rings_pp which consists of point locations for 52 archaeological sites, and elev which is a DEM converted to a spatstat pixel image. I am trying to evaluate the fit of a Gibbs point process model with an Area Interaction parameter (AreaInter) and elev

I fit the model to the data using the following code:
rr1 <- data.frame(r=seq(100, 2000, by=50)) p1 <- profilepl(rr1, AreaInter, rings_pp~elev, aic=T) ppm5 <- as.ppm(p1)

Everything seems to be working fine (e.g., other diagnostics show the model is a reasonable fit to the data) except when I try to use effectfun and parres to evaluate the effect of the elevation covariate, they wont work.

parres gives the error Error in m[, d$relevant, drop = FALSE] : (subscript) logical subscript too long When I traceback the error I get the following:

`> trace(parres(ppm5, "elev"))
Error in m[, d$relevant, drop = FALSE] : 
  (subscript) logical subscript too long
> traceback()
4: effectFun.can(x)
3: effectFun(xxx)
2: parres(ppm5, "elev")
1: trace(parres(ppm5, "elev"))`

effectfun works when se.fit=FALSE but will not work when se.fit=TRUE and gives the following error: Error in quadform(mm, vc) : ncol(x) == nrow(v) is not TRUE When I traceback the error I get the following:

`> trace(plot(effectfun(ppm5, "elev", se.fit=T)))
Error in quadform(mm, vc) : ncol(x) == nrow(v) is not TRUE
> traceback()
8: stop(simpleError(msg, call = sys.call(-1)))
7: stopifnot(ncol(x) == nrow(v))
6: quadform(mm, vc)
5: predict.ppm(orig.model, locations = fakeloc, covariates = fakecov, 
   se = se.fit)
4: predict(orig.model, locations = fakeloc, covariates = fakecov, 
   se = se.fit)
3: effectfun(ppm5, "elev", se.fit = T)
2: plot(effectfun(ppm5, "elev", se.fit = T))
1: trace(plot(effectfun(ppm5, "elev", se.fit = T)))`

When I fit an inhomogeneous model ppm1 <- ppm(rings_pp, ~elev) both effectfun and parres work fine and show a good fit (though using the residual K function Kres suggest the model isn't accounting for clustering). So it seems to be something with how I've fit the Gibbs AreaInter model.

Any advice would be much appreciated.

dinapoli
  • 11
  • 1
  • Thanks for a good question with many details. Can you either supply your data or try to recreate the problem with a simulated dataset or using one of the built in datasets? It's much easier to debug a problem which can be reproduced locally. – Ege Rubak Apr 25 '19 at 06:37
  • Hi Dr. Rubak. Thank you for your reply. I was unable to recreate the problem with one of the example datasets (bei), and I don't think it's possible to upload my data to stackoverflow. The DEM covariate is pretty large (380 mb). However, I would be happy to email you the data and the R script if that is okay with you. – dinapoli Apr 25 '19 at 16:26
  • If at all possible post a public link to e.g. Dropbox for anybody to download and test the example. I will not have time to run the code and debug it the next week or more, but I will probably get around to it sometime. I cannot receive 380 MB in an email. I think my limit is 10 MB. – Ege Rubak Apr 25 '19 at 22:26

1 Answers1

0

Thanks again for a clear question.

This is a bug, which can be demonstrated by the simple example

fit <- ppm(cells ~ x, AreaInter(0.07))
plot(effectfun(fit, se.fit=TRUE))
plot(parres(fit))

I will investigate and fix it soon.

Postscript: this has now been fixed in the development version of spatstat (version 1.59-0.020) available on the GitHub repository

.

Adrian Baddeley
  • 2,534
  • 1
  • 5
  • 8
  • Hi Dr. Baddeley. Okay thanks so much! Your R book with Rubak and Turner is fantastic. Any plans for a 2nd edition? – dinapoli Apr 26 '19 at 14:56