model
{
	Tau.noninformative <- 1.0E-2
	P.gamma <- 1.0E-2
	for (i in 1:N.sample) {
		# observation
		Area[i] ~ dnorm(area[i], tau.err)
		Bno[i] ~ dpois(lambda[i])
		area[i] <- exp(
			betaB[1]
			+ betaS[Sp[i], 1]
			+ betaT[No[i], 1]
		)
		lambda[i] <- exp(
			betaB[2]
			+ betaS[Sp[i], 2]
			+ betaT[No[i], 2]
			+ exp(betaB[3] + betaS[Sp[i], 3]) * LogH[i]
		)

	}
	# Parameters and hyper parameters
	tau.err ~ dgamma(P.gamma, P.gamma)
	for (k in 1:N.betaB) {
		betaB[k] ~ dnorm(0, Tau.noninformative)
	}
	for (j in 1:N.sp) {
		betaS[j, 1:N.betaS] ~ dmnorm(V.zero[], inv.vc[,])
	}
	inv.vc[1:N.betaS,1:N.betaS] ~ dwish(R[,], Deg.freedom)
	vc[1:N.betaS,1:N.betaS] <- inverse(inv.vc[,])
	for (x in 1:N.betaS) {
		for (y in 1:N.betaS) {
			rho[x, y] <- vc[x, y] / sqrt(vc[x, x] * vc[y, y])
		}
	}
	for (k in 1:N.betaT) {
		for (j in 1:N.tree) {
			betaT[j, k] ~ dnorm(0, tau.betaT[k])
		}
		tau.betaT[k] ~ dgamma(P.gamma, P.gamma)
	}
}