geese.R

library(runjags)
library(R2WinBUGS)
library(R2jags)
library(BRugs)
library(gtools)
library(gdata)
library(siar)

# Load in the geese data from SIAR
u = as.matrix(sourcesdemo[,c(2,4)] + correctionsdemo[,c(2,4)])
sigma2 = sourcesdemo[,c(3,5)]^2 + correctionsdemo[,c(3,5)]^2
X = geese1demo
N = dim(X)[1]
num.prey = dim(u)[1]
num.iso = dim(u)[2]

#These are the parameters that need to be set by the user
  mcmc.chainLength <- as.integer(10000)  # post-burn
  mcmc.burn <- as.integer(5000)
  mcmc.thin = 1;
  mcmc.chains = 3       # needs to be at least 2 for DIC
  alpha=rep(1,num.prey)
  
  jags.data = list("u", "sigma2", "N", "num.prey", "num.iso", "num.pop","X","alpha")  

  model.loc = "MixSIR.txt"
  jags.params = c("p")
  jags.1 = jags(jags.data, inits = NULL, parameters.to.save= jags.params, model.file=model.loc, n.chains = mcmc.chains, n.burnin = mcmc.burn, n.thin = mcmc.thin, n.iter = mcmc.chainLength, DIC = TRUE) 
  attach.jags(jags.1)
  p.mixSIR = p
  print(jags.1$BUGSoutput$DIC)
  
  model.loc = "SIAR.txt"
  jags.params = c("p","resid.sigma")  
  jags.2 = jags(jags.data, inits = NULL, parameters.to.save= jags.params, model.file=model.loc, n.chains = mcmc.chains, n.burnin = mcmc.burn, n.thin = mcmc.thin, n.iter = mcmc.chainLength, DIC = TRUE) 
  attach.jags(jags.2)  
  p.SIAR = p
   print(jags.2$BUGSoutput$DIC)
  
  model.loc = "individual.txt"
  jags.params = c("p","p.mean","ind.sig")  
  jags.3 = jags(jags.data, inits = NULL, parameters.to.save= jags.params, model.file=model.loc, n.chains = mcmc.chains, n.burnin = mcmc.burn, n.thin = mcmc.thin, n.iter = mcmc.chainLength, DIC = TRUE) 
  attach.jags(jags.3)  
  p.ind.mean = p.mean
  p.ind = p
   print(jags.3$BUGSoutput$DIC)
  
  model.loc = "individual-residual.txt"
  jags.params = c("p","p.mean", "ind.sig", "resid.sig")  
  jags.4 = jags(jags.data, inits = NULL, parameters.to.save= jags.params, model.file=model.loc, n.chains = mcmc.chains, n.burnin = mcmc.burn, n.thin = mcmc.thin, n.iter = mcmc.chainLength, DIC = TRUE) 
  attach.jags(jags.4)  
  p.indresid.mean = p.mean
  p.indresid = p  
  print(jags.4$BUGSoutput$DIC)
  
  # Make density plots of the 4 sources
  par(mfrow = c(2,2))
  cols = c("blue1","tomato1","green3","purple")
  for(i in 1:4) {
    d = density(p.mixSIR[,i], from=0, to = 1)	
  	d$y = d$y/sum(d$y)
  	plot(d$x,d$y, xlab=sourcesdemo[i,1], ylab="Density",main="", lwd = 2, col = cols[1], type="l")
    d = density(p.SIAR[,i], from=0, to = 1)	
  	d$y = d$y/sum(d$y)  	
  	lines(d$x,d$y,lwd=2,col=cols[2])
    d = density(p.ind.mean[,i], from=0, to = 1)	
  	d$y = d$y/sum(d$y)  	
  	lines(d$x,d$y,lwd=2,col=cols[3])  	
    d = density(p.indresid.mean[,i], from=0, to = 1)	
  	d$y = d$y/sum(d$y)  	
  	lines(d$x,d$y,lwd=2,col=cols[4]) 
  	if(i==1) {
  		legend('topleft',legend=c("SIAR","MixSIR","Ind","Ind-resid"),lwd=2,lty=1,col=cols,inset=0.05)	
  	} 	
  }