Hi,
I'm working on the codes below however every time I run them when they get
to OpenBUGS I keep getting the error message: array index is greater than
array upper bound for hab.
Any help would be greatly appreciated,
Suzie
Codes:
ungulate <- read.csv(file.choose ()) #ungulate
ungulate <-
as.matrix(ungulate);colnames(ungulate)<-NULL;rownames(ungulate)<-NULL
habitat <- read.csv(file.choose ()) #ungulate habitat
habitat <- habitat[,3]
site.data <- read.csv(file.choose ()) #ungulate site
site <- site.data$SiteId
visit <- site.data$visit
date <- site.data$date
date <- matrix(date,nrow=4,ncol=5,byrow=TRUE)
S <- dim(ungulate)[1] # set dimension forungulate i.e. number of species
m <- 6 # number of augmented species
G <- read.csv(file.choose ()) # ungulate group
G <- G[,2]
g <- rep(NA,length=m)
G <- c(G,g)
g <- length(table(G))# number of groups
# habitat id
Hab <- cbind(seq(1,4),habitat)
hab <- matrix(NA,nrow=3,ncol=2) #make a matrix with 4 rows and two colums
for the four sites with 2 different habitat types has NA values??
for(i in 1:2){
hab[i,] <- Hab[habitat==i,1]
}
# arrange bird data
S <- dim(ungulate)[1]# number of collected species in our case
Y <- array(NA,dim=c(S,4,5))
for(i in 1:S){
Y[i,,] <- matrix(ungulate[i,],nrow=4,ncol=5,byrow=TRUE)
}
AY <- array(NA,dim=c((S+m),4,5))# Y of augmented
for(i in 1:4){
y <- matrix(0,nrow=m,ncol=5)
AY[,i,] <- rbind(Y[,i,],y)
}
## bugs code
library(R2OpenBUGS)
sink("ungulate.txt")
cat("
model{
# hyperparameters
# habitat effects for each functional group
for(i in 1:g){ # number of functional group
for(j in 1:2){ # number of habitat
type
mu.h[i,j] ~ dnorm(0,0.0001) I(-2,2) # filling mu.h with values based
on a random distribution
sigma.h[i,j] ~ dunif(0,6)
tau.h[i,j] <- 1/(sigma.h[i,j]*sigma.h[i,j])
}
}
# detectability
mu.r ~ dnorm(0,0.0001) I(-3,3)
sigma.r ~ dunif(0,6)
tau.r <- 1/(sigma.r*sigma.r)
psi ~ dunif(0,1)# inclusion rate that generates wi
# proportion of number of species among groups
for(i in 1:g){
prop[i] ~ dgamma(1,1)
prob[i] <- prop[i]/sum(prop[])
}
for(i in 1:(S+m)){
r[i] ~ dnorm(mu.r,tau.r) I(-2,2)# generating parameters related to
detectability
p[i] <- 1/(1+exp(-(r[i])))# individual-level detection probability
w[i] ~ dbern(psi)# indicator variable whether each species is exposed to
sampling or not
G[i] ~ dcat(prob[1:g])# group identity
for(h in 1:2){# habitat effects
habitat.eff[i,h] ~ dnorm(mu.h[G[i],h],tau.h[G[i],h]) I(-2,2)
}
for(j in 1:4){# fitting process
# ecological process model
lambda[i,j] <- exp(habitat.eff[i,habitat[j]])
Z[i,j] ~ dpois(lambda[i,j])# latent abundance of each species at each site
at each visit
A[i,j] <- Z[i,j]*w[i]# latent abundance only for species exposed to sampling
for(v in 1:5){
# detection process model
AY[i,j,v] ~ dbin(p[i],A[i,j])
}
}
# group identity (indicator variable)
G1[i] <- equals(G[i],1)
G2[i] <- equals(G[i],2)
}
for(i in 1:(S+m)){
for(j in 1:4){
A1[i,j] <- A[i,j]*G1[i]
A2[i,j] <- A[i,j]*G2[i]
O[i,j] <- step(A[i,j]-1)# latent occupancy of each species for each site
O1[i,j] <- O[i,j]*G1[i]
O2[i,j] <- O[i,j]*G2[i]
}
}
for(j in 1:4){
AB0[j] <- sum(A[,j])
AB1[j] <- sum(A1[,j])
AB2[j] <- sum(A2[,j])
SpR0[j] <- sum(O[,j])
SpR1[j] <- sum(O1[,j])
SpR2[j] <- sum(O2[,j])
}
for(i in 1:2){
for(j in 1:4){
HabAB0[i,j] <- AB0[hab[i,j]]
HabAB1[i,j] <- AB1[hab[i,j]]
HabAB2[i,j] <- AB2[hab[i,j]]
HabSpR0[i,j] <- SpR0[hab[i,j]]
HabSpR1[i,j] <- SpR1[hab[i,j]]
HabSpR2[i,j] <- SpR2[hab[i,j]]
}
HAB0[i] <- mean(HabAB0[i,])
HAB1[i] <- mean(HabAB1[i,])
HAB2[i] <- mean(HabAB2[i,])
HSpR0[i] <- mean(HabSpR0[i,])
HSpR1[i] <- mean(HabSpR1[i,])
HSpR2[i] <- mean(HabSpR2[i,])
}
R <- sum(w[1:(S+m)])# estimating unknown number of species that occupy any
sites
}
",fill=TRUE)
sink()
data <- list("m","S",
"habitat",
"G","g",
#"date",
"hab",
"AY")
inits <- function()list(
mu.h=matrix(rnorm(g*2),nrow=2),sigma.h=matrix(runif(g*2),nrow=2),
mu.r=rnorm(1),sigma.r=runif(1),
r=rnorm(S+m),
prop=runif(n=g,min=0,max=5),
Z=array(rpois(n=(S+m)*4,lambda=2),dim=c((S+m),4)),
w=rep(1,(S+m)),psi=runif(1))
parameters <- c(
"mu.h","sigma.h",
"habitat.eff",
"mu.r","sigma.r",
"r",
"R","psi",
"A",
"O",
"HAB0","HAB1","HAB2",
"HSpR0","HSpR1","HSpR2"
)
out <- bugs(data,inits,parameters,"ungulate.txt",
OpenBUGS.pgm="C:/Program Files (x86)/OpenBUGS/OpenBUGS322/OpenBUGS.exe",
n.thin=10,n.burnin=10,n.chains=1,n.iter=1000,debug=TRUE)
--
View this message in context:
http://r.789695.n4.nabble.com/Openbugs-Array-Index-tp4647587.html
Sent from the R help mailing list archive at Nabble.com.
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
