data{
 int<lower=1> I;               //number of areas      
 int<lower=1> X;               //number of age groups 
 int<lower=1> S;               //number of groups     

 matrix[I,X] deaths[S];        //observed deaths
 matrix[I,X] pop[S];           //observed population
 matrix[I,I] A_sp;             //spatial adjacency matrix
 matrix[I,I] A_age;            //random walk(1) adjacency matrix
}


transformed data {
 matrix[I,I] D_sp;          // diagonal matrix with d_ii = sum(c_ij)
 matrix[I,I] D_age;         // diagonal matrix with d_ii = sum(c_ij)
 matrix[2,2] Q;             // scale matrix for Wishart
 vector[I] zeros_sp;        // zeros for mean of MVN
 vector[X] zeros_age;       // zeros for mean of MVN

 for (i in 1:I){
  zeros_sp[i] <- 0.0;
  for (j in 1:I)
   D_sp[i,j] <- if_else(i==j, sum(row(A_sp, i)), 0.0);
 }
 for (i in 1:X){
  zeros_age[i] <- 0.0;
  for (j in 1:X)
   D_age[i,j] <- if_else(i==j, sum(row(A_age, i)), 0.0);
 }
 for (i in 1:2)
  for (j in 1:2)
   Q[i,j] <- if_else(i==j, 1.0, 0.0);
}


parameters{
 vector[S] b0;                 // fixed gender effect
 vector[X] f1[S];              // unconstrained age-group random effect
 vector[I] f2[S];              // unconstrained spatial random effect
 vector[X] f3[S];              // unconstrained Gamma(1,1) parameters
 real<lower=-1,upper=1> p_sp;  // strength of spatial correlation
 real<lower=-1,upper=1> p_age; // strength of age-group correlation
 matrix[2,2] Lambda_sp         // spatial MVCAR scale matrix
 matrix[2,2] Lambda_age        // age-group MVCAR scale matrix
}


derived parameters{
 real mrate[S,I,X];            // estimated mortality rate
 vector[X] b1[S];              // sum-to-0 constrained age-group random effect
 vector[I] b2[S];              // sum-to-0 constrained spatial random effect        
 vector[X] b3[S];              // sum-to-1 constrained Gamma(1,1) interaction-parameters

 for (s in 1:S){
  b1[s] <- f1[s]-mean(f1[s]);  // b1 ~ MVCAR
  b2[s] <- f2[s]-mean(f2[s]);  // b2 ~ MVCAR
  b3[s] <- f3[s]/sum(f3[s]);
 }
 for (s in 1:S)
  for (i in 1:I)
   for (x in 1:X)
    mrate[s,i,x] <- exp( b0[s] + b1[s,x] + b2[s,i]*b3[s,x] ); // required for life expectancy calculations   
}


model{
 matrix[2N,2N] prec_sp;          // precision matrix for spatial MVCAR
 matrix[2X,2X] prec_age;         // precision matrix for RW(1)   MVCAR
 
 Lambda_sp  ~ wishart(2,Q);
 Lambda_age ~ wishart(2,Q);

 // Kronecker products 
 prec_age  <- Kronecker_product(D_age - p_age * A_age) , Lambda_age);
 prec_sp   <- Kronecker_product(D_sp  - p_sp  * A_sp ) , Lambda_sp );

 f1 ~ multi_normal_prec(zeros_age, prec_age);
 f2 ~ multi_normal_prec(zeros_sp , prec_sp );
 f3 ~ gamma(1,1);

 // likelihood
 for (s in 1:S)
  deaths[s] ~ dpois( pop[s].*mrate[s] );     	
}