Spline in JAGS mixing badly

288 Views Asked by At

I have a model that calculates a spline for Mark-recapture data with survival data. The model is working fine, but the parameters that calculates the spline are mixing super badly.

                        mean        2.5%      97.5%   That n.eff
...
m[1]             1.667899656   -0.555606    4.18479 2.8829     4
m[2]             1.293023680   -0.951046    3.90294 2.8476     4
m[3]             1.717855378   -0.484097    4.23105 2.8690     4
m[4]             1.723899423   -0.474260    4.23869 2.8686     4
m[5]             1.747050770   -0.456455    4.26314 2.8578     4
...

Basically, I'm calculating a recapture rate p composed of a species specific effect p.sp and the sampling effort p.effort. I also calculate a fitness component phi with a species specific term phi.sp, the effect of year phi.year, a climate factor phi.sum.preci and the spline m.

run.model <- function(d, ## incoming data (packaged up in src/analyses.R)
                      ni=1100, ## number of iterations to run ## number of draws per chain
                      nt=10, ## thinning rate ##to save space on computer disk space see p.61 Kéry
                      nb=100, ## burn in ## should be large enough to discard initial part of Markov chains that have not yet converged
                      nc=3, ## number of chains to run ## multiple chain to check the convergence 
                      n.cluster = 3) { 

  model.jags <- function() {
    ## Priors ------------------------------------------------------------------
    ## Random effect species-specific intercept (survival)
    mu.phi.sp ~ dnorm(0,0.01)
    sigma.phi.sp ~ dunif(0,10)
    tau.phi.sp <- 1/(sigma.phi.sp)^2
    ## Random effect for recapture rate 
    mu.p.sp ~ dnorm(0,0.01)

    ## Random effect of year and fixed effect of precipitation & abundance
    sigma.phi.year ~ dunif(0,10)
    tau.phi.year <- 1/(sigma.phi.year)^2

    ## fixed effect of effort
    p.effort ~ dnorm(0, 0.01) ## fixed effect

    ## Fixed precipitation per year 
    phi.sum.preci ~ dnorm(0, 0.01) ## fixed effect

    # Prior spline ------------------------------------------------------------
    ###BEGIN SPLINE### 
    # prior distribution for the fixed effects parameters
    for (l in 1:3) {
      beta[l] ~ dnorm(0,0.1)
    }
    prior.scaleeps <- 1
    xi ~ dnorm(0, tau.xi)
    tau.xi <- pow(prior.scaleeps, -2)

    for (k in 1:nknotsb) {
      b[k] <- yi*etab[k]
      etab[k] ~ dnorm(0, tau.etab) # hierarchical model for theta
    } # closing k

    prior.scaleb <- 1
    yi ~ dnorm (0, tau.yi)
    tau.yi <- pow(prior.scaleb, -2)
    tau.etab ~ dgamma(.5, .5) # chi^2 with 1 d.f.
    sigmab <- abs(xi)/sqrt(tau.etab) # cauchy = normal/sqrt(chi^2)
    ###END SPLINE###

    for(sp in 1:nsp) {
      ## Random species-specific intercept
      phi.sp[sp]  ~ dnorm(mu.phi.sp, tau.phi.sp)
      ## Random recapture rate
      p.sp[sp] <- mu.p.sp # Changed from a comment from Luke Jan. 9 2017
    }

    for (yr in 1:nyear) {
      ## random year
      phi.year[yr] ~ dnorm(0, tau.phi.year)
    }

    ## Likelihood!
    for(sp in 1:nsp) { ## per species
      ## Rates -------------------------------------------------------------------
      ## recapture rate
      for (yr in 1:nyear) {
        logit(p[sp,yr]) <- # added logit here
          p.sp[sp]  +
          p.effort*effort[yr]
      } ## closing for (year in 1:nyear)
    } ## closing for (sp in 1:nsp)

    ## Each ID  ----------------------------------------------------------------    
    ## Likelihood!
    for(ind in 1:nind) { ## nind = nrow(d$X)
      ### BEGIN SPLINE ###
      ## mean function model
      m[ind] <-mfe[ind] + mre1[ind] + mre2[ind]

      # fixed effect part
      mfe[ind] <- beta[1] * Xfix[ind,1] +beta[2] * Xfix[ind,2] + beta[3] * Xfix[ind,3]

      mre1[ind] <- b[1]*Z[ind,1] + b[2]*Z[ind,2] + b[3]*Z[ind,3] + b[4]*Z[ind,4] + b[5]*Z[ind,5] + b[6]*Z[ind,6] + b[7]*Z[ind,7] + b[8]*Z[ind,8] + b[9]*Z[ind,9] + b[10]*Z[ind,10]
      mre2[ind] <- b[11]*Z[ind,11] + b[12]*Z[ind,12] + b[13]*Z[ind,13] + b[14]*Z[ind,14] + b[15]*Z[ind,15]
      ###END SPLINE###
    }

    ## for each individual
    for(ind in 1:nind) { ## nind = nrow(d$X)
      for(yr in 1:nyear) {
        logit(phi[ind,yr]) <-
          phi.sp[species[ind]] + ## effect of species 
          phi.year[yr] + ## effect of year
          # Effect of the traits on survival values
          m[ind]+ # spline
          phi.sum.preci*sum.rainfall[yr]  # effect of precipitation per sampling event
      } ## (yr in 1:nyear)

      ## First occasion
      for(yr in 1:first[ind]) { 
        z[ind,yr] ~ dbern(1)
      } ## (yr in 1:first[ind])

      ## Subsequent occasions
      for(yr in (first[ind]+1):nyear) { # (so, here, we're just indexing from year "first+1" onwards).
        mu.z[ind,yr] <- phi[ind,yr-1]*z[ind,yr-1]
        z[ind,yr] ~ dbern(mu.z[ind,yr])
        ## Observation process
        sight.p[ind,yr] <- z[ind,yr]*p[species[ind],yr] ## sightp probability of something to be seen 
        X[ind,yr] ~ dbern(sight.p[ind,yr]) ## X matrix : ind by years
      } ## yr
    } ## closing for(ind in 1:nind)
  } ## closing model.jags function 

  ## Calling Jags ------------------------------------------------------------    
  jags.parallel(data = d$data,
                inits = d$inits,
                parameters.to.save = d$params,
                model.file = model.jags,
                n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb, 
                working.directory = NULL,
                n.cluster = n.cluster) 
} ## closing the run.model function

# Monitored parameters ----------------------------------------------------
get.params <- function() 
  c('phi.sp','mu.phi.sp','sigma.phi.sp','mu.p.sp','sigma.p.sp','phi.year','phi','p', 'phi.sum.preci','p.sp','p.effort','z',
    # Spline parameters
    "m","sigmab","b","beta")
0

There are 0 best solutions below