[R] Openbugs- Array Index
Jeff Newmiller
jdnewmil at dcn.davis.ca.us
Sat Oct 27 07:36:02 CEST 2012
You seem to be lost. Try http://mathstat.helsinki.fi/openbugs/community/CommBBS.html
---------------------------------------------------------------------------
Jeff Newmiller The ..... ..... Go Live...
DCN:<jdnewmil at dcn.davis.ca.us> Basics: ##.#. ##.#. Live Go...
Live: OO#.. Dead: OO#.. Playing
Research Engineer (Solar/Batteries O.O#. #.O#. with
/Software/Embedded Controllers) .OO#. .OO#. rocks...1k
---------------------------------------------------------------------------
Sent from my phone. Please excuse my brevity.
SuzieK <suzanne.f.kenny at hotmail.co.uk> wrote:
>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 at 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.
More information about the R-help
mailing list