wtitCounts.Bayes1 = function(){ # input data from files data<- read.table("wtmatrix.csv",header=TRUE,sep=",",na.strings=c("NA")) ymat<-as.matrix(data[,c("c.1","c.2","c.3")]) elev<-as.vector(data[,"elev"]) forest<-as.vector(data[,"forest"]) length<-as.vector(data[,"length"]) forest<-as.vector(data[,"forest"]) elev = (elev-mean(elev,na.rm=T))/sd(elev, na.rm=T) forest = (forest-mean(forest,na.rm=T))/sd(forest,na.rm=T) date<-as.matrix(data[,c("day.1","day.2","day.3")]) mdate<-mean(date,na.rm=TRUE) sddate<-sqrt(var(date[1:length(date)],na.rm=TRUE)) date<-(date-mdate)/sddate date2<-date*date elev2<-elev*elev dur<-as.matrix(data[,c("dur.1","dur.2","dur.3")]) intensity<-dur/length mint<-mean(intensity,na.rm=TRUE) sdint<-sqrt(var(intensity[1:length(intensity)],na.rm=TRUE)) intensity<-(intensity-mint)/sdint length<- (length-mean(length))/sqrt(var(length)) # vectorize matrices to eliminate missing values nsites = dim(ymat)[1] nrepls = dim(ymat)[2] site = rep(1:nsites, nrepls) y = as.vector(ymat) ymax = apply(ymat,1,max, na.rm=T) intensity = as.vector(intensity) date = as.vector(date) date2 = as.vector(date2) notmiss = !is.na(y) site = site[notmiss] y = y[notmiss] intensity = intensity[notmiss] date = date[notmiss] date2 = date2[notmiss] # model specification in WinBUGS modelFilename = "wtitBayes1.txt" cat(" model { # priors for (k in 1:5) { beta[k] ~ dnorm(0,.01) } alpha1Prec <- 1/pow(1.6,2) alpha[1] ~ dnorm(0,alpha1Prec) for (k in 2:4) { alpha[k] ~ dnorm(0,.01) } # abundance model for (i in 1:nsites) { N[i] ~ dpois(lambda[i]) log(lambda[i]) <- beta[1] + beta[2]*elev[i] + beta[3]*elev2[i] + beta[4]*forest[i] + beta[5]*length[i] } # observation model for (j in 1:nobs) { eta[j] <- alpha[1] + alpha[2]*intensity[j] + alpha[3]*date[j] + alpha[4]*date2[j] p[j] <- 1 / (1 + exp(-eta[j])) y[j] ~ dbin(p[j], N[site[j]]) } } ", fill=TRUE, file=modelFilename) # arguments for bugs() data = list(nsites=nsites, nobs=length(y), y=y, site=site, intensity=intensity, date=date, date2=date2, elev=elev, elev2=elev2, forest=forest, length=length) params = list('alpha','beta') inits = function() { list(alpha = rnorm(4), beta=rnorm(5), N=ymax+1) } # call to bugs() library(R2WinBUGS) fit = bugs(data, inits, params, model.file=modelFilename, n.chains=1, n.iter=60000, n.burnin=10000, n.thin=5, bugs.seed=sample(1:9999,size=1), debug=FALSE, DIC=FALSE) fit }