JAGS模型末行报语法错误的原因排查求助

运行rjags时触发语法错误：

Error parsing model file: syntax error on line 115 near "#"

但第115行是模型末尾的注释行，无异常内容。自查括号配对未发现问题，附上完整模型代码请求协助定位错误原因。

原始模型代码 #

model{
#-----------------------------------------------------------------------------------------------------
# 1) Priors
#-----------------------------------------------------------------------------------------------------
  
  ######## Abundance ########
  mean_alpha0.psi ~ dnorm(0, 0.01)   # Mean intercept for occupancy probability for all species
  sd_alpha0 ~ dunif(0, 10)      # Standard deviation on mean abundance for all species
  prec_alpha0 <- 1 / (sd_alpha0 ^ 2)
  
  for (h in 1:nspecies) {
     alpha0.psi[h] ~ dnorm(mean_alpha0.psi, prec_alpha0) # Mean log abundance per point for species h
    
    # Uninformative priors on covariates
    alpha1[h] ~ dnorm(mu_alpha1, tau_alpha1) # Effect of fire frequency for each species
    alpha2[h] ~ dnorm(mu_alpha2, tau_alpha2) # Effect of ground cover for each species
    alpha3[h] ~ dnorm(mu_alpha3, tau_alpha3) # Effect of vegetation height for each species
    alpha4[h] ~ dnorm(mu_alpha4, tau_alpha4) # Effect of total basal area for each species
    alpha5[h] ~ dnorm(mu_alpha5, tau_alpha5) # Effect of % deciduous area for each species
  } # close h
  
  sd_eps ~ dexp(1) # Standard deviation for random property effect
  prec_eps <- 1 / (sd_eps ^ 2) # Precision for random property effect
  
  for (i in 1:nprops) {
    eps[i] ~ dnorm(0, prec_eps) # Random property effect
  } # close i
  
  mu_alpha1 ~ dnorm(0, 0.001)
  tau_alpha1 ~ dgamma(0.01, 0.01)
  
  mu_alpha2 ~ dnorm(0, 0.001)
  tau_alpha2 ~ dgamma(0.01, 0.01)
  
  mu_alpha3 ~ dnorm(0, 0.001)
  tau_alpha3 ~ dgamma(0.01, 0.01)
  
  mu_alpha4 ~ dnorm(0, 0.001)
  tau_alpha4 ~ dgamma(0.01, 0.01)
  
  mu_alpha5 ~ dnorm(0, 0.001)
  tau_alpha5 ~ dgamma(0.01, 0.01)
  
  ######## Availability #########
  for (h in 1:nspecies) {
    mean_phi0[h] ~ dnorm(0, 0.001)    # Mean availability for each species
    sd_phi0[h] ~ dexp(1)              # Standard deviation of availability for each species
    prec_phi0[h] <- 1 / (sd_phi0[h] ^ 2)
    
    for (i in 1:nprops) {
      for (j in 1:npoints_site[i]) {
        phi0[h, i, j] ~ dnorm(mean_phi0[h], prec_phi0[h]) # Mean availability for species h at property i and point j
      } # close j
    } # close i
  } # close h
  
  phi1 ~ dnorm(0, 0.001) # coefficient of day of year in availiabity linear model
  phi2 ~ dnorm(0, 0.001) # coefficient of day of year^2 in availiabity linear model
  phi3 ~ dnorm(0, 0.001) # coefficient of time after sunrise in availiabity linear model
  phi4 ~ dnorm(0, 0.001) # coefficient of temp in availiabity linear model
  phi5 ~ dnorm(0, 0.001) # coefficient of wind in availiabity linear model
  
  
  
#-----------------------------------------------------------------------------------------------------
# 2) Ecological process model
#-----------------------------------------------------------------------------------------------------

  for (h in 1:nspecies) {
    for (i in 1:nprops) {
      for (j in 1:npoints_site[i]) {
        z[h,i,j] ~ dbern(psi[h,i,j]) # True occupancy based on occupancy probability
        psi[h,i,j] <- 1 / (1 + exp(-lpsi.lim[h,i,j]))
        lpsi.lim[i] <- min(999, max(-999, lpsi[i]))
        lpsi[i] <- alpha0.psi[h] + alpha1[h] * fire[i, j] + alpha2[h] * ground[i, j] + alpha3[h] *
          height[i, j] + alpha4[h] * tba[i, j] + alpha5[h] * decid[i, j] + eps[i]
        
        
        

#-----------------------------------------------------------------------------------------------------
# 3) Observational process model
#-----------------------------------------------------------------------------------------------------

        ###### Availability with time removal ######
        for (v in 1:vh[i, j]) {
          # loop over each visit (v) to each site
          
          # Availability is a function of intercept + covariates
          
          y[h,i,j,v] ~ dbern(mu.p[h,i,j,v]) # Detection/non-detection
          mu.p[h,i,j,v] <- z[h,i,j] * pa[h,i,j,v]
          
          logit(p[h,i,j,v]) <- phi0[h, i, j] + phi1 * obs_date[i, j, v] + phi2 * obs_date[i, j, v] * obs_date[i, j, v] + phi3 * obs_time[i, j, v] + phi4 * obs_temp[i, j, v] + phi5 * obs_wind[i, j, v]
          
          for (z in 1:ntbins) {
            # Probability of being available in each time bin
            # pi_pa in time bin >1 conditional on not being detected before
            pi_pa[h, i, j, z, v] <- p[h, i, j, v] * pow(1 - p[h, i, j, v], (z - 1))
            pi_pa_normalized[h, i, j, z, v] <- pi_pa[h, i, j, z, v] / pa[h, i, j, v]
          } # close z
          
          # Probability ever available
          pa[h, i, j, v] <- sum(pi_pa[h, i, j, 1:ntbins, v])
          
          # Time data - when were detections mentally removed? NOT SURE ABOUT THIS LINE
          ytb[h, i, j, 1:ntbins, v] ~ dbern(pi_pa_normalized[h, i, j, 1:ntbins, v])
          
        } # close v
      } # close j
    } # close i
  } # close h
} # end model

