library("BRugs")
library("R2WinBUGS")
###################################
# STEP 1: SET PARAMETERS
###################################
numClones<-50; # must be at least 2
chainLength <- 2000; # min 2000, any statistic = function of sample size
burn <- 4000; # min 4000
###################################
# STEP 2: READ IN THE DATA
###################################
d <- read.csv("dat.csv",header=TRUE)
# replace 0s with NAs
for(j in 1:dim(d)[2]) {
ind <- which(d[,j] < 0);
d[ind,j] <- NA;
d[-ind,j] <- log(d[-ind,j])
}
N <- dim(d)[2]; # number of sites
numY <- dim(d)[1]; # number of years
DATA <- matrix(0,numY,N);
for(i in 1:N) DATA[,i] <- d[,i];
#####################################################################
# USER NEEDS TO CHANGE THE FLAGS (F2, F3) AND THE GROUP VECTORS BELOW
# To summarize, there are N sites being monitored, and there may be M subpopulations (1 <= M <= N).
# The vector 'groups' assigns sites to subpopulations, and the remaining vectors allow for parameters
# to be shared among subpopulations. For example, if 2 subpopulations have the same process error vars,
# but different growth rates and observation errors, U.groups = c(1,2), Q.groups = c(1,1), R.groups = c(1,2)
#####################################################################
f1 <- 4; # FLAG FOR B (1, 2, 3, 4), interaction matrix
#case 1 % no interactions (identity matrix)
#case 2 % diagonal matrix, equal
#case 3 % diagonal matrix, unequal
#case 4 % unconstrained interaction matrix
f2 <- 1; # FLAG FOR Q (1, 2, 3), process error
f3 <- 2; # FLAG FOR R (1, 2, 3), observation error
#case 1 %unconstrained
#case 2 %diagonal covariance matrix, variances can be equal or not
#case 3 %equal variances and covariances; the same but not independent
# THESE NEXT 4 VECTORS CAN BE MODIFIED FOR DIFFERENT STRUCTURES
groups <- c(1,2,3,4); # this is an Nx1 vector describing which site belongs to which group
M <- max(groups); # M is the clusters / groups / signals, DO NOT CHANGE THIS LINE
if(M==1) f2 = 2; # this line prevents an error from being thrown
U.groups <- rep(1,M); # this is an Mx1 vector describing which group-specific growth rates are unique
#Q.groups <- seq(1,M);
Q.groups <- rep(1,M);# this is an Mx1 vector describing which group-specific process errors are unique
#R.groups <- seq(1,N);
R.groups <- rep(1,N); # this is an Mx1 vector describing which group-specific observation errors are unique
############################################################
# USER DOES NOT CHANGE THESE VALUES
############################################################
U.G <- max(U.groups); # this is the number of unique Us being estimated
Q.G <- max(Q.groups); # this is the number of unique Qs being estimated
R.G <- max(R.groups); # this is the number of unique Rs being estimated
source("writeModel.R");
# note: B can't be in the parameter list if it is fixed
if(f1 == 1) parameters <- c("U","Q","R","A","states1")
if(f1 > 1) parameters <- c("B","U","Q","R","A","states1")
mod <- bugs(data,inits=NULL,parameters,DIC=TRUE,"model.txt",n.chains=1,n.burnin=burn,
n.thin=1,n.iter=(burn+chainLength),program="openbugs",debug=TRUE)
attach.bugs(mod, overwrite=TRUE)
dev <- deviance/numClones;
nPar <- mean(dev)-min(dev)
#bestDev <- which(dev==min(dev))
# calculate the deviance at the mean posterior
pD <- mod$pD;
Dbar <- mean(deviance);
Dhat <- (Dbar - pD)/numClones;
Dhat <- (mod$DIC - 2*pD)/numClones;
z<- rep(0,N);
z[1] <- min(dev);
z[2] <- mean(dev);
z[3] <- nPar;
z[4] <- pD;
z[5] <- Dbar;
z[6] <- Dhat;
z[7] <- mod$DIC;
# calculate the distribution of the states
mat <- matrix(0,dim(d)[1],dim(d)[2]);
for(i in 1:dim(d)[1]) {
for(j in 1:dim(d)[2]) {
mat[i,j]<-(states1[1,i,j]);
}
}
if(f2 == 1) {
# unconstrained
Qmat <- matrix(0,dim(d)[2],dim(d)[2]);
if(length(Q) == chainLength) {
den <- density(Q);
Qmat[1,1] <- den$x[which(den$y==max(den$y))];
}
if(length(Q) > chainLength) {
for(i in 1:U.G) {
for(j in 1:U.G) {
den <- density(Q[,i,j]);
Qmat[i,j] <- den$x[which(den$y==max(den$y))];
}
}
}
}
if(f2 == 2) {
# constrained
Qmat <- matrix(0,dim(d)[2],dim(d)[2]);
if(length(Q) == chainLength) {
den <- density(Q);
Qmat[1,1] <- den$x[which(den$y==max(den$y))];
}
if(length(Q) > chainLength) {
for(i in 1:U.G) {
den <- density(Q[,i]);
Qmat[i,i] <- den$x[which(den$y==max(den$y))];
}
}
}
if(f3 == 1) {
# unconstrained
Rmat <- matrix(0,dim(d)[2],dim(d)[2]);
for(i in 1:dim(d)[2]) {
for(j in 1:dim(d)[2]) {
den <- density(R[,i,j]);
Rmat[i,j] <- den$x[which(den$y==max(den$y))];
}
}
}
if(f3 == 2) {
# constrained
Rmat <- matrix(0,dim(d)[2],dim(d)[2]);
for(i in 1:dim(d)[2]) {
den <- density(R[,i]);
Rmat[i,i] <- den$x[which(den$y==max(den$y))];
}
}