# Oikos 			oik-00233
# Sundell, J., O’Hara, R. B., Helle, P., Hellstedt, P., Henttonen, H. and Pietiäinen, H. 2013. 
# Numerical response of small mustelids to vole abundance: delayed or not? – Oikos 000: 000-000.

# Appendix 1
# OpenBugs’ model code

model {
	for(s in 1:NSeries) {
		for(t in 1:NYears) {
			MeanDate[s,t] ~ dnorm(0,1)
			MeanSnow[s,t] ~ dnorm(0,1)
			MeanTemp[s,t] ~ dnorm(0,1)
		}
		for(sp in 1:NSp) {			mu[sp,s,1] ~ dnorm(0, 0.01)	}
		for(t in 1:NYears) {
			VoleSpring[s,t] ~ dnorm(mu[1,s,t], tau[1,s])		# trapped May
			VoleAutumn[s,t] ~ dnorm(mu[2,s,t], tau[2,s])
			Weasel[s,t] ~ dnorm(mu[3,s,t], tau[3,s])			# trapped Jan/Feb
			Stoat[s,t] ~ dnorm(mu[4,s,t], tau[4,s])			# trapped Jan/Feb
		}
		for(t in 2:NYears) {
			for(sp in 1:NSp) {
				mean[sp,s,t] <- beta.Date[sp]*MeanDate[s,t] + mu.00[sp] + 
						 beta.Snow[sp]*MeanSnow[s,t] + beta.Temp[sp]*MeanTemp[s,t] + W[sp,s]
			}
			mu[1,s,t] <- (mean[1,s,t] + rho[1,1]*VoleSpring[s,t-1] + rho[1,2]*VoleAutumn[s,t-1] + 
									rho[1,3]*Weasel[s,t] + rho[1,4]*Stoat[s,t] + 
									rho[1,5]*Weasel[s,t-1]) + rho[1,6]*Stoat[s,t-1]
			mu[2,s,t] <- (mean[2,s,t] + rho[2,1]*VoleSpring[s,t] + rho[2,2]*VoleAutumn[s,t-1] + 
									rho[2,3]*Weasel[s,t] + rho[2,4]*Stoat[s,t] + 
									rho[2,5]*Weasel[s,t-1]) + rho[2,6]*Stoat[s,t-1]
			mu[3,s,t] <- (mean[3,s,t] + rho[3,1]*VoleSpring[s,t-1] + rho[3,2]*VoleAutumn[s,t-1] + 
									rho[3,3]*Weasel[s,t-1]) + rho[3,4]*Stoat[s,t-1]
			mu[4,s,t] <- (mean[4,s,t] + rho[4,1]*VoleSpring[s,t-1] + rho[4,2]*VoleAutumn[s,t-1] + 
									rho[4,3]*Weasel[s,t-1]) + rho[4,4]*Stoat[s,t-1]
		}

		for(sp in 1:NSp) {
			sd[sp,s] ~ dunif(0,100);	tau[sp,s] <- pow(sd[sp,s],-2)
		}
		VarPopSeries[1,s] <- pow(sd(VoleSpring[s,]),2)
		VarPopSeries[2,s] <- pow(sd(VoleAutumn[s,]),2)
		VarPopSeries[3,s] <- pow(sd(Weasel[s,]),2)
		VarPopSeries[4,s] <- pow(sd(Stoat[s,]),2)
	}
	for(sp in 1:NSp) {
		mu.00[sp]  ~ dnorm(0, 0.001)
#		mean.00[sp] <- mu.00[sp] + mean(W[sp,])

		sd0[sp] ~ dunif(0,100);	tau0[sp] <- pow(sd0[sp],-2)
		for(s in 1:NSeries) {
			W[sp,s] ~ dnorm(mu.00[sp], tau0[sp])
			mean.pop[sp,s] <- W[sp,s] - mean(W[sp,])	
		}
		for(spp in 1:NSp+2) {
			rho[sp,spp] ~ dnorm(0, 0.001)
			rhoST[sp,spp] <- rho[sp,spp]*sqrt(VarSeries[spp])
		}
		
		VarSeries[sp] <- mean(VarPopSeries[sp,])
	}
	VarSeries[5] <- VarSeries[3];			VarSeries[6] <- VarSeries[4];	

	for(sp in 1:2) {
# No model selection
		beta.Date[sp+2] ~ dnorm(0, 0.01);	beta.Date[sp] <- 0
		beta.Snow[sp+2] ~ dnorm(0, 0.01);	beta.Snow[sp] <- 0
		beta.Temp[sp+2] ~ dnorm(0, 0.01);	beta.Temp[sp] <- 0
	}
}