model
{
	Tau.noninformative <- 1.0E-3
	Hyper.gamma <- 1.0E-2

	# stem number
	for (quadrat in 1:N.quadrat) {
		for (tspc in 1:N.tree) {
			Tree.stem[quadrat, tspc] ~ dpois(lambda[quadrat, tspc])
			lambda[quadrat, tspc] <- exp(log.lambda[quadrat, tspc])
		}
	}

	# lambda and mean.logstem
	for (quadrat in 1:N.quadrat) {
		for (tspc in 1:N.tree) { # random effects of quadrat x tspc
			log.lambda[quadrat, tspc] ~ dnorm(mean.logstem[quadrat, tspc], tau[1])
			mean.logstem[quadrat, tspc] <- (
				bb[1] + bs[1, tspc]
				+ (bb[2] + bs[2, tspc]) * shading[quadrat]
				+ re.slope.spc[Slope[quadrat], tspc]
				+ re.slope[Slope[quadrat]]
				+ re.quadrat[quadrat]
			)
		}
	}

	# bb and bs
	for (k in 1:N.b) {
		bb[k] ~ dnorm(0.0, Tau.noninformative)
		for (tspc in 1:N.tree) {
			bs[k, tspc] ~ dnorm(0.0, tau.b[k])
		}
		tau.b[k] ~ dgamma(Hyper.gamma, Hyper.gamma)
	}

	# shading and w.herb
	for (quadrat in 1:N.quadrat) {
		shading[quadrat] <- 1 - exp(-inprod(w.herb[], Herb[quadrat,]))
	}
	w.herb[1] <- 1.0
	for (hspc in 2:N.herb) {
		w.herb[hspc] <- exp(log.w.herb[hspc])
	}
	for (hspc in 1:N.herb) {
		log.w.herb[hspc] ~ dnorm(0.0, tau[2])
	}

	# re.slope.spc and re.slope
	for (slope in 1:N.slope) {
		re.slope.spc[slope, 1:N.tree] ~ dmnorm(Zero[], tau.slope.spc[,])
		re.slope[slope] ~ dnorm(0.0, tau[3])
	}
	tau.slope.spc[1:N.tree, 1:N.tree] ~ dwish(R[,], Deg.freedom)
	var.cov[1:N.tree, 1:N.tree] <- inverse(tau.slope.spc[,])
	for (i in 1:N.tree) {
		for (j in 1:N.tree) {
			rho[i, j] <- var.cov[i, j] / sqrt(var.cov[i, i] * var.cov[j, j])
		}
	}

	# re.quadrat
	for (quadrat in 1:N.quadrat) {
		re.quadrat[quadrat] ~ dnorm(0.0, tau[4])
	}

	# tau = 1 / variance
	for (k in 1:N.tau) {
		tau[k] ~ dgamma(Hyper.gamma, Hyper.gamma)
	}
}