error in cpo recompute inla.cpo mac m1

[Jonathan Simkin]

Hi there,

I'm running a bym2 model. I get no problem in terms of running the model, getting the summary and plotting the model output. A few months ago (May) I could easily recompute the CPOs with inla.cpo but I tried re-running my code and suddenly the inla.cpo crashes.

I tried re-installing RStudio, R, and stable and testing versions of INLA... not sure if this is an issue with my local set up.

Thanks!

J

Here's my code and error follows:

> prior <- list(
  prec = list(
    prior = "pc.prec",
    param = c(0.2 / 0.31, 0.01),
    initial = 5),
  phi = list(
    prior = "pc",
    param = c(0.5, 2 / 3),
    initial = -3)
)

> res_bym2 <- inla(formula_bym,
            family = "poisson", data = chsadf_inla,
            E = exp, control.predictor = list(compute = TRUE),
            control.compute = list(dic=TRUE, mlik=TRUE,cpo=TRUE, config = T, waic = T),
            control.inla = list(strategy = "laplace", npoints = 21))

Here's the error I get

> improved_res <- inla.cpo(res_bym2)
Compute new CPO/PIT values manually, for 30 cases...
Error in xx$cpo : $ operator is invalid for atomic vectors
In addition: Warning message:
In parallel::mclapply(..., mc.cores = mc.cores) :
  all scheduled cores encountered errors in user code

I also tried with another model result

