'Attempt to redefine node p[1,1] When adding covariates to beta-binomial mixture model
I keep getting the above error when I try adding detection covariates on a beta-binomial N-mixture model in rags. According to Royle(2004). A binomial N mixture model can be used to model abundance data arising from repeat count surveys. The number of individuals on a site can be modeled by a Poisson [For simplicity I will stick to the Poisson model only] such that;
Ni ~ Poisson(λi)
yit ~ bin(pit,Ni)
Ni - is the number of animals available in site i
yit - is the count of observed animals at site i visit t
λi - is the average number of animals at site i
pit - is the detection probability.
Covariate effects can be modeled as:
Abudance:
log(λi)= B0+ B1xi1 +...+Brxir where 1...r are covariates
detection:
logit(pit)= B0+ B1xi1 +...+Brxir where 1...r are covariates
Probabilty of detection pit is assumed to be contant for all present animals.
The ** beta binomial model ease this assumption* by letting pit follow a stochastic distribution instead, such that
pi't ~ Beta(pit(1-δ2)/δ2,(1-pit)(1-δ2)/δ2
for 0<δ<1
I tried implementing the model with a simulated data:
20 sites, 5 visits, site covariate = Location, and 2 observed covariates.
simulated data
library(modelr)
Location<-c("A","B","C","D")
Location<-data.frame(Location=rep(Location,5))
location=Location%>%model_matrix(~Location)%>%select(-1)
set.seed(100)
y<-matrix(rpois(100,0.5),ncol=5)
# Cov1
set.seed(100)
cov1<-matrix(rnorm(100,100,5),ncol=5)
# cov 2
set.seed(100)
cov2<-matrix(rnorm(100,50,2),ncol=5)
data<-list(y=y,
nSites=20,
nOcc=5,
location=location,
cov1=cov1,
cov2=cov2)
if i try estimating this model in rjags without the covariate effects for detection it works.
nx<-"
model{
for(i in 1:nSites) {
# Biological model
N[i] ~ dpois(lambda[i])
log(lambda[i])<-alphao+inprod(alpha, location[i,])
# Observation model
for(j in 1:nOcc) {
C[i,j] ~ dbin(p[i,j], N[i])
p[i,j]~ dbeta(1,1)
}
}
# Priors
alphao ~ dnorm(5,1)
for(i in 1:nA){
alpha[i] ~ dnorm(2,.5)
}
pit_h ~ dunif(0,1)
sigma ~ dunif(0,1)
fac <-(1-sigma^2)/sigma^2
A<- pit_h*fac
B<-(1-pit_h)*fac
}"
writeLines(nx,con="mod.txt")
watch=c("alpha0","alpha","lambda","pit_h")
mod<-jagsUI::jags(data,
parameters.to.save=watch,
model.file="C:/Users/user/Documents/mod.txt",n.iter=3,
n.chains=2)
yielding
JAGS output for model 'C:/Users/user/Documents/mod.txt', generated by jagsUI.
Estimates based on 2 chains of 3 iterations,
adaptation = 100 iterations (sufficient),
burn-in = 0 iterations and thin rate = 1,
yielding 6 total samples from the joint posterior.
MCMC ran for 0.166 minutes at time 2022-04-02 21:11:21.
mean sd 2.5% 50% 97.5% overlap0 f Rhat n.eff
alpha[1] 2.355 1.476 0.785 1.962 4.319 FALSE 1 0.990 6
alpha[2] 2.800 1.441 0.766 2.974 4.435 FALSE 1 0.854 6
alpha[3] 2.919 1.568 0.810 2.941 4.757 FALSE 1 0.838 6
lambda[1] 195.074 131.799 41.549 208.461 333.312 FALSE 1 0.844 6
lambda[2] 3060.224 3984.560 326.627 1431.883 9905.889 FALSE 1 1.494 5
lambda[3] 9518.864 11955.073 111.455 5098.236 28641.383 FALSE 1 0.958 6
lambda[4] 8610.148 8704.268 116.832 6850.301 20351.308 FALSE 1 0.896 6
lambda[5] 195.074 131.799 41.549 208.461 333.312 FALSE 1 0.844 6
lambda[6] 3060.224 3984.560 326.627 1431.883 9905.889 FALSE 1 1.494 5
lambda[7] 9518.864 11955.073 111.455 5098.236 28641.383 FALSE 1 0.958 6
lambda[8] 8610.148 8704.268 116.832 6850.301 20351.308 FALSE 1 0.896 6
lambda[9] 195.074 131.799 41.549 208.461 333.312 FALSE 1 0.844 6
lambda[10] 3060.224 3984.560 326.627 1431.883 9905.889 FALSE 1 1.494 5
lambda[11] 9518.864 11955.073 111.455 5098.236 28641.383 FALSE 1 0.958 6
lambda[12] 8610.148 8704.268 116.832 6850.301 20351.308 FALSE 1 0.896 6
lambda[13] 195.074 131.799 41.549 208.461 333.312 FALSE 1 0.844 6
lambda[14] 3060.224 3984.560 326.627 1431.883 9905.889 FALSE 1 1.494 5
lambda[15] 9518.864 11955.073 111.455 5098.236 28641.383 FALSE 1 0.958 6
lambda[16] 8610.148 8704.268 116.832 6850.301 20351.308 FALSE 1 0.896 6
lambda[17] 195.074 131.799 41.549 208.461 333.312 FALSE 1 0.844 6
lambda[18] 3060.224 3984.560 326.627 1431.883 9905.889 FALSE 1 1.494 5
lambda[19] 9518.864 11955.073 111.455 5098.236 28641.383 FALSE 1 0.958 6
lambda[20] 8610.148 8704.268 116.832 6850.301 20351.308 FALSE 1 0.896 6
pit_h 0.434 0.340 0.024 0.491 0.881 FALSE 1 0.967 6
**WARNING** Rhat values indicate convergence failure.
Rhat is the potential scale reduction factor (at convergence, Rhat=1).
For each parameter, n.eff is a crude measure of effective sample size.
overlap0 checks if 0 falls in the parameter's 95% credible interval.
f is the proportion of the posterior with the same sign as the mean;
i.e., our confidence that the parameter is positive or negative.
>
if i ignore the stochastic distribution of pi't it works actually this is the constant pit i mentioned. everybody is doing it on every tutorial
The code
nx1<-"
model{
for(i in 1:nSites) {
# Biological model
N[i] ~ dpois(lambda[i])
log(lambda[i])<-alphao+inprod(alpha, location[i,])
# Observation model
for(j in 1:nOcc) {
C[i,j] ~ dbin(p[i,j], N[i])
#p[i,j]~ dbeta(1,1)
logit(p[i,j])<-beta0+inprod(beta1, location[i,])+inprod(beta2,c(cov1[i,j],cov2[i,j])) # beta1 are coefficients for loc, beta 2 are coef for other covs
}
}
# Priors
alphao ~ dnorm(5,1)
for(i in 1:3){
alpha[i] ~ dnorm(2,.5)
}
beta0 ~ dunif(0,0.5)
for(i in 1:3){
beta1[i] ~ dnorm(2,.5)
}
for(i in 1:2){
beta2[i] ~ dnorm(0.7,1)
}
pit_h ~ dunif(0,1)
sigma ~ dunif(0,1)
fac <-(1-sigma^2)/sigma^2
A<- pit_h*fac
B<-(1-pit_h)*fac
}"
writeLines(nx1,con="mod1.txt")
watch=c("alpha0","alpha","beta0","beta1","beta2","lambda","pit_h")
mod<-jagsUI::jags(data,
parameters.to.save=watch,
model.file="C:/Users/user/Documents/mod1.txt",n.iter=3,
n.chains=2)
yielding
JAGS output for model 'C:/Users/user/Documents/mod1.txt', generated by jagsUI.
Estimates based on 2 chains of 3 iterations,
adaptation = 100 iterations (sufficient),
burn-in = 0 iterations and thin rate = 1,
yielding 6 total samples from the joint posterior.
MCMC ran for 0.029 minutes at time 2022-04-02 21:19:02.
mean sd 2.5% 50% 97.5% overlap0 f Rhat n.eff
alpha[1] 1.833 1.142 -0.050 2.039 2.809 TRUE 0.833 1.085 6
alpha[2] 2.328 1.648 0.461 2.241 4.634 FALSE 1.000 0.931 6
alpha[3] 2.308 0.591 1.882 2.075 3.306 FALSE 1.000 2.012 3
beta0 0.234 0.140 0.041 0.256 0.391 FALSE 1.000 0.984 6
beta1[1] 1.679 1.872 -0.019 0.828 4.194 TRUE 0.833 2.415 3
beta1[2] 1.983 0.572 1.226 2.054 2.724 FALSE 1.000 0.924 6
beta1[3] 1.184 2.221 -1.716 2.022 3.333 TRUE 0.667 1.533 4
beta2[1] 1.931 0.736 1.301 1.592 2.949 FALSE 1.000 1.909 3
beta2[2] 0.079 0.789 -0.988 0.272 0.871 TRUE 0.667 0.855 6
lambda[1] 283.702 327.675 96.299 166.285 853.559 FALSE 1.000 1.382 5
lambda[2] 2215.821 1754.254 190.467 2134.896 4482.595 FALSE 1.000 1.303 5
lambda[3] 20986.924 48155.609 319.250 1187.748 104778.837 FALSE 1.000 1.306 6
lambda[4] 2488.915 1964.987 924.619 1971.019 5782.935 FALSE 1.000 1.050 6
lambda[5] 283.702 327.675 96.299 166.285 853.559 FALSE 1.000 1.382 5
lambda[6] 2215.821 1754.254 190.467 2134.896 4482.595 FALSE 1.000 1.303 5
lambda[7] 20986.924 48155.609 319.250 1187.748 104778.837 FALSE 1.000 1.306 6
lambda[8] 2488.915 1964.987 924.619 1971.019 5782.935 FALSE 1.000 1.050 6
lambda[9] 283.702 327.675 96.299 166.285 853.559 FALSE 1.000 1.382 5
lambda[10] 2215.821 1754.254 190.467 2134.896 4482.595 FALSE 1.000 1.303 5
lambda[11] 20986.924 48155.609 319.250 1187.748 104778.837 FALSE 1.000 1.306 6
lambda[12] 2488.915 1964.987 924.619 1971.019 5782.935 FALSE 1.000 1.050 6
lambda[13] 283.702 327.675 96.299 166.285 853.559 FALSE 1.000 1.382 5
lambda[14] 2215.821 1754.254 190.467 2134.896 4482.595 FALSE 1.000 1.303 5
lambda[15] 20986.924 48155.609 319.250 1187.748 104778.837 FALSE 1.000 1.306 6
lambda[16] 2488.915 1964.987 924.619 1971.019 5782.935 FALSE 1.000 1.050 6
lambda[17] 283.702 327.675 96.299 166.285 853.559 FALSE 1.000 1.382 5
lambda[18] 2215.821 1754.254 190.467 2134.896 4482.595 FALSE 1.000 1.303 5
lambda[19] 20986.924 48155.609 319.250 1187.748 104778.837 FALSE 1.000 1.306 6
lambda[20] 2488.915 1964.987 924.619 1971.019 5782.935 FALSE 1.000 1.050 6
pit_h 0.251 0.183 0.086 0.184 0.504 FALSE 1.000 0.843 6
**WARNING** Rhat values indicate convergence failure.
Rhat is the potential scale reduction factor (at convergence, Rhat=1).
For each parameter, n.eff is a crude measure of effective sample size.
overlap0 checks if 0 falls in the parameter's 95% credible interval.
f is the proportion of the posterior with the same sign as the mean;
i.e., our confidence that the parameter is positive or negative.
But if i include both the covariate effect for detection and the stochastic distribution for detection probability. things go south. see the code below
nx2<-"
model{
for(i in 1:nSites) {
# Biological model
N[i] ~ dpois(lambda[i])
log(lambda[i])<-alphao+inprod(alpha, location[i,])
# Observation model
for(j in 1:nOcc) {
C[i,j] ~ dbin(p[i,j], N[i])
p[i,j]~ dbeta(1,1)
logit(p[i,j])<-beta0+inprod(beta1, location[i,])+inprod(beta2,c(cov1[i,j],cov2[i,j])) # beta1 are coefficients for loc, beta 2 are coef for other covs
}
}
# Priors
alphao ~ dnorm(5,1)
for(i in 1:3){
alpha[i] ~ dnorm(2,.5)
}
beta0 ~ dunif(0,0.5)
for(i in 1:3){
beta1[i] ~ dnorm(2,.5)
}
for(i in 1:2){
beta2[i] ~ dnorm(0.7,1)
}
pit_h ~ dunif(0,1)
sigma ~ dunif(0,1)
fac <-(1-sigma^2)/sigma^2
A<- pit_h*fac
B<-(1-pit_h)*fac
}"
writeLines(nx2,con="mod2.txt")
watch=c("alpha0","alpha","beta0","beta1","beta2","lambda","pit_h")
mod<-jagsUI::jags(data,
parameters.to.save=watch,
model.file="C:/Users/user/Documents/mod2.txt",n.iter=3,
n.chains=2)
this is the error.
Error in jags.model(file = model.file, data = data, inits = inits, n.chains = n.chains, :
RUNTIME ERROR:
Compilation error on line 13.
Attempt to redefine node p[1,1]
I understand it is telling me that; p[i,j]~ dbeta(1,1) is being overwritten by logit(p[i,j])<-beta0+inprod(beta1, location[i,])+inprod(beta2,c(cov1[i,j],cov2[i,j])) but how exactly should this model be implemented. here the model is implemented without detection covariates. Its not what I am looking for.
Sources
This article follows the attribution requirements of Stack Overflow and is licensed under CC BY-SA 3.0.
Source: Stack Overflow
| Solution | Source |
|---|