poissonconsulting / mcmcr

An R package to manipulate MCMC samples
https://poissonconsulting.github.io/mcmcr/
Other
17 stars 2 forks source link

Rhat function returns unexpected error #28

Closed EcologyTom closed 5 years ago

EcologyTom commented 5 years ago

I have run a model using the jagsUI package, with the function jags.basic. I have a list of mcmc objects as output;

List of 8
  $ : 'mcmc' num [1:80, 1:2227] 76.2 75.6 75.7 75.5 75.5 ...
   ..- attr(*, "dimnames")=List of 2
   .. ..$ : NULL
   .. ..$ : chr [1:2227] "SS1[1]" "SS1[2]" "SS1[3]" "SS1[4]" ...
   ..- attr(*, "mcpar")= num [1:3] 206 601 5
  $ : 'mcmc' num [1:80, 1:2227] 76.3 75.9 75.7 75.7 75.7 ...
   ..- attr(*, "dimnames")=List of 2
   .. .. etc etc

I now want to use the rhat function, from my reading of the manual I ought to be able to do so directly on the lists of mcmc objects? However, I run into the following error;

Error in array(data = as.vector(x), dim = c(1, niters(x), pdims)) :    
negative length vectors are not allowed 
In addition: There were 33 warnings (use warnings() to see them)
--
Warning messages: 
1: In lapply(x, as.integer) : NAs introduced by coercion
2: In lapply(x, as.integer) : NAs introduced by coercion 
etc (all 33 are the same)

Here is the output from my model run; mcmc_test_output.zip

Perhaps I need to convert the mcmc list into another class first? However, I get the same error when I try as.mcmcr. Maybe I am just doing something really daft? If you can shed any light on this I would be really grateful. Many thanks.

Session Info ```r R version 3.6.0 (2019-04-26) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows >= 8 x64 (build 9200) Matrix products: default locale: [1] LC_COLLATE=English_United Kingdom.1252 [2] LC_CTYPE=English_United Kingdom.1252 [3] LC_MONETARY=English_United Kingdom.1252 [4] LC_NUMERIC=C [5] LC_TIME=English_United Kingdom.1252 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] mcmcr_0.2.0 tidyr_0.8.3 dplyr_0.8.3 [mcmc_test_output.zip](https://github.com/poissonconsulting/mcmcr/files/3443010/mcmc_test_output.zip) loaded via a namespace (and not attached): [1] Rcpp_1.0.1 lattice_0.20-38 crayon_1.3.4 assertthat_0.2.1 [5] grid_3.6.0 R6_2.4.0 magrittr_1.5 coda_0.19-3 [9] checkr_0.5.0 pillar_1.4.2 rlang_0.4.0 rjags_4-8 [13] err_0.2.0 tools_3.6.0 glue_1.3.1 purrr_0.3.2 [17] jagsUI_1.5.0 abind_1.4-5 parallel_3.6.0 compiler_3.6.0 [21] pkgconfig_2.0.2 tidyselect_0.2.5 tibble_2.1.3 --   > | > > ```
joethorley commented 5 years ago

The problem was "." in the term names. I've updated the dev version so it now handles these.

mcmc <- readRDS("mcmc_test_output.rds")
rhat(mcmc)
# [1] 64.534

If you want to use the current CRAN release then you'll need to remove the "."s in the parameter names.

Note I'm afraid it's really slow to calculate the rhat on an object this big. I've filed an issue to optimize but I don't have time right now. FWIW coda::gelman.diag() errors on the object.

coda::gelman.diag(mcmc)
# Error in chol.default(W) : 
#  the leading minor of order 474 is not positive definite

I'm also attaching a csv of the rhat values by term in case its useful.

df <- rhat(mcmc, by = "term", as_df = TRUE)
write.csv(df, "rhat.csv")

rhat.csv.zip

EcologyTom commented 5 years ago

Great, thanks a lot. I should probably avoid using "." in my object names anyway, and this is timely reminder.