ompute new CPO/PIT values manually, for 88 cases...


    GitID: file: error-handler.c  16bac7359467ef92a59050760b3d6da476e1846d - Thu Jul 21 18:20:46 2022 +0200
    Error:21 Reason: This should not happen
    Message: Condition `density->Pinv' is not TRUE
    Line:912 Function: GMRFLib_density_Pinv

sh: line 1: 71009 Abort trap: 6           '/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/INLA/bin/mac.arm64/inla.run' -s -v -t1:1 -B0 -P classic '/private/var/folders/jf/tsgz8r717w1br8rb664b19d00000gn/T/RtmpNkmipx/file765d924a824/Model.ini' > '/private/var/folders/jf/tsgz8r717w1br8rb664b19d00000gn/T/RtmpNkmipx/file765d924a824/Logfile.txt'

and this error ^ iterates for a while...

My session info:

R version 4.2.1 (2022-06-23)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Monterey 12.1

Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib

locale:
[1] en_CA.UTF-8/en_CA.UTF-8/en_CA.UTF-8/C/en_CA.UTF-8/en_CA.UTF-8

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

other attached packages:
 [1] viridis_0.6.2      viridisLite_0.4.0  RColorBrewer_1.1-3 tmaptools_3.1-1    gridExtra_2.3      hrbrthemes_0.8.0  
 [7] INLA_22.07.21-2    sp_1.5-0           foreach_1.5.2      Matrix_1.4-1       tmap_3.3-3         magrittr_2.0.3    
[13] sf_1.0-7           forcats_0.5.1      stringr_1.4.0      dplyr_1.0.9        purrr_0.3.4        readr_2.1.2      
[19] tidyr_1.2.0        tibble_3.1.7       ggplot2_3.3.6      tidyverse_1.3.1  

loaded via a namespace (and not attached):
 [1] fs_1.5.2           lubridate_1.8.0    httr_1.4.3         tools_4.2.1        backports_1.4.1    utf8_1.2.2        
 [7] R6_2.5.1           KernSmooth_2.23-20 DBI_1.1.2          colorspace_2.0-3   raster_3.5-9       withr_2.5.0      
[13] tidyselect_1.1.2   leaflet_2.1.1      compiler_4.2.1     extrafontdb_1.0    leafem_0.2.0       cli_3.3.0        
[19] rvest_1.0.2        xml2_1.3.3         scales_1.2.0       classInt_0.4-3     proxy_0.4-26       systemfonts_1.0.4
[25] digest_0.6.29      rmarkdown_2.14     base64enc_0.1-3    dichromat_2.0-0.1  pkgconfig_2.0.3    htmltools_0.5.2  
[31] extrafont_0.18     dbplyr_2.1.1       fastmap_1.1.0      htmlwidgets_1.5.4  rlang_1.0.2        readxl_1.4.0      
[37] rstudioapi_0.13    generics_0.1.2     jsonlite_1.8.0     crosstalk_1.2.0    s2_1.0.7           Rcpp_1.0.8.3      
[43] munsell_0.5.0      fansi_1.0.3        gdtools_0.2.4      abind_1.4-5        lifecycle_1.0.1    terra_1.4-22      
[49] stringi_1.7.6      leafsync_0.1.0     crayon_1.5.1       lattice_0.20-45    stars_0.5-5        haven_2.5.0      
[55] splines_4.2.1      hms_1.1.1          knitr_1.39         pillar_1.7.0       codetools_0.2-18   wk_0.6.0          
[61] reprex_2.0.1       XML_3.99-0.9       glue_1.6.2         evaluate_0.15      modelr_0.1.8       png_0.1-7        
[67] vctrs_0.4.1        tzdb_0.3.0         MatrixModels_0.5-0 Rttf2pt1_1.3.10    cellranger_1.1.0   gtable_0.3.0      
[73] assertthat_0.2.1   xfun_0.31          lwgeom_0.2-8       broom_0.8.0        e1071_1.7-9        class_7.3-20      
[79] iterators_1.0.14   units_0.8-0        ellipsis_0.3.2    

[Helpdesk]

to check, I would need data++ so I can rerun it here.


however, there is work in progress on a new (group-)cpo implementation
(paper coming on arxiv soon), that do cpo as default. the approximations
are done differently, and more stable and more accurate.

its only for experimental mode, like this:

can you check if this would work better for you?

> r=inla(y ~ 1, data =data.frame(y=rnorm(10)),
control.compute=list(control.gcpo=list(enable=TRUE)),
inla.mode="experimental")
> r$gcpo$gcpo
[1] 0.509624652604 0.366757093028 0.323252764109 0.003874241164
0.483853663926
[6] 0.553997740777 0.461789820581 0.236727841468 0.488620518390
0.420246564313

> --
> You received this message because you are subscribed to the Google
> Groups "R-inla discussion group" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to r-inla-discussion...@googlegroups.com.
> To view this discussion on the web, visit
> https://groups.google.com/d/msgid/r-inla-discussion-group/1d8a395f-adf8-4944-aa1a-c28841259ff1n%40googlegroups.com
> .

--
Håvard Rue
he...@r-inla.org

[Jonathan Simkin]

My small area data is sensitive but I will tweak it and attach a copy soon with the graph structure for bym2. Thanks!

I ran the gcpo and that works! I also ran the same model with cpo=T (not in experimental mode), and inla.cpo works. (output below). 

> r<-inla(y ~ 1, data =data.frame(y=rnorm(10)),
+ control.compute=list(control.gcpo=list(enable=TRUE)),
+ inla.mode="experimental")
> r$gcpo$gcpo
 [1] 0.26262786 0.28616346 0.28385820 0.22408453 0.27350315 0.08868468 0.17671519 0.05909643 0.02096582 0.21480555
>
> r2<-inla(y ~ 1, data =data.frame(y=rnorm(10)),
+ control.compute=list(
+   cpo = T
+   )
+ )
>
> r2$cpo$failure[1] = 1      
> r2$cpo$cpo
 [1] 0.22542982 0.25743535 0.14716585 0.16791578 0.31492531 0.36915148 0.22855609 0.29569158 0.38738110 0.03415789
> inla.cpo(r2)
Compute new CPO/PIT values manually, for 1 cases...
     index cpo.old  cpo.new  pit.old  pit.new
[1,]     1 0.22543 0.225427 0.154653 0.154618

The retured result contain the new values.
Call:
   c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, quantiles =
   quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, Ntrials = Ntrials, strata =
   strata, ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = verbose, ", " lincomb =
   lincomb, selection = selection, control.compute = control.compute, ", " control.predictor =
   control.predictor, control.family = control.family, ", " control.inla = control.inla, control.fixed =
   control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " control.hazard =
   control.hazard, control.lincomb = control.lincomb, ", " control.update = control.update, control.lp.scale
   = control.lp.scale, ", " control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, ", "
   inla.call = inla.call, inla.arg = inla.arg, num.threads = num.threads, ", " blas.num.threads =
   blas.num.threads, keep = keep, working.directory = working.directory, ", " silent = silent, inla.mode =
   inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = .parent.frame)")
Time used:
  Pre = 1.33, Running = 0.0908, Post = 0.00631, Total = 1.43
> 

[Jonathan Simkin]

Actually I could reproduce the error with another dataset, a public scotland lip cancer dataset. I've attached the data (scotlip_shiny_input.csv), the R script (scot_test.R) and the shapefiles (scotlip_test dbf, prj, shp, shx) but you can also just use the graph structure scot.adj.

I get the same error here as well working with BYM2, same session info as before. 

    GitID: file: error-handler.c  16bac7359467ef92a59050760b3d6da476e1846d - Thu Jul 21 18:20:46 2022 +0200
    Error:21 Reason: This should not happen
    Message: Condition `density->Pinv' is not TRUE
    Line:912 Function: GMRFLib_density_Pinv

