OpenBUGS错误未定义变量 [英] OpenBUGS error undefined variable

问题描述

我正在使用OpenBUGS和R包R2OpenBUGS构建二项式混合模型.我已经成功构建了更简单的模型，但是一旦我添加了另一个用于不完善检测的级别，我就会持续收到错误variable X is not defined in model or in data set.我尝试了许多不同的方法，包括更改数据结构并将数据直接输入OpenBUGS.我发布此消息的目的是希望其他人对此错误有经验，并且也许知道为什么OpenBUGS不能识别变量X，即使据我所知它是明确定义的.

I'm working on a binomial mixture model using OpenBUGS and R package R2OpenBUGS. I've successfully built simpler models, but once I add another level for imperfect detection, I consistently receive the error variable X is not defined in model or in data set. I've tried a number of different things, including changing the structure of my data and entering my data directly into OpenBUGS. I'm posting this in the hope that someone else has experience with this error, and perhaps knows why OpenBUGS is not recognizing variable X even though it is clearly defined as far as I can tell.

我也遇到了错误expected the collection operator c error pos 8-这不是我以前遇到的错误，但是我同样感到困惑.

I've also gotten the error expected the collection operator c error pos 8 - this is not an error I've been getting previously, but I am similarly stumped.

模型和数据模拟功能均来自Kery的《面向生态学家的WinBUGS简介》(2010年).我会注意到这里的数据集代替了我自己的数据，这是相似的.

Both the model and the data-simulation function come from Kery's Introduction to WinBUGS for Ecologists (2010). I will note that the data set here is in lieu of my own data, which is similar.

我包括构建数据集和模型的功能.很抱歉的长度.

I am including the function to build the dataset as well as the model. Apologies for the length.

# Simulate data: 200 sites, 3 sampling rounds, 3 factors of the level 'trt', 
# and continuous covariate 'X'

data.fn <- function(nsite = 180, nrep = 3, xmin = -1, xmax = 1, alpha.vec = c(0.01,0.2,0.4,1.1,0.01,0.2), beta0 = 1, beta1 = -1, ntrt = 3){
  y <- array(dim = c(nsite, nrep))  # Array for counts
  X <- sort(runif(n = nsite, min = xmin, max = xmax))   # covariate values, sorted
  # Relationship expected abundance - covariate
  x2 <- rep(1:ntrt, rep(60, ntrt)) # Indicator for population
  trt <- factor(x2, labels = c("CT", "CM", "CC"))
  Xmat <- model.matrix(~ trt*X)
  lin.pred <- Xmat[,] %*% alpha.vec # Value of lin.predictor
  lam <- exp(lin.pred)
  # Add Poisson noise: draw N from Poisson(lambda)
  N <- rpois(n = nsite, lambda = lam)
  table(N)                # Distribution of abundances across sites
  sum(N > 0) / nsite          # Empirical occupancy
  totalN <- sum(N)  ;  totalN
  # Observation process
  # Relationship detection prob - covariate
  p <- plogis(beta0 + beta1 * X)
  # Make a 'census' (i.e., go out and count things)
  for (i in 1:nrep){
    y[,i] <- rbinom(n = nsite, size = N, prob = p)
  }
  # Return stuff
  return(list(nsite = nsite, nrep = nrep, ntrt = ntrt, X = X, alpha.vec = alpha.vec, beta0 = beta0, beta1 = beta1, lam = lam, N = N, totalN = totalN, p = p, y = y, trt = trt))
}

data <- data.fn()

这是模型:

