Error in Holmes-Held update method

  • Jacob Socolar

    Jacob Socolar - 2014-06-13

    I'm still a relative neophyte to both JAGS and this forum. I am trying to fit a large multi-species occupancy model, running JAGS from R using R2jags.

    After the model graph is compiled, JAGS consistently crashes and throws the error

    Invalid precision in Holmes-Held update method.
    This is a known bug and we are working on it.
    Please bear with us

    Is there any available work-around? For example, is there a way to tell JAGS to use the Albert-Chib sampler instead of Holmes-Held, and would this even solve the problem? Alternatively, is there a hope that I could solve the problem with inconsequential modifications to the code (pasted below) or with different choices of initial values?

    I'm perfectly happy to run the model using rjags instead of R2jags if that makes any difference. I just don't know whether it's possible in either rjags or R2jags to specify that I don't want JAGS to use the Holmes-Held update method.

    My code is pasted after my signature, in case that's of any help. I apologize if the code itself is amateurish--like I said, I am a neophyte here. The model estimates the true occupancy probability for 404 bird species at 256 sites, and also estimates the probability that the bird participates in mixed species flock at a site, given that it occupies the site.

    Best, and thanks for JAGS--it's wonderful.

    Jacob Socolar
    PhD Candidate
    Dept of Ecology and Evolutionary Biology
    Princeton University


    # This is a multi-species occupancy model, to run on the data contained in X and Fl, two arrays generated by Format_for_occupancy.R
    # Priors    
    hm.alpha.occ ~ dunif(-100,100)
     ## Mean intercept for occupancy (logit scale)
    hv.alpha.occ ~ dunif(0,100)
    hm.beta1.occ ~ dunif(-100,100)
     ## Mean linear effect of elevation (logit scale)
    hv.beta1.occ ~ dunif(0,100)
    hm.beta2.occ ~ dunif(-100,100)
     ## Mean quadratic effect of elevation (logit scale)
    hv.beta2.occ ~ dunif(0,100)    
    hm.alpha.p ~ dunif(-100,100)
     ## Intercept for detection (logit scale)
    hv.alpha.p ~ dunif(0,100)
     ## Variance
    hm.beta6.p ~ dunif(-100,100)
     ## Mean linear effect of elevation
    hv.beta6.p ~ dunif(0,100)
     ## Variance
    hm.alpha.flocc ~ dunif(-100,100)
    ## Mean intercept for floccupancy
    hv.alpha.flocc ~ dunif(0,100)
    ## Variance
    hm.beta1.flocc ~ dunif(-100,100)
    ## Mean effect of elevation
    hv.beta1.flocc ~ dunif(0,100)
    ## Variance
    hm.beta2.flocc ~ dunif(-100,100)
    ## Mean effect of elevation squared
    hv.beta2.flocc ~ dunif(0,100)
    ## Variance
    hm.alpha.fp ~ dunif(-100,100)
    ## mean intercept for detection in flock given detection
    hv.alpha.fp ~ dunif(0,100)
    ## variance for detection in flock
    hm.beta6.fp ~ dunif(-100,100)
    ## mean linear effect of elevation
    hv.beta6.fp ~ dunif(0,100)
    ## variance in linear effect of elevation
        # Model structure
    for(h in 1:404){
    alpha.occ[h] ~ dnorm(hm.alpha.occ, hv.alpha.occ)
    beta1.occ[h] ~ dnorm(hm.beta1.occ, hv.beta1.occ)
    beta2.occ[h] ~ dnorm(hm.beta2.occ, hv.beta2.occ)
    alpha.p[h] ~ dnorm(hm.alpha.p, hv.alpha.p)
    beta6.p[h] ~ dnorm(hm.beta6.p, hv.beta6.p)
    alpha.flocc[h] ~ dnorm(hm.alpha.flocc, hv.alpha.flocc)
    beta1.flocc[h] ~ dnorm(hm.beta1.flocc, hv.beta1.flocc)
    beta2.flocc[h] ~ dnorm(hm.beta2.flocc, hv.beta2.flocc)
    alpha.fp[h] ~ dnorm(hm.alpha.fp, hv.alpha.fp)
    beta6.fp[h] ~ dnorm(hm.beta6.fp, hv.beta6.fp)
        for(i in 1:R){
        # R is the number of points, see code near line 43: <- ...
        # occ[h,i] holds the true occupancy status of each site i for species h
        # lp.occ[h,i] holds the linear predictor for the occupancy status of each site
        logit(lp.occ[h,i]) <- alpha.occ[h] + beta1.occ[h]*elev1[i] + beta2.occ[h]*elev2[i]
        # Note that the covariate must be a funciton of i
        flocc[h,i] ~ dbern(occ[h,i]*lp.flocc[h,i])
        logit(lp.flocc[h,i]) <- alpha.flocc[h] + beta1.flocc[h]*elev1[i] + beta2.flocc[h]*elev2[i]
        for(j in 1:T){
        # T is the max number of visits per point, see code near line 43: <- ...
        X[i,j,h] ~ dbern(eff.p[h,i,j])
        # We're working here with the hth species.
        # The actual detection history given by the true occupancy state times the probability of detection occ[h,i]*p[h,i,j]
        eff.p[h,i,j] <- occ[h,i]*p[h,i,j]
        # eff.p is the true probability of detection on the jth visit to site i
        # p[h,i,j] gives the probability of detection given occupancy on the jth visit to site i
        logit(p[h,i,j]) <- alpha.p[h] + beta6.p[h]*elev1[i]
        # note that the covariates may be functions of i, j, or both
        Fl[i,j,h] ~ dbern(f.eff.p[h,i,j])
        f.eff.p[h,i,j] <- flocc[h,i]*flocc.p[h,i,j]*X[i,j,h]
        logit(flocc.p[h,i,j]) <- alpha.fp[h] + beta6.fp[h]*elev1[i]
        Presi[h,i,j] <- abs(X[i,j,h]-p[h,i,j])
        # Presi[h,i,j] gives the absolute residual[h,i,j] ~ dbern(eff.p[h,i,j])[h,i,j] <- abs([h,i,j]-p[h,i,j])
        fit <- sum(Presi[,,]) <- sum([,,])
    #    occ.fs <- sum(occ[,])
        ", fill=TRUE)
    sink() <- list(X=X, Fl=Fl, elev1=as.numeric(covariates[,1,29]), elev2=as.numeric(covariates[,1,30]), hps1=covariates[,,31], hps2=covariates[,,32], obs=covariates[,,35], R=dim(X)[1], T=dim(X)[2])
    zst <- mat.or.vec(404,256)
    fzst <- mat.or.vec(404,256)
    for(h in 1:404){
      zst[h,] <- apply(X[,,h],1,max,na.rm=T)
    fzst[h,] <- apply(Fl[,,h],1,max,na.rm=T)
    inits <- function(){list(occ=zst, flocc=fzst, hm.alpha.occ=runif(1,-10,10), hv.alpha.occ=runif(1,0,10), alpha.occ=runif(404,-1,1), hm.beta1.occ=runif(1,-10,10), hv.beta1.occ=runif(1,0,10), beta1.occ=runif(404,-1,1), hm.beta2.occ=runif(1,-10,10), hv.beta2.occ=runif(1,0,10), beta2.occ=runif(404,-1,1), hm.alpha.p=runif(1,-10,10), hv.alpha.p=runif(1,0,10), alpha.p=runif(404,-1,1), hm.beta6.p=runif(1,-10,10), hv.beta6.p=runif(1,0,10), beta6.p=runif(404,-1,1), hm.alpha.flocc=runif(1,-10,10), hv.alpha.flocc=runif(1,0,10), alpha.flocc=runif(404,-1,1), hm.beta1.flocc=runif(1,-10,10), hv.beta1.flocc=runif(1,0,10), beta1.flocc=runif(404,-1,1), hm.alpha.fp=runif(1,-10,10), hv.alpha.fp=runif(1,0,10), alpha.fp=runif(404,-1,1), hm.beta6.fp=runif(1,-10,10), hv.beta6.fp=runif(1,0,10), beta6.fp=runif(404,-1,1))}
    params <- c("hm.alpha.occ", "hm.beta1.occ", "hm.beta2.occ", "hm.alpha.p", "hm.beta6.p", "hm.alpha.flocc", "hm.beta1.flocc", "hm.beta2.flocc", "hm.alpha.fp", "hm.beta6.fp", "hv.alpha.occ", "hv.beta1.occ", "hv.beta2.occ", "hv.alpha.p", "hv.beta6.p", "hv.alpha.flocc", "hv.beta1.flocc", "hv.beta2.flocc", "hv.alpha.fp", "hv.beta6.fp", "alpha.occ[27]", "beta1.occ[27]", "beta2.occ[27]", "alpha.p[27]", "beta6.p[27]", "alpha.flocc[27]", "beta1.flocc[27]", "beta2.flocc[27]", "alpha.fp[27]", "beta6.fp[27]", "fit", "")
    nc <- 3
    nb <- 5000
    ni <- 20000
    nt <- 20
    MS_flock_model <- jags(, inits, params, "model.txt", n.chains=nc, n.iter=ni, n.burnin=nb, n.thin=nt,
    Last edit: Jacob Socolar 2014-06-13
  • Martyn Plummer

    Martyn Plummer - 2014-06-13

    I fixed this in JAGS 3.4.0 by adding a small regularization term that was just big enough to overcome these numerical problems without perturbing the posterior.

    Upgrade if you can. If you are on a Mac then I appreciate that a binary for 3.4.0 is still not avaiable, so you can turn off the Holmes-Held factory with this:

    set.factory("glm::Holmes-Held", type="sampler", state=FALSE)

    Check the result:

    > list.factories("sampler")
                      factory status
    1         glm::ConjugateF   TRUE
    2        glm::Holmes-Held  FALSE
    3        glm::Albert-Chib   TRUE
    4  glm::Albert-Chib-Gibbs   TRUE
    5  glm::Auxiliary-Mixture   TRUE
    6             glm::Linear   TRUE
    7        glm::LinearGibbs   TRUE
    8               glm::IWLS   TRUE
    9              bugs::DSum   TRUE
    10        bugs::Conjugate   TRUE
    11        bugs::Dirichlet   TRUE
    12          bugs::MNormal   TRUE
    13           base::Finite   TRUE
    14            base::Slice   TRUE
  • Jacob Socolar

    Jacob Socolar - 2014-06-13

    Thanks so much, both for JAGS and for your phenomenal helpfulness on this site.

    Yes, I am on a Mac, but the set.factory command seems to have worked properly.

    I should know for sure in another couple of hours once the model is compiled and running that JAGS is behaving as expected (i.e. not attempting Holmes-Held sampling).

    To your knowledge, does there exist a book/article/document, accessible to a programming dilettante, that motivates a rough understanding of what these different update methods are, and also what "adaptation" does? I'd love a vague sketch of the under-the-hood workings of JAGS so that I can justify (mostly to myself) my acceptance of the default preferences.

  • Martyn Plummer

    Martyn Plummer - 2014-06-16

    I try to explain what adaptation is in the user manual. For the various methods, you need to look in the source code, where you will find the references. There is no overview document, although I probably should write one.


Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:

No, thanks