State space model for population time series - general worries and fit

Question

First time poster, long time reader. I am also totally new to Bayesian analyses, so please be gentle with me.

I am trying to fit a model to observations of adult organisms that are semivoltine--that is, they take two years to reach adulthood and then die shortly after mating.

Because I worry about the observation error, I'm trying to fit a Bayesian state space model. My questions are these:

• Does this seem like a somewhat reasonable approach? I'm particularly curious if anyone takes offense to the use of negative binomial error distribution in the process model---or the poisson distribution in the observation model.

• When I run the provided code and look at the posterior estimates of adult counts, the estimates almost entirely match the data, which makes me think that I'm doing something very wrong. However, the parameter estimates are similar to those I get when i run a GLM without the observation error.

# Load Packages
library(R2jags)
library(rjags)

#Data
Counts <- c(156, 62, 62, 300, 415, 429, 366, 188, 142, 201, 375, 
  438, 508, 386, 105, 293, 167, 105, 146, 129, 153, 220, 
  162, 74, 94, 202, 217, 246, 120, 182)

Data <- c(list(TYr = length(Counts),
               N1 = Counts[1],
               N2 = Counts[2],
               Counts = Counts))

#State space Ricker model
mod<- function(){
  #Process model. 
  for (t in 3:TYr){ #data# TYr
    N[t] ~ dnegbin(p[t],size) 
    p[t]<- size/(size+lambda[t])
    log(lambda[t])<- log(N[t-2]) + a + b*N[t-2] 
  }
  
  #Observation model
  for (i in 1:TYr){ 
    Counts[i] ~ dpois(N[i]) #data# Counts
  }
  
  #Priors
  N[1] = (N1)  #data# N1
  N[2] = (N2)  #data# N2
  a~dnorm(0,0.001)
  b~dnorm(0,0.0001)
  size ~ dunif(0.001, 50)
  
  #Derived values, I'd like to be able to report the estimated carrying capacity, K
  K = -exp(a)/b
}

#Initial values. b has to start negative or I run into issues
Ricker_inits_1<-list(a = 0.1, b = -0.0001)
Ricker_inits_2<-list(a = 0.15, b = -0.001)
Ricker_inits_3<-list(a = 0.3, b = -0.01)
Total_Ricker_inits<-list(Ricker_inits_1,Ricker_inits_2,Ricker_inits_3)

#fit the model
jags_fit <- jags(data = Data,
                     inits = Total_Ricker_inits,
                     parameters.to.save = c('a','b',"K", 'size'),
                     model.file = mod,
                     n.chains=3,
                     n.burnin=3000,
                     n.iter = 20000,
                     n.thin = 1)

#get posterior estimates
samples <- jags.samples(jags_fit$model, 
                           c("N"), 
                           type = "mean", 
                           n.iter = 20000,
                           n.burnin = 1000,
                           n.thin = 1)

#crazy high goodness of fit makes me question everything 
#removing the first two points because these were explicit in the model
plot(Data$Counts[-c(1,2)]~summary(samples$N, FUN = mean)$stat[-c(1,2)])
abline(0,1)

#similar parameter estimates for a and b
require(MASS)
glm.nb(Counts[-c(1,2)]~offset(log(Counts[-c(29,30)]))+Counts[-c(29,30)])