model
{
	for (i in 1:N.site) {
		# Observation
		Fish[i] ~ dbern(p.fish[i])
		BiomassZP[i] ~ dnorm(bzp[i], tau.err[1]) 
		Chl.a[i] ~ dnorm(chl.a[i], tau.err[2])
		log.din[i] ~ dnorm(log.DIN[i], tau.err[3])
		log.dtp[i] ~ dnorm(log.DTP[i], tau.err[4])
		W.temp[i] ~ dnorm(w.temp[i], tau.err[5]) 
		w.temp[i] ~ dnorm(0.0, tau.err[6])
		# log.density -> mean[i]
		p.fish[i] <- 1 - exp(-exp(log.density[i, 1]))
		log(bzp[i]) <- log.density[i, 2]
		log(chl.a[i]) <- log.density[i, 3]
		log.density[i, 1:N.level] ~ dmnorm(mu[i,], inv.vc[,])
		# Means
		mu[i,1] <- ( # log(fish density)
			beta.f[1]
			+ beta.f[2] * w.temp[i]
			+ beta.f[3] * log.Area[i]
		)
		mu[i,2] <- ( # log(biomass of zoo plankton)
			beta.z[1]
			+ beta.z[2] * w.temp[i]
		)
		mu[i,3] <- ( # log(density of Chl-a)
			beta.c[1]
			+ beta.c[2] * w.temp[i]
			+ beta.c[3] * log.din[i]
			+ beta.c[4] * log.dtp[i]
		)
	}
	# parameters
	for (k in 1:N.beta.f) {
		beta.f[k] ~ dnorm(0, Tau.noninformative)
	}
	for (k in 1:N.beta.z) {
		beta.z[k] ~ dnorm(0, Tau.noninformative)
	}
	for (k in 1:N.beta.c) {
		beta.c[k] ~ dnorm(0, Tau.noninformative)
	}
	# hyper priors
	for (k in 1:N.err) {
		tau.err[k] ~ dgamma(Err.gamma, Err.gamma)
	}
	inv.vc[1:N.level,1:N.level] ~ dwish(R[,], Deg.freedom)
	vc[1:N.level,1:N.level] <- inverse(inv.vc[,])
	for (l1 in 1:N.level) {
		for (l2 in 1:N.level) {
			rho[l1, l2] <- vc[l1,l2] / sqrt(vc[l1, l1] * vc[l2, l2])
		}
	}
}