sink("nmix1.txt")
cat("
    model {

    # Priors
    for (i in 1:3){     # 3 treatment levels (factor)   
    alpha0[i] ~ dnorm(0, 0.01)       
    alpha1[i] ~ dnorm(0, 0.01)       
    }
    beta0 ~ dnorm(0, 0.01)       
    beta1 ~ dnorm(0, 0.01)

    # Likelihood
    for (i in 1:180) {      # 180 sites
    C[i] ~ dpois(lambda[i])
    log(lambda[i]) <- log.lambda[i]
    log.lambda[i] <- alpha0[trt[i]] + alpha1[trt[i]]*X[i]

    for (j in 1:3){     # each site sampled 3 times
    y[i,j] ~ dbin(p[i,j], C[i])
    lp[i,j] <- beta0 + beta1*X[i]
    p[i,j] <- exp(lp[i,j])/(1+exp(lp[i,j]))
    }
    }

    # Derived quantities

    }
    ",fill=TRUE)
sink()

# Bundle data
trt <- data$trt
y <- data$y
X <- data$X
ntrt <- 3

# Standardise covariates
s.X <- (X - mean(X))/sd(X)

win.data <- list(C = y, trt = as.numeric(trt), X = s.X)

# Inits function
inits <- function(){ list(alpha0 = rnorm(ntrt, 0, 2), 
                          alpha1 = rnorm(ntrt, 0, 2),
                beta0 = rnorm(1,0,2), beta1 = rnorm(1,0,2))}

# Parameters to estimate
parameters <- c("alpha0", "alpha1", "beta0", "beta1")

# MCMC settings
ni <- 1200
nb <- 200
nt <- 2
nc <- 3

# Start Markov chains
out <- bugs(data = win.data, inits, parameters, "nmix1.txt", n.thin=nt, 
            n.chains=nc, n.burnin=nb, n.iter=ni, debug = TRUE)

推荐答案

注意:在我注意到代码的另一个问题之后，该答案已进行了重大修订.

Note: This answer has gone through a major revision, after I noticed another problem with the code.

如果我正确理解了您的模型，则您是从模拟数据中混合了 y 和 N ，以及作为 C 传递的内容虫子.您正在将 y 变量(矩阵)传递给Bugs模型中的C变量，但这是作为向量访问的.从我可以看到， C 代表您的二项式抽取(实际丰度)中的试验"数量，即数据集中的 N .变量 y (矩阵)在模拟数据和Bugs模型中都被称为同一事物.

If I understand your model correctly, you are mixing up the y and N from the simulated data, and what is passed as C to Bugs. You are passing the y variable (a matrix) to the C variable in the Bugs model, but this is accessed as a vector. From what I can see C is representing the number of "trials" in your binomial draw (actual abundances), i.e. N in your data set. The variable y (a matrix) is called the same thing in both the simulated data and in the Bugs model.

据我所知，这是您模型的重新构造，并且运行正常:

This is a reformulation of your model, as I understand it, and this runs ok:

sink("nmix1.txt")
cat("
    model {

    # Priors
    for (i in 1:3){     # 3 treatment levels (factor)   
    alpha0[i] ~ dnorm(0, 0.01)       
    alpha1[i] ~ dnorm(0, 0.01)       
    }
    beta0 ~ dnorm(0, 0.01)       
    beta1 ~ dnorm(0, 0.01)

    # Likelihood
    for (i in 1:180) {      # 180 sites
    C[i] ~ dpois(lambda[i])
    log(lambda[i]) <- log.lambda[i]
    log.lambda[i] <- alpha0[trt[i]] + alpha1[trt[i]]*X[i]

    for (j in 1:3){     # each site sampled 3 times
        y[i,j] ~ dbin(p[i,j], C[i])
        lp[i,j] <- beta0 + beta1*X[i]
        p[i,j] <- exp(lp[i,j])/(1+exp(lp[i,j]))
    }
    }

    # Derived quantities

    }
    ",fill=TRUE)
sink()

# Bundle data
trt <- data$trt
y <- data$y
X <- data$X
N<- data$N
ntrt <- 3

# Standardise covariates
s.X <- (X - mean(X))/sd(X)

win.data <- list(y = y, trt = as.numeric(trt), X = s.X, C= N)

# Inits function
inits <- function(){ list(alpha0 = rnorm(ntrt, 0, 2), 
                          alpha1 = rnorm(ntrt, 0, 2),
                beta0 = rnorm(1,0,2), beta1 = rnorm(1,0,2))}

# Parameters to estimate
parameters <- c("alpha0", "alpha1", "beta0", "beta1")

# MCMC settings
ni <- 1200
nb <- 200
nt <- 2
nc <- 3

# Start Markov chains
out <- bugs(data = win.data, inits, parameters, "nmix1.txt", n.thin=nt, 
            n.chains=nc, n.burnin=nb, n.iter=ni, debug = TRUE)

总体而言，该模型的结果看起来还不错，但是beta0和beta1的自相关滞后时间很长. beta1的估计值似乎也有些偏离(〜= -0.4)，因此您可能需要重新检查Bugs模型规范，以使其与仿真模型匹配(即您正在拟合正确的统计模型).目前，我不确定是否可以，但是我现在没有时间进一步检查.

Overall, the results from this model looks ok, but there are long autocorrelation lags for beta0 and beta1. The estimate of beta1 also seems a bit off(~= -0.4), so you might want to recheck the Bugs model specification, so that it is matching the simulation model (i.e. that you are fitting the correct statistical model). At the moment, I'm not sure that it does, but I don't have the time to check further right now.

这篇关于OpenBUGS错误未定义变量的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！