sgo - 5 months ago 65

R Question

I'm working on a binomial mixture model using OpenBUGS and R package

`R2OpenBUGS`

`variable X is not defined in model or in data set`

I've also gotten the error

`expected the collection operator c error pos 8`

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()

And here is the model:

`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)

Answer

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

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)
```

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.