错误定位与修正方案 #

以下是导致JAGS解析失败的核心问题：

1. 变量索引不匹配
   在生态过程模块中，lpsi.lim[i]和lpsi[i]仅使用了地块索引i，但后续的psi[h,i,j]包含物种h和生境点j维度，索引维度不匹配会导致JAGS解析混乱。需修正为：

   lpsi.lim[h,i,j] <- min(999, max(-999, lpsi[h,i,j]))
   lpsi[h,i,j] <- alpha0.psi[h] + alpha1[h] * fire[i, j] + alpha2[h] * ground[i, j] + alpha3[h] *
     height[i, j] + alpha4[h] * tba[i, j] + alpha5[h] * decid[i, j] + eps[i]
   

2. 循环变量与状态变量重名
   观测模块内层循环使用z作为循环变量，但外层已经定义了状态变量z[h,i,j]，JAGS不允许循环变量与已有变量重名，会触发语法解析错误。需将循环变量改为其他名称（如tb）：

   for (tb in 1:ntbins) {
     pi_pa[h, i, j, tb, v] <- p[h, i, j, v] * pow(1 - p[h, i, j, v], (tb - 1))
     pi_pa_normalized[h, i, j, tb, v] <- pi_pa[h, i, j, tb, v] / pa[h, i, j, v]
   } # close tb
   

3. 变量定义顺序错误
   pi_pa_normalized的计算依赖pa，原代码中pa的定义在pi_pa循环之后，但pi_pa_normalized的计算在pa定义之前，违反JAGS变量先定义后使用的规则。需调整顺序：

   for (tb in 1:ntbins) {
     pi_pa[h, i, j, tb, v] <- p[h, i, j, v] * pow(1 - p[h, i, j, v], (tb - 1))
   } # close tb
   # Probability ever available
   pa[h, i, j, v] <- sum(pi_pa[h, i, j, 1:ntbins, v])
   for (tb in 1:ntbins) {
     pi_pa_normalized[h, i, j, tb, v] <- pi_pa[h, i, j, tb, v] / pa[h, i, j, v]
   } # close tb
   

修正后的核心模块示例 #

#-----------------------------------------------------------------------------------------------------
# 2) Ecological process model
#-----------------------------------------------------------------------------------------------------

  for (h in 1:nspecies) {
    for (i in 1:nprops) {
      for (j in 1:npoints_site[i]) {
        z[h,i,j] ~ dbern(psi[h,i,j]) # True occupancy based on occupancy probability
        psi[h,i,j] <- 1 / (1 + exp(-lpsi.lim[h,i,j]))
        lpsi.lim[h,i,j] <- min(999, max(-999, lpsi[h,i,j]))
        lpsi[h,i,j] <- alpha0.psi[h] + alpha1[h] * fire[i, j] + alpha2[h] * ground[i, j] + alpha3[h] *
          height[i, j] + alpha4[h] * tba[i, j] + alpha5[h] * decid[i, j] + eps[i]
        
        
        

#-----------------------------------------------------------------------------------------------------
# 3) Observational process model
#-----------------------------------------------------------------------------------------------------

        ###### Availability with time removal ######
        for (v in 1:vh[i, j]) {
          y[h,i,j,v] ~ dbern(mu.p[h,i,j,v]) # Detection/non-detection
          mu.p[h,i,j,v] <- z[h,i,j] * pa[h,i,j,v]
          
          logit(p[h,i,j,v]) <- phi0[h, i, j] + phi1 * obs_date[i, j, v] + phi2 * obs_date[i, j, v]^2 + phi3 * obs_time[i, j, v] + phi4 * obs_temp[i, j, v] + phi5 * obs_wind[i, j, v]
          
          for (tb in 1:ntbins) {
            pi_pa[h, i, j, tb, v] <- p[h, i, j, v] * pow(1 - p[h, i, j, v], (tb - 1))
          } # close tb
          
          pa[h, i, j, v] <- sum(pi_pa[h, i, j, 1:ntbins, v])
          
          for (tb in 1:ntbins) {
            pi_pa_normalized[h, i, j, tb, v] <- pi_pa[h, i, j, tb, v] / pa[h, i, j, v]
          } # close tb
          
          ytb[h, i, j, 1:ntbins, v] ~ dbern(pi_pa_normalized[h, i, j, 1:ntbins, v])
          
        } # close v
      } # close j
    } # close i
  } # close h

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