Random effect of year not estimated well

[Cassidy Waldrep]

Hello, 

I'm fitting a model in NIMBLE and though I got it to converge, I'm worried that the year random effect is not being estimated well. 

My lab mate ran the exact same model in brms and while the sd of state was exactly the same, the sd of year in the brms model was much smaller and had a tighter estimate. The estimate of the sd of year in the nimble model is around 2.3 whereas in the brms model it's 0.45. Again, the sd of state in both models is the same. 

I'm wondering if anyone sees an issue in this code, specifically with year as a random effect. Is there a reason I'm not estimating year well? There are four years in my model. The general idea of the model is to look at how bander, weight, temperature, and day all influence daily survival after banding mallards. 

Here's my model: 

model_code <- nimbleCode({
 
  ## priors -------------------------------------------------------------------
 
  ### Beta coefficients  ------------------------------------------------------

 
  for(b in 1:nbander) {
    beta_bander[b] ~ dnorm(0, sd = 1.5)
  }

  beta_mean_mintemp_prev5 ~ dnorm(0, sd = 1.5) # CLW added new parameter
  beta_weight ~ dnorm(0, sd = 1.5) # CLW added new parameter
 
  for(d in 1:nday) {
    beta_day[d] ~ dnorm(0, sd = 1.5)
  }
 
  ### Random Effects ---------------------------------------------------------

  for(y in 1:nyear) {
    eps_year[y] ~ dnorm(mu_year, sd = sd_year)
  }
 
  mu_year ~ dnorm(0, 1.5)
  sd_year ~ dexp(1) # Prior for standard deviation of years

  for(s in 1:nstate) {
    eps_state[s] ~ dnorm(mu_state, sd = sd_state)
  }
 
  mu_state ~ dnorm(0, 1.5)
  sd_state ~ dexp(1)
 
  ## likelihoods --------------------------------------------------------------
  # figuring out indexing for day
 
  for (m in 1:nobs) {
   
    # random part
   
    survival[m] ~ dbern(p[m])
   
    # link function
   
    logit(p[m]) <- beta_bander[bander_index[m]] +
      beta_day[day_index[m]] +
      beta_mean_mintemp_prev5*mean_mintemp_prev5[m] +
      beta_weight*weight[m] +
      eps_year[year_index[m]] +
      eps_state[state_index[m]]
  }
})


# set up NIMBLE model ---------------------------------------------------------

# constants
nimble_constants <- list(
  nobs = nobs,
  nyear = nyear,
  nday = nday,
  nbander = nbander,
  nstate = nstate,
  bander_index = data_model$bander_index,
  year_index = data_model$year_index,
  day_index = data_model$day_index,
  state_index = data_model$state_index,
  mean_mintemp_prev5 = as.numeric(data_model$mean_mintemp_prev5_scaled),
  weight = as.numeric(data_model$weight_scaled)
)

# data
nimble_data <- list(
  survival = data_model$survived
)


# run NIMBLE model in series ---------------------------
nc <- 3

# nimble_inits <- inits_function()

# build model
model <- nimbleModel(code = model_code,
                     constants = nimble_constants,  
                     dat =  nimble_data,
                     buildDerivs = T)

# configure MCMC
mcmc_Conf  <- configureHMC(model, print = F)

# build MCMC
modelMCMC  <- buildMCMC(mcmc_Conf)

# compile model and MCMC
Cmodel     <- compileNimble(model)
CmodelMCMC <- compileNimble(modelMCMC, project = model)

# run model
samples <- runMCMC(CmodelMCMC,
                   niter = 5000,
                   nburnin = 1000,
                   thin = 10,
                   nchains = nc,
                   samplesAsCodaMCMC = T)

Let me know if any questions.

Best, 

Cassidy 

[Perry de Valpine]

Hi Cassidy,

Thanks for the question. From looking through your code, I see that some of the dnorm priors have sd=1.5 but others have just 1.5, but without naming the parameter as "sd", the default will be precision (tau), and you'll end up with sd = 1/sqrt(1.5) = 0.82. If you intended those priors to be sd=1.5, does fixing that resolve the discrepancy between the posteriors from nimble and Stan?

If not, I'd suggest sending also your code for brms so the two model specifications can be compared more directly. Are you confident other components of the priors match between the two?

HTH

Perry

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/nimble-users/313ddcf7-e9c2-490f-a099-20aad28173b4n%40googlegroups.com.