Модель пространства состояний для временных рядов населения – общие опасения и соответствие

Question

Впервые постер, давний читатель. Я также совершенно новичок в байесовском анализе, поэтому, пожалуйста, будьте со мной осторожны.

Я пытаюсь подогнать модель к наблюдениям за взрослыми организмами, которые являются семивольтинными, то есть им требуется два года, чтобы достичь взрослой жизни, а затем они умирают вскоре после спаривания.

Поскольку меня беспокоит ошибка наблюдения, я пытаюсь подогнать байесовскую модель пространства состояний. Мои вопросы таковы:

• Кажется ли это несколько разумным подходом? Мне особенно любопытно, не обидится ли кто-нибудь на использование отрицательного биномиального распределения ошибок в модели процесса или распределения Пуассона в модели наблюдения.

• Когда я запускаю предоставленный код и смотрю на апостериорные оценки количества взрослых, оценки почти полностью совпадают с данными, что заставляет меня думать, что я делаю что-то очень неправильно. Однако оценки параметров аналогичны тем, которые я получаю, когда запускаю GLM без ошибки наблюдения.

      # 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)])