model {
	# Likelihood
	for (i in 1:Nareas) {
		for (k in 1:Ndiseases) {
			Y[i, k] ~ dpois(mu[i, k])
			log(mu[i, k]) <- log(E[i, k]) + alpha[k] + S[i,k]
  			RR[i,k] <- exp(alpha[k] + S[i,k])
		}
  	}

	for (i in 1:Nareas){
		for(k in 1:Ndiseases){
			S[i,k]<-inprod(tDelta[,i],tCholDisRotated[,k])
		}
	}
  	for(j in 1:Ndiseases){
		tDelta[j, 1:Nareas]~car.proper(ceros[],C[],adj[],num[],M[],1,gamma.ord[j])
		ordena[j]<-Ndiseases+1-j
		gamma.ord[j]<-ranked(gamma[],ordena[j])
		gamma[j]~dunif(gamma.inf,gamma.sup)
	} 		
	for(i in 1:Nareas){ceros[i]<-0}
	gamma.inf<-min.bound(C[],adj[],num[],M[])
	gamma.sup<-max.bound(C[],adj[],num[],M[])

	#Rotation of CholDis by means of an orthogonal transformation:
	#Orthogonal matrix definition:
	P[1,1]<-cos(theta12)*cos(theta13)
	P[1,2]<-sin(theta12)*cos(theta23)-cos(theta12)*sin(theta13)*sin(theta23)
	P[1,3]<-sin(theta12)*sin(theta23)+cos(theta12)*sin(theta23)*cos(theta23)
	P[2,1]<--sin(theta12)*cos(theta13)
	P[2,2]<-cos(theta12)*cos(theta23)+sin(theta12)*sin(theta13)*sin(theta23)
	P[2,3]<-cos(theta12)*sin(theta23)-sin(theta12)*sin(theta13)*cos(theta23)	
	P[3,1]<--sin(theta13)
	P[3,2]<--cos(theta13)*sin(theta23)
	P[3,3]<-cos(theta13)*cos(theta23)
	for(i in 1:3){for(j in 1:3){tCholDisRotated[j,i]<-inprod(CholDis[i,1:i],P[1:i,j])}}
  	theta12~dunif(0,1.5708)
	theta13~dunif(0,1.5708)
	theta23~dunif(0,1.5708)	

	#Between-diseases structure
	#Cholesky triangle of a general correlation matrix
	#diagonal
	CholDis[1,1]<-sigma[1]
	for(j in 2:Ndiseases){
		diag[j]<-pow(sigma[j],2)-inprod(CholDis[j,1:(j-1)],CholDis[j,1:(j-1)])
		CholDis[j,j]<-sqrt(abs(diag[j]))
	}

	#Lower triangle
	#First column
	for(i in 2:Ndiseases){CholDis[i,1]<-sigma[i]*roDis[i,1]}
	#Rest of columns 
	for(j in 2:(Ndiseases-1)){for(i in (j+1):Ndiseases){
		CholDis[i,j]<-(roDis[i,j]*sigma[i]*sigma[j]-inprod(CholDis[i,1:(j-1)],CholDis[j,1:(j-1)]))/CholDis[j,j]
	}}

	for(j in 1:Ndiseases){roDis[1,j]<-0}
	for(i in 2:Ndiseases){
		for(j in 1:(i-1)){roDis[i,j]~dunif(-1,1)}
		for(j in i:Ndiseases){roDis[i,j]<-0}
	}

	#Condition of being positive defined the covariance matrix
	uno<-1 
	uno~dbern(condition)
	condition<-step(sum(subcondition[2:Ndiseases])-(Ndiseases-1))
	for(j in 2:Ndiseases){subcondition[j]<-step(diag[j])}

	# Other priors
	for (k in 1:Ndiseases){
		alpha[k] ~ dflat()
		sigma[k] ~ dunif(0,10)
	}
}