分年龄组生存与恢复状态空间模型参数估计偏差求助

我使用状态空间模型（State-Space Model）模拟估计分年龄组的恢复概率(p)和生存概率(phi)，但得到的参数估计结果与设定值偏差较大。

估计值： #

beta [1] = 0.331
beta [2] = 0.634 
delta [1] = 0.469
delta [2] = 0.529

模拟设定的真实值： #

beta [1] = 0.8
beta [2] = 0.8
delta [1] = 0.2
delta [2] = 0.1

以下是我的代码：

### Finding first capture year
get.first <- function(x)min(which(x!=0))

known.state.mr <- function (mr) {
  state <- matrix (NA, nrow=dim(mr)[1], ncol= dim(mr)[2])
  rec <- which(rowSums(mr)==2)
  for (i in 1: length(rec)){
    n1<- min(which(mr[rec[i],]==1)) #n1 time of capture
    n2<- max(which(mr[rec[i],]==1)) # Recovery time
    state[rec[i],n1:n2]<-1 
    state[rec[i],n1]<-NA 
    state[rec[i],n2:dim(mr)[2]]<- 0 
  }
  return(state)
}


mr.init.z<-function(mr){
  ch<-matrix(NA, nrow=dim(mr)[1],ncol=dim(mr)[2])
  rec<- which(rowSums(mr)==1) 
  for(i in 1:length(rec)){ 
    n1<-which(mr[rec[i],]==1)
    ch[rec[i], n1:dim(mr)[2]]<-0
    ch[rec[i],n1] <- NA 
  }
  return (ch)
}

# Define function to simulate mark-recovery data
simul.mr <- function(S, R, marked){
  n.occasions <- dim(S)[2]
  MR <- matrix(NA, ncol = n.occasions+1, nrow = sum(marked))
  # Define a vector with the occasion of marking
  mark.occ <- rep(1:n.occasions, marked)
  # Fill the CH matrix
  for (i in 1:sum(marked)){
    MR[i, mark.occ[i]] <- 1    
    for (t in mark.occ[i]:n.occasions){
      # Bernoulli trial: has individual survived occasion? 
      sur <- rbinom(1, 1, S[i,t])
      if (sur==1) next     
       
      rp <- rbinom(1, 1, R[i,t])
      if (rp==0){
        MR[i,t+1] <- 0
        break
      }
      if (rp==1){
        MR[i,t+1] <- 1
        break
      }
    } #t
  } #i
  # Replace the NA in the file by 0
  MR[which(is.na(MR))] <- 0
  return(MR)
}


########### Simulate Data #########################

# Juvenile Matrix
n.occasions <- 14                 # Number of release occasions
marked <- rep(50, n.occasions)    # Annual number of newly marked individuals
as <- rep(0.8, n.occasions)
ar <- rep(0.2, n.occasions)

# Define matrices with survival and recovery probabilities
aS <- matrix(as, ncol = n.occasions, nrow = sum(marked))
aR <- matrix(ar, ncol = n.occasions, nrow = sum(marked))


# Execute function
MR.A <- simul.mr(aS, aR, marked)


### Juvenile Matrix #####
js <- rep(0.8, n.occasions)
jr <- rep(0.1, n.occasions)

# Define matrices with survival and recovery probabilities
jS <- matrix(js, ncol = n.occasions, nrow = sum(marked))
jR <- matrix(jr, ncol = n.occasions, nrow = sum(marked))

# Execute function
MR.J <- simul.mr(jS, jR, marked)

get.first <- function(x) min(which(x!=0))
f.j <- apply(MR.J, 1, get.first)
f.a <- apply(MR.A, 1, get.first)

# Create matrices X indicating age classes
x.j <- matrix(NA, ncol = dim(MR.J)[2]-1, nrow = dim(MR.J)[1])
x.a <- matrix(NA, ncol = dim(MR.A)[2]-1, nrow = dim(MR.A)[1])

