R / JAGS code to implement the model
	The following code was run in R v3.0.2 (R core development team, 2013) and JAGS 3.1.0 using R2Jags 0.03-11 (Su and Yajima, 2013). 

rm(list=ls())

library(R2jags)

### Bayesian hierarchical model
sink('Ndep-model-matrix-corr-nomat.txt')
cat("
model { 
	# following Gelman and Hill p378-9 for multiple varying coefficients with multiple group-level predictors
	# using matrix notation
		# y = anpp data
		# y.hat = expectation
		# 
		# B = matrix of parameters including slope and intercept (K x J)
		# K = vector of named parameters (including intercept)
		# J = vector of groups
		# X = individual (plot) data matrix including intercept indicator (n x K)
		# G = group level parameters Iincluding intercept
		# U = group level (site) data matrix including intercept indicator
		# xi = 'multiplicative factor' for each parameter to scale inverse-Wishart covariance matrix
		# B.raw = matrix of random, group-varying plot-level parameters
		# B.raw.hat = matrix of expected group level of each parameter, based on G.raw x U
		# G.raw = hyperpriors drawn from dnorm(0,0.0001)
		# Tau.B.raw = inverse scaled Wishart matrix of covariance
		# W = identity matrix; diag(K)
		# rho.B = matrix (K x K) of parameter correlations 
		# sigma.B = vector of scaled standard deviations of parameters
		# Sigma.B.raw = variance-covariance matrix of group parameters
		
		
	# Likelihood of observations
		for (i in 1:n){
			y[i] ~ dnorm(y.hat[i], tau.y)
			y.hat[i] <- inprod(B[site[i],],X[i,]) # B x X is parameters times predictors
			err.y[i] <- y[i] - y.hat[i]	
			y.new[i] ~ dnorm(y.hat[i], tau.y)
			err.new[i] <- y.new[i] - y.hat[i]
			D[i] <- pow(err.y[i],2)
			D.new[i] <- pow(err.new[i],2) 		
		}
		fit <- sum(D[])
		fit.new <- sum(D.new[])
		
		# variance parameters for observations
		tau.y <- pow(sigma.y,-2)
		sigma.y ~ dunif(0,10)

	# Creating B[site[]]: plot-level predictors
			for (k in 1:K){
				for (j in 1:J){
				B[j,k] <- xi[k]*B.raw[j,k]
				}
			xi[k] ~ dunif(0,10)	
			}
			
	# hyperparameters for plot-level predictors drawn from multivariate normal with corr errors
			for(j in 1:J){
			B.raw[j,1:K] ~ dmnorm(B.raw.hat[j,],Tau.B.raw[,]) # parameters are corrleated random around expected
				for(k in 1:K){
				B.raw.hat[j,k] <- inprod(G.raw[k,],U[j,])
				E.B[j,k] <- B.raw[j,k] - B.raw.hat[j,k]
				}
			}
	
	# incorporating group-level predictor data
		for (k in 1:K){
			for (l in 1:L){
				G[k,l] <- xi[k]*G.raw[k,l]
				G.raw[k,l] ~ dnorm(0,0.01)
			}
		}
		
	# inverse Wishart covariance matrix
	Tau.B.raw[1:K,1:K] ~ dwish(W[,],df)
	df <- K+1
	
	## Calculating standard deviations of all variables
	s.y <- sd(err.y[])
	s.N <- sd(B[,2])
	s.P <- sd(B[,3])
	s.pH <- sd(B[,4])
	s.C <- sd(B[,5])
	s.ndep <- abs(G[1,2])*sd(U[,2])
	s.map <- abs(G[1,3])*sd(U[,3])
	s.pet <- abs(G[1,4])*sd(U[,4])
	s.elev <- abs(G[1,5])*sd(U[,5])
	
	# estimating covariance (rho) among group predictors
	Sigma.B.raw[1:K,1:K] <- inverse(Tau.B.raw[,])
	for(k in 1:K){
		for (k.prime in 1:K){
			rho.B[k,k.prime] <- Sigma.B.raw[k,k.prime]/sqrt(Sigma.B.raw[k,k]*Sigma.B.raw[k.prime,k.prime])
		}
		sigma.B[k] <- abs(xi[k])*sqrt(Sigma.B.raw[k,k])	# scaled SD of parameter k
	}	
}		
	", fill=TRUE)
sink()
file.show('Ndep-model-matrix-corr-nomat.txt')




### Preparing data and initial parameter inputs

#data
y <- as.numeric(log(data$live_mass))
site <- as.numeric(factor(data$site_code))
n <- length(y)
J <- length(unique(data$site_code))	# number groups
# group parameter matrix from scaled site climate / Ndep predictors
U <- matrix(data=c(rep(1,J),ndep.z,map.z,pet.z,elev.z),nrow=J)	
L <- dim(U)[2]
# plot parameter matrix from scaled plot soil predictors
X <- matrix(data=c(rep(1,n),N.z,P.z,pH.z,sC.z),nrow=n)
K <- dim(X)[2]
W <- diag(K)

# bundled data for JAGS
jags.data <- list(y=y,site=site,n=n,J=J,U=U,L=L,X=X,K=K,W=W)

# need to create an inverse Rwishart function to generate initial values
rwish <- function(df,sigma){
	rWishart(1,df,sigma)[,,1]
}

# values for chain initiation
inits <- function(){
	list(
	sigma.y=runif(1),
	xi = runif(K),
	B.raw = array(rnorm(J*K),c(J,K)), # matrix of initial slopes 
	G.raw = array(rnorm(K*L),c(K,L)), # matrix of initial group-level predictors
	Tau.B.raw = rwish(K,W))
}

# parameters to track
params <- c('B','G','sigma.y','err.y','fit','fit.new','rho.B','E.B','B.raw','s.y','s.N','s.P','s.pH','s.C','s.ndep','s.map','s.pet','s.elev') 

# MCMC parameters – these are FINAL parameters (used much shorter for tuning)
# number iterations
ni <- 100000
# burn-in period to discard
nb <- 5000
# number of steps between samples in MCMC chain
nt <- 200
# number of chains
nc <- 3

# number of total draws from model for inference
(reps <-(ni-nb)/nt)

# set it off
out.model <- jags(data=jags.data, inits=inits, parameters=params, 
	model.file="Ndep-model-matrix-corr-nomat.txt", n.thin=nt, n.chains=nc, 
	n.burnin=nb, n.iter=ni)

print(out.model,dig=3)

References

R Core Team. 2013. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL http://www.R-project.org/.

Su, Y. and M. Yajima. (2013). R2jags: A Package for Running jags from R. R package version 0.03-11.  http://CRAN.R-project.org/package=R2jags