sh: line 1: 75581 Abort trap: 6           '/Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/library/INLA/bin/mac.arm64/inla.run' -s -v -t1:1 -B0 -P classic '/private/var/folders/jf/tsgz8r717w1br8rb664b19d00000gn/T/RtmpfpF8vW/filed5a130986f87/Model.ini' > '/private/var/folders/jf/tsgz8r717w1br8rb664b19d00000gn/T/RtmpfpF8vW/filed5a130986f87/Logfile.txt'

Error in xx$cpo : $ operator is invalid for atomic vectors
In addition: Warning message:
In parallel::mclapply(..., mc.cores = mc.cores) :
  all scheduled cores encountered errors in user code

[Helpdesk]

+ read.csv("scotlip_shiny_input.csv")

> map <-
+ st_read("scot_map") %>% st_make_valid()
Error: Cannot open "scot_map"; The file doesn't seem to exist.

On Fri, 2022-07-22 at 13:35 -0700, Jonathan Simkin wrote:
> Actually I could reproduce the error with another dataset, a public
> scotland lip cancer dataset. I've attached the data
> (scotlip_shiny_input.csv), the R script (scot_test.R) and the
> shapefiles (scotlip_test dbf, prj, shp, shx) but you can also just use
> the graph structure scot.adj.
>
> I get the same error here as well working with BYM2, same session info
> as before. 
>
>     GitID: file: error-handler.c
>  16bac7359467ef92a59050760b3d6da476e1846d - Thu Jul 21 18:20:46 2022
> +0200
>     Error:21 Reason: This should not happen
>     Message: Condition `density->Pinv' is not TRUE
>     Line:912 Function: GMRFLib_density_Pinv
>
> sh: line 1: 75581 Abort trap: 6          

> '/Library/Frameworks/R.framework/Versions/4.2-
> arm64/Resources/library/INLA/bin/mac.arm64/inla.run' -s -v -t1:1 -B0 -
> P classic

> '/private/var/folders/jf/tsgz8r717w1br8rb664b19d00000gn/T/RtmpfpF8vW/f
> iled5a130986f87/Model.ini' >
> '/private/var/folders/jf/tsgz8r717w1br8rb664b19d00000gn/T/RtmpfpF8vW/f
> iled5a130986f87/Logfile.txt'

[Jonathan Simkin]

oh apologies for that - that's because I was reading it in from a folder called "scot_map" that contained the shapefile... 

I think you can ignore the lines that create the graph structure since I've  attached the graph (scot.adj). So that would ignore:

lines10-11

map <-

  st_read("scot_map") %>% st_make_valid()

lines 25-26
set.ZeroPolicyOption(TRUE)scot
get.ZeroPolicyOption()

lines 28-29
nb <- poly2nb(map)
nb2INLA(file = "scot.adj", nb)

Otherwise, you could throw the shapefiles into a folder and use st_read() to the folder path. Let me know if you want me to re-write the script!

J

[Helpdesk]

it runs into overflow due to large uncertainty...

I would recommend using the default strategy (just remove the 'laplace'
line), and use


res <- inla(formula,

family = "poisson", data = chsadf_inla,
E = exp, control.predictor = list(compute = TRUE),

control.compute = list(dic = TRUE, cpo = T, config = T,
return.marginals.predictor=TRUE))

plot(res)
inla.cpo(res)



we have worked on a new approach for removing more than one points, and
then the computations are more stable for the cpo (which is a special
case), so I guess in the future, we'll redefine 'cpo' to be

res2 <- inla(formula,

family = "poisson", data = chsadf_inla,
E = exp, control.predictor = list(compute = TRUE),

inla.mode = "experimental",
control.compute = list(control.gcpo = list(enable = TRUE),
return.marginals.predictor=TRUE))

and then 'res2$gcpo$gcpo' is the cpo values (which is the default).



plot(res$cpo$cpo, res2$gcpo$gcpo)

> > > > > > mpNk
> > > > > > mipx/f
> > > > > > ile765d924a824/Model.ini' >
> > > > > > '/private/var/folders/jf/tsgz8r717w1br8rb664b19d00000gn/T/Rt

[Jonathan Simkin]

Fantastic! I can confirm it's working well :) Thanks so much. 

Looking forward to the implementation of gcpo - will that also have failure check and inla.cpo recomputation.

Appreciate all your help!!

Jonathan