model
{
	Tau.noninformative <- 1.0E-2
	P.gamma <- 1.0E-2
	for (i in 1:N.leaf) {
		# observation
		Area[i] ~ dnorm(mean.area[i], tau.err[1])
		mean.area[i] <- exp(
			betaB[1] + betaS[SpL[i], 1]
			+ retree[Tree[i], 1]
		)
	}
	for (i in 1:N.tree) {
		Bno[i] ~ dpois(lambda[i])
		lambda[i] <- exp(
			betaB[2] + betaS[Sp[i], 2]
			+ retree[i, 2]
			+ exp(phiB[1] + phiS[Sp[i], 1]) * log.h[i]
		)
		LogH[i] ~ dnorm(log.h[i], tau.err[2])
		log.h[i] <- (
			betaB[3] + betaS[Sp[i], 3]
			+ retree[i, 3]
			+ exp(phiB[2] + phiS[Sp[i], 2]) * log.d10[i]
		)
		LeafT[i] ~ dpois(lambda.leaf[i])
		lambda.leaf[i] <- exp(
			betaB[4] + betaS[Sp[i], 4]
			+ retree[i, 4]
			+ exp(phiB[3] + phiS[Sp[i], 3]) * log.d10[i]
		)
		log.d10[i] ~ dnorm(LogD10[i], tau.err[3])
	}
	# Parameters and hyper parameters
	for (k in 1:N.beta) {
		betaB[k] ~ dnorm(0, Tau.noninformative)
	}
	for (j in 1:N.sp) {
		betaS[j, 1:N.beta] ~ dmnorm(V.zero[], inv.vc[,])
	}
	inv.vc[1:N.beta,1:N.beta] ~ dwish(R[,], Deg.freedom)
	vc[1:N.beta,1:N.beta] <- inverse(inv.vc[,])
	for (x in 1:N.beta) {
		for (y in 1:N.beta) {
			rho[x, y] <- vc[x, y] / sqrt(vc[x, x] * vc[y, y])
		}
	}
	for (k in 1:N.phi) {
		phiB[k] ~ dnorm(0, Tau.noninformative)
		for (j in 1:N.sp) {
			phiS[j, k] ~ dnorm(0, tau.phi[k])
		}
		tau.phi[k] ~ dgamma(P.gamma, P.gamma)
	}
	for (k in 1:N.retree) {
		for (j in 1:N.tree) {
			retree[j, k] ~ dnorm(0, tau.retree[k])
		}
		tau.retree[k] ~ dgamma(P.gamma, P.gamma)
	}
	for (k in 1:N.tau.err) {
		tau.err[k] ~ dgamma(P.gamma, P.gamma)
	}
}