4

OpenBUGS と R パッケージを使用して二項混合モデルに取り組んでいますR2OpenBUGS。より単純なモデルを正常に構築できましたが、不完全な検出のために別のレベルを追加すると、一貫してエラーが発生しますvariable X is not defined in model or in data set。データの構造を変更したり、データを OpenBUGS に直接入力したりするなど、さまざまなことを試しました。他の誰かがこのエラーを経験したことがあり、OpenBUGS が変数 X を認識していない理由を知っていることを願ってこれを投稿しています。

エラーも発生しました。expected the collection operator c error pos 8これは以前に発生したエラーではありませんが、同様に困惑しています。

モデルとデータ シミュレーション機能の両方は、Kery の「Introduction to WinBUGS for Ecoologists (2010)」に由来します。ここでのデータ セットは、同様の私自身のデータの代わりであることに注意してください。

データセットとモデルを構築する関数を含めています。長々と申し訳ありません。

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

2 に答える 2

0

因子を OpenBUGS に渡そうとすると、同じメッセージが表示されました。そのようです、

Ndata <- list(yrs=N$yrs, site=N$site), ... )

変数 "site" が "bugs" 関数によって渡されませんでした。OpenBUGS に渡されたリストになかっただけです

サイトを数値として渡すことで問題を解決しました。

Ndata <- list(yrs=N$yrs, site=as.numeric(N$site)), ... )
于 2014-12-19T14:31:48.663 に答える