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

	
	for (i in 1:Nareas){
		
		#### if M is a square matrix define Nsp (Number of spatial underlying patterns) as Ndiseases
		for(k in 1:Ndiseases){
			
			Theta[i,k]<-inprod2(tPhi[,i],M[,k])
		
		}
	
	}
  	

	#### Matrix of spatially correlated random effects:
	for(j in 1:Nsp){
		
		tPhi[j, 1:Nareas]~car.proper(ceros[],C[],adj[],num[],M.car[],1,gamma[j])
		
		gamma[j]~dunif(gamma.inf,gamma.sup)
	
	} 		
	

	for(i in 1:Nareas){ceros[i]<-0}
	
	gamma.inf<-min.bound(C[],adj[],num[],M.car[])
	
	gamma.sup<-max.bound(C[],adj[],num[],M.car[])

	

	#### M-matrix
	for (i in 1:Ndiseases){
		
		for (j in 1:Ndiseases){
	
			### This corresponds to the fixed effects model		
			M[i,j] ~ dflat() 

			### In case of implementing the random effects model comment the previous line
			### and uncomment the next one
			###M[i,j] ~ dnorm(0,prec)		
		}
	
	}

	

	### In case of implementing the random effects model uncomment the next two lines
	###prec<-pow(sdstruct,-2)
	
	###sdstruct~dunif(0,100)

	

	# Other priors
	
	for (k in 1:Ndiseases){
		
		mu[k] ~ dflat()
	
	}


}