for (i in 1:dim(MR.J)[1]){
  for (t in f.j[i]:(dim(MR.J)[2]-1)){
    x.j[i,t] <- 2
    x.j[i,f.j[i]] <- 1   
  } #t
} #i

for (i in 1:dim(MR.A)[1]){
  for (t in f.a[i]:(dim(MR.A)[2]-1)){
    x.a[i,t] <- 2
  } #t
} #i

MR <- rbind(MR.J, MR.A)
f <- c(f.j, f.a)
x <- rbind(x.j, x.a)


# Specify model in BUGS language
sink("Model_MR")
cat("
model {
# Priors and constraints
for (i in 1:nind){
   for (t in f[i]:(n.occasions-1)){
      
      phi[i,t] <- beta[x[i,t]]
      p[i,t] <- delta[x[i,t]]
      
      } #t
   } #i
   
for (u in 1:2){
   
   beta[u] ~ dunif(0,1) # Priors for age-specific survival
   delta[u] ~ dunif(0,1)
   
}
   
mean.p ~ dunif(0, 1)                  # Prior for mean recapture
# Likelihood 
for (i in 1:nind){
   
   # Define latent state at first capture
   
   z[i,f[i]] <- 1
   for (t in (f[i]+1):n.occasions){
      
      # State process
      
      z[i,t] ~ dbern(mu1[i,t])
      mu1[i,t] <- phi[i,t-1] * z[i,t-1]
      
      # Observation process
      
      y[i,t] ~ dbern(mu2[i,t])
      mu2[i,t] <- p[i,t-1] * (z[i,t-1]-z[i,t])
      
      } #t
   } #i
}
",fill = TRUE)
sink()

# Bundle data
bugs.data <- list(y = MR, 
                  f = f, 
                  nind = dim(MR)[1],
                  n.occasions = dim(MR)[2],
                  z = known.state.mr(MR), 
                  x = x)

# Initial values
inits <- function(){list(z = mr.init.z(MR), 
                         beta = runif(2, 0, 1), 
                         mean.p = runif(1, 0, 1))}  

# Parameters monitored
parameters <- c("beta", "mean.p","delta")

# MCMC settings
ni <- 2000
nt <- 3
nb <- 1000
nc <- 3

# Call WinBUGS from R (BRT 3 min)
cjs.age <- jags(bugs.data, inits, parameters, "Model_MR", 
n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb)

print(cjs.age, digits = 3)

问题排查方向： #

1. 初始值函数逻辑错误
   mr.init.z函数中，对仅被标记无恢复记录的个体，直接将标记后所有时刻的状态z设为0，但这类个体可能存活到研究结束，初始值设为0会误导MCMC收敛，应改为将标记后到最后时刻的z初始化为1。

2. 年龄组矩阵赋值冗余
   构建x.j时，循环内每次都重新赋值x.j[i,f.j[i]] <- 1，虽然最终结果正确，但逻辑冗余，可调整为仅在t=f.j[i]时赋值1，其余时刻赋值2：

   for (i in 1:dim(MR.J)[1]){
     x.j[i,f.j[i]] <- 1
     if(f.j[i] < dim(MR.J)[2]-1){
       x.j[i,(f.j[i]+1):(dim(MR.J)[2]-1)] <- 2
     }
   }
   

3. BUGS模型时间索引不匹配
   模拟数据中，MR的列数是n.occasions+1，对应标记时刻1到恢复时刻n.occasions+1；但模型中n.occasions设为dim(MR)[2]，观测方程里t对应恢复时刻，而phi[i,t-1]和p[i,t-1]对应t-1时刻的生存/恢复概率，需确认时间步是否与模拟逻辑完全对应。

4. MCMC迭代次数不足
   当前设置ni=2000，nb=1000，有效样本量可能不足，导致估计值不稳定。建议增加迭代次数至ni=10000，nb=5000，并检查Rhat值是否接近1，确认模型收敛。

5. 冗余参数干扰
   模型中定义了mean.p但未使用，应删除该参数的先验定义，减少MCMC不必要的计算。

内容的提问来源于stack exchange，提问作者ariggs154