2D Integration points with region & domain producing NA indicies

[jdselwyn]

When trying to create integration points using ipoints when I include both the domain and region for the integration points I get the following error:

Error in intI(j, n = x@Dim[2], dn[[2]], give.dn = FALSE) : 
  'NA' indices are not (yet?) supported for sparse Matrices

I don't get the error with the gorillas dataset so it seems to be something about the combination of my domain and region but I can't tell what. It appears that the specific error message is triggered by this line A <- proj_cover$proj$A[, mesh_cover$idx$loc, drop = FALSE] within the function integration_weight_construction called internally. I can't figure out why but it looks like two of the values in mesh_cover$idx$loc are NA for unknown reasons.

I've attached a zip file with an rds file for both the domain (mesh.rds) and region (site.rds).

library(inlabru)

mesh <- readRDS("mesh.rds")
site <- readRDS("site.rds")

integration_points <- ipoints(region = site, domain = mesh)

#Both of the below produce an output
integration_points <- ipoints(domain = mesh)
integration_points <- ipoints(region = site)

This also happens when I try to use the samplers argument in lgcp.

Thank you for any help/advice you have. mesh_and_site.zip

[finnlindgren]

Thanks for the report! Can you also provide your sessionInfo() output?

[jdselwyn]

The sessionInfo() for the computer that I used to create the code above is here:

> sessionInfo()
R version 4.0.2 (2020-06-22)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19041)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252    LC_MONETARY=English_United States.1252 LC_NUMERIC=C                           LC_TIME=English_United States.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] inlabru_2.1.13 ggplot2_3.3.2  sp_1.4-2      

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.5            INLA_19.09.03         compiler_4.0.2        pillar_1.4.6          tools_4.0.2           rpart_4.1-15          goftest_1.2-2         lifecycle_0.2.0       tibble_3.0.3         
[10] nlme_3.1-148          gtable_0.3.0          lattice_0.20-41       mgcv_1.8-31           pkgconfig_2.0.3       rlang_0.4.7           Matrix_1.2-18         rstudioapi_0.11       parallel_4.0.2       
[19] rgdal_1.5-16          terra_0.8-11          spatstat.data_1.4-3   withr_2.2.0           dplyr_1.0.2           raster_3.3-16         generics_0.0.2        vctrs_0.3.4           spatstat.utils_1.17-0
[28] tidyselect_1.1.0      grid_4.0.2            glue_1.4.2            R6_2.4.1              polyclip_1.10-0       purrr_0.3.4           deldir_0.1-25         tensor_1.5            magrittr_1.5         
[37] scales_1.1.1          codetools_0.2-16      ellipsis_0.3.1        MASS_7.3-51.6         splines_4.0.2         abind_1.4-5           spatstat_1.64-1       colorspace_1.4-1      munsell_0.5.0        
[46] crayon_1.3.4  

I also get the error with a computer with this sessionInfo()

> sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux Server release 6.10 (Santiago)

Matrix products: default
BLAS/LAPACK: /home/apps/openblas/0.3.3/lib/libopenblas_sandybridgep-r0.3.3.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
 [9] LC_ADDRESS=C               LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods
[8] base

other attached packages:
 [1] doParallel_1.0.15    iterators_1.0.12     foreach_1.5.0
 [4] patchwork_1.0.1      rgeos_0.5-3          INLA_99.99.9999
 [7] Matrix_1.2-18        inlabru_2.1.13       sf_0.9-6
[10] magrittr_1.5         forcats_0.5.0        stringr_1.4.0
[13] dplyr_1.0.2          purrr_0.3.4          readr_1.3.1
[16] tidyr_1.1.1          tibble_3.0.3         ggplot2_3.3.2
[19] tidyverse_1.3.0.9000 terra_0.8-6          raster_3.3-13
[22] sp_1.4-2

loaded via a namespace (and not attached):
 [1] httr_1.4.2            jsonlite_1.7.0        splines_3.5.1
 [4] modelr_0.1.8          assertthat_0.2.1      blob_1.2.1
 [7] cellranger_1.1.0      pillar_1.4.6          backports_1.1.9
[10] lattice_0.20-41       glue_1.4.1            polyclip_1.10-0
[13] rvest_0.3.6           colorspace_1.4-1      pkgconfig_2.0.3
[16] broom_0.7.0           haven_2.3.1           scales_1.1.1
[19] tensor_1.5            spatstat.utils_1.17-0 mgcv_1.8-32
[22] generics_0.0.2        ellipsis_0.3.1        withr_2.2.0
[25] cli_2.0.2             deldir_0.1-28         crayon_1.3.4
[28] readxl_1.3.1          fs_1.5.0              fansi_0.4.1
[31] nlme_3.1-149          MASS_7.3-52           xml2_1.3.2
[34] class_7.3-17          tools_3.5.1           hms_0.5.3
[37] lifecycle_0.2.0       munsell_0.5.0         reprex_0.3.0
[40] compiler_3.5.1        e1071_1.7-3           rlang_0.4.7
[43] classInt_0.4-3        units_0.6-0           grid_3.5.1
[46] rstudioapi_0.11       goftest_1.2-2         gtable_0.3.0
[49] codetools_0.2-16      abind_1.4-5           DBI_1.1.0
[52] R6_2.4.1              splancs_2.01-40       lubridate_1.7.9
[55] rgdal_1.5-18          KernSmooth_2.23-17    stringi_1.4.6
[58] spatstat.data_1.4-3   spatstat_1.64-1       Rcpp_1.0.5
[61] vctrs_0.3.2           rpart_4.1-15          dbplyr_1.4.4
[64] tidyselect_1.1.0

I believe the second session showing INLA_99.99.9999 is related to some issues with installing it on RedHat6 (but that's well beyond my understanding). It was installed ~Feb 2019 and appears to run without errors that I've noticed.

[finnlindgren]

This appears to be a very strange bug in the fmesher code itself in the INLA package. I wonder if it might be related to the code merging etc that had to be done when moving INLA from a mercurial repo on bitbucket to a git repo on github. I have a separate codebase for fmesher to compare with so should be easy to fix in that case; otherwise I'll need to do some more fmesher debugging, and/or modify the method in inlabru to work around it. Edit: if your INLA is from 2019, it's not a code merge problem, so I'll go ahead and do fmesher debugging.

[finnlindgren]

I believe I've solved it in the inlabru devel branch by replacing some code with simpler code that doesn't involve creating the specific extra internal meshes that triggered the fmesher problem. The devel branch is going through a major update, so if you try it you may need to tweak some code a bit. The main thing that you'll notice is that effectname(map = covariate, ...) should now be either effectname(covariate, ... or effectname(main = covariate, ...) in the component specifications. The map= argument should however still work for backwards compatibility, but will give you a warning.

[jdselwyn]

Thanks for looking into this and solving the issue!

Unfortunately I got a new error Error: 'inla.sp_get_crs' is not an exported object from 'namespace:INLA'

This is the new sessionInfo() after updating packages to install the development version of inlabru

R version 4.0.2 (2020-06-22)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19041)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252    LC_MONETARY=English_United States.1252 LC_NUMERIC=C                           LC_TIME=English_United States.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] inlabru_2.1.14.900 sp_1.4-4           ggplot2_3.3.2     

loaded via a namespace (and not attached):
 [1] INLA_19.09.03         compiler_4.0.2        pillar_1.4.6          tools_4.0.2           rpart_4.1-15          goftest_1.2-2         lifecycle_0.2.0       tibble_3.0.4          gtable_0.3.0         
[10] nlme_3.1-148          lattice_0.20-41       mgcv_1.8-31           pkgconfig_2.0.3       rlang_0.4.8           Matrix_1.2-18         rstudioapi_0.11       spatstat.data_1.4-3   withr_2.3.0          
[19] dplyr_1.0.2           generics_0.0.2        vctrs_0.3.4           spatstat.utils_1.17-0 grid_4.0.2            tidyselect_1.1.0      glue_1.4.2            R6_2.5.0              polyclip_1.10-0      
[28] purrr_0.3.4           deldir_0.1-25         tensor_1.5            magrittr_1.5          scales_1.1.1          ellipsis_0.3.1        MASS_7.3-51.6         splines_4.0.2         abind_1.4-5          
[37] spatstat_1.64-1       colorspace_1.4-1      munsell_0.5.0         crayon_1.3.4

[finnlindgren]

Due to recent changes in how rgdal and sp handle CRS information (PROJ6 instead of PROJ4), the new inlabru version needs a newer inla version (20.06.18 or later), so on your R4.0.02 Windows system you should be ok by upgrading INLA to a recent/current version. On the older system you'd likely need to go through a similar procedure to what you mentioned to build a compatible version. But I believe there are some tricks to bypass that that might work as an alternative (see INLA::inla.binary.install()) to building on some systems.

But the amount of code in inlabru that needs these helper function is relatively small, so I'm checking if I can duplicate those in inlabru, at least temporarily, to make it possible to use older versions of INLA in cases where the PROJ4 information still works.

[jdselwyn]

Thank you! updating to the testing version of INLA solved the issue with not finding the function from INLA.

Thank you for looking into duplicating the functions in inlabru. I believe I should be able to update INLA on the older machine, thanks for pointing out the binary installation

[finnlindgren]

The latest devel version now has its own functions for dealing with CRS information, so it might work with older INLA versions.

[finnlindgren]

Small note: when crs information is available for the mesh and data, inlabru constructs integration weights for "per km squared" units. Your case only has around 250 square metres, so the integration weights are really small, and in some cases I think it may lead to floating point approximation issues; in the future, I plan to add an option to allow the user to override the units of integration. Currently, that can only be done by setting the crs information on the mesh to CRS(NA_character_) which leads it to use raw coordinates for everything.

[jdselwyn]

Ok thank you. I had wondered if something like that was happening because I had an lgcp model which had previously been able to fit reasonably well the total number of points ~300 but was now giving estimates of the total number much less than 1. If I set the the CRS to NA when the CRS for the region and coordinates originally had units in "m" then the output will be "per m squared" units if I'm interpreting using the raw coordinates right?

[finnlindgren]

Yes, that's correct. The estimated intensity should differ by a factor 10^6. And to be clear, the "per km^2" and "CRS vs no CRS" is an old feature, so should be the same in old and current inlabru versions. The initial inlabru version was geared towards point processes covering large parts of the Pacific ocean, so km^2 was the natural choice, since we had to take into account the curvature of the earth and not use non-equal-area projections for the calculations. The current version transforms the domain to a lambert equal area projection when constructing the integration weights, when crs information is available. This conversion seems to contributed to the original error in the fmesher code (in a way I have yet to figure out; the changed code bypasses that problem).