mikemeredith / ahm_code Goto Github PK

Code from the book 'Applied Hierarchical Modeling in Ecology'

R 100.00%

ahm_code's Introduction

AHM code

The two-volume work Applied Hierarchical Modeling in Ecology: Analysis of distribution, abundance and species richness in R and BUGS by Marc Kéry and Andy Royle contains lots of R and BUGS code.

The R package AHMbook, available on CRAN, has all the data sets and the custom functions used in the books. Commented code for the functions is on GitHub here.

This repository has all the code in the printed books, plus code referred to as "available on the website" but not printed. The aim is to have code which works with current versions of R, JAGS and contributed R packages. The code is regularly tested and updated code inserted, with the original printed code retained but commented out with #. Please open an issue if you find other code which does not work.

In addition to these updates, some code has been inserted:

Code added at the top of the script to recreate or reload objects from previous sections; each script is self-contained.
Some changes have been made to facilitate automated checking of scripts, in particular reductions in the number of iterations for simulations, bootstraps and MCMC runs.
After long runs of unmarked, JAGS or WinBUGS, I've inserted code to save the results to RData files.

Additional code and comments are marked off with twiddly lines like this:

#~~~~ oldfunction has been replaced with newfunction ~~~~~~~
# oldfunction(foo)
newfunction(foo)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Avoiding WinBUGS

WinBUGS is not essential to work through the code. In most cases, jagsUI::jags is a drop-in replacement for R2WinBUGS::bugs. JAGS does not do spatial autocorrelation (CAR) models, as used in AHM2 chapters 3 and 9; for those, the nimble package can be used, and alternative scripts are provided.

Volume 1 (AHM1)

The code for the first volume (AHM1), with updates up to 2017, is available as a single huge text file on the main book web page.

Volume 2 (AHM2)

The book appeared in October 2020 (with copyright dated 2021). The code here is based on the final proofs.

ahm_code's People

Contributors

Stargazers

Watchers

Forkers

gpandradep primulalhd anhnguyendepocen rhine3 smandujanor elahi jjherlan ecustwy trashbirdecology bmourgui kls229 ecsflsf jeff-wagner moniquepm memo1986 xuletajr farcego aldridgecaleb

ahm_code's Issues

NIMBLE model for AHM2 section 9.5 does not converge

The code currently implemented in the script "AHM2_09.5_nimble.R" does not converge, as can be seen in the attached diagnostic plots.

(The plots are for chains of 6000 iterations thinned by 6, after 6000 iterations burn-in.)

The corresponding WinBUGS model does converge, though mixing is poor, especially for the variance of the CAR component.

~~Two of the chains are stuck, presumably at their starting values, and the third is moving very slowly. ...~~

After some more poking around, it seems that the green chain has gone rogue, and the black and red chains have converged pretty well. We can proceed by removing the rogue chain from the mcmc.list output and using just the rest. We may be able to avoid the problem by judicious selection of RNG seeds.

calculation of chi2 in GOF test for DYNOCC models

Dear Mike,

First of all, thanks a lot for sharing the codes from the AHM-book. It has been incredibly helpful for my work.
I have recently tried to do a model check of a dynamic occupancy model based on your code AHM_code/AHM2_ch04/AHM2_04.08.R. On line 197-205 in the code, chi2 is calculated. I apologize if there is something I'm missing, but it seems to me like the chi2 is estimated as: (observed - expected)^2 / observed, instead of the typical way of estimating chi2 :(observed - expected)^2 / expected. Is it a reason for this or am I missing something?

Best Regards,
Eivind

non-ASCII characters in comments in BUGS/JAGS code

I've just hit a weird problem when loading saved R workspace images. Example:

Error in load("AHM2_03.4.4.Rdata") : 
  embedded nul in string: '  # Note b\0\031open indexb\0\031 in pr[t,s,] comes last'

The culprit seems to be non-ASCII characters - specifically smart quotes - in a comment in the BUGS/JAGS code. That code should
go straight to the model file as text (via cat) and should not be in the R workspace. The jagUI output object includes the compiled JAGS model, which does include the model text, but this happens with R2WinBUGS output too.

Replacing the smart quotes with straight quotes solves the problem.

Error message for updated JagsUI/Coda MCMC object specification: Error in mcmc.list(x) : Arguments must be mcmc objects

Thank you for reporting a bug in our scripts.

The code here is based on the book, and we update it to work with current versions of JAGS, R and R packages. Please make sure you have the latest version of the script and software (run update.packages()) before reporting an issue.

File name: AHM1_ch11

Lines which produced the error or unexpected output (copy/paste from the console):
op <- par(mfrow = c(4,4)) ; traceplot(out5) ; print(out5, dig = 3)

error message: Error in mcmc.list(x) : Arguments must be mcmc objects

Do you have a suggested update for the code:

Coda specifies that in the updated version, older mcmc objects will result in an error message. "n previous releases of CODA, an mcmc object could be a single or multiple chains. A new class mcmc.list has now been introduced to deal with multiple chains and mcmc objects can only have data from a single chain."

mcmcUpgrade(x) can be used to update outdated code for mcmc objects; or objects could be modified to specified correctly under the updated version of the coda package

Run sessionInfo() and paste the output here:

R version 4.0.3 (2020-10-10)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.6

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

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

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

other attached packages:
[1] AHMbook_0.2.0 coda_0.19-4 jagsUI_1.5.1 lattice_0.20-41

loaded via a namespace (and not attached):
[1] Rcpp_1.0.5 mvtnorm_1.1-1 codetools_0.2-16 unmarked_1.0.1
[5] MASS_7.3-53 grid_4.0.3 plyr_1.8.6 RandomFields_3.3.8
[9] raster_3.3-13 sp_1.4-4 rjags_4-10 RColorBrewer_1.1-2
[13] tools_4.0.3 plotrix_3.7-8 parallel_4.0.3 compiler_4.0.3
[17] RandomFieldsUtils_0.5.3

Other information that may be helpful to pinpoint the issue, eg, screenshots:

In which package can one find the function 'whiskerplot' and 'denstrip' ?

First, warm thanks for your excellent opus. Being more field ecologist than statistician, it gives me an opportunity to gradually understand (and move ?) from the frequencist concepts to Bayesians (a wish I wanted to fulfil since long but I was bogged down by various kind of bureaucracies before retiring ;-) ).

My question is very simple and not exactly a relate to a bug : pages 157-158 volume 1, you use the functions 'whiskerbox' and 'denstrip'. I have not been able to find whiskerbox in any package. A denstrip does exist in the package 'denstrip'.

Can you tell me where to find whiskerbox and if the denstrip you use is the one of the package 'denstrip' ?

Best,

Patrick

Object not found for 11.7.2

Hi, in section 11.7.2, I am running into this error when trying to enter this code:

str(win.data <- list(y = yaug, nsite = dim(y)[1], nrep = dim(y)[2],
nspec = dim(y)[3], nz = nz, M = nspec + nz, ele = ele, forest = forest,
DAT = DAT, DUR = DUR) )

error message:

Error in str(win.data <- list(y = yaug, nsite = dim(y)[1], nrep = dim(y)[2], :
object 'DAT' not found

I have the AHM book package loaded, I'm not sure if the object name was changed. When I look at the MHB2014 dataset, I see "date" and "dur" listed as objects, but not DAT and DUR. I actually tried substituting these for DAT and DUR but received the error message:

Error in str(win.data <- list(y = yaug, nsite = dim(y)[1], nrep = dim(y)[2], :
object 'dur' not found

I'm not sure what I'm missing, any help appreciated! Thanks!!