# Tumor
# han: Xt1, kan: Xt2, kes: Xt3, hyo:Xt4
Xt1 <- c(1,1,1,0,0,0,0,0,0,0,0,0)
Xt2 <- c(0,0,0,1,1,1,0,0,0,0,0,0)
Xt3 <- c(0,0,0,0,0,0,1,1,1,0,0,0)
Xt4 <- c(0,0,0,0,0,0,0,0,0,1,1,1)

# Region
# te: Xr1, to: Xr2, do: Xr3
Xr1 <- c(1,0,0,1,0,0,1,0,0,1,0,0)
Xr2 <- c(0,1,0,0,1,0,0,1,0,0,1,0)
Xr3 <- c(0,0,1,0,0,1,0,0,1,0,0,1)

# Number of incidence
Y   <- c(10, 22, 2, 28, 11, 17, 73, 19, 33, 115, 16, 54)

# Raw Data
RawData <- data.frame(cbind
           (Y,
            Xt1, Xt2, Xt3, Xt4,
            Xr1, Xr2, Xr3)
           )

# List
data <- list(Y   = RawData$Y,
             Xt1 = RawData$Xt1,
             Xt2 = RawData$Xt2,
             Xt3 = RawData$Xt3,
             Xt4 = RawData$Xt4,
             Xr1 = RawData$Xr1,
             Xr2 = RawData$Xr2,
             Xr3 = RawData$Xr3,
             N   = nrow(RawData)
             )
stancode <- "
data {
  int N;
  int Y[N];
  int Xt1[N];
  int Xt2[N];
  int Xt3[N];
  int Xt4[N];
  int Xr1[N];
  int Xr2[N];
  int Xr3[N];

}

parameters {
  real mu;
  real a1;
  real a2;
  real a3;
  real b1;
  real b2;
  real a1b1;
  real a1b2;
  real a2b1;
  real a2b2;
  real a3b1;
  real a3b2;
}

transformed parameters {
  real Lmbd[N];
  real LogLmbd[N];
  real a4;
  real b3;
  real a4b1;
  real a4b2;
  real a1b3;
  real a2b3;
  real a3b3;
  real a4b3;

  a4 = -(a1 + a2 + a3);
  b3 = -(b1 + b2);
  a4b1 = -(a1b1 + a2b1 + a3b1);
  a4b2 = -(a1b2 + a2b2 + a3b2);
  a1b3 = -(a1b1 + a1b2);
  a2b3 = -(a2b1 + a2b2);
  a3b3 = -(a3b1 + a3b2);
  a4b3 = -(a4b1 + a4b2);

  for (i in 1:N)
     LogLmbd[i] = mu + a1 * Xt1[i] + a2 * Xt2[i] + a3 * Xt3[i] + a4 * Xt4[i]
                     + b1 * Xr1[i] + b2 * Xr2[i] + b3 * Xr3[i]
                     + a1b1 * Xt1[i] * Xr1[i]
                     + a1b2 * Xt1[i] * Xr2[i]
                     + a1b3 * Xt1[i] * Xr3[i]
                     + a2b1 * Xt2[i] * Xr1[i]
                     + a2b2 * Xt2[i] * Xr2[i]
                     + a2b3 * Xt2[i] * Xr3[i]
                     + a3b1 * Xt3[i] * Xr1[i]
                     + a3b2 * Xt3[i] * Xr2[i]
                     + a3b3 * Xt3[i] * Xr3[i]
                     + a4b1 * Xt4[i] * Xr1[i]
                     + a4b2 * Xt4[i] * Xr2[i]
                     + a4b3 * Xt4[i] * Xr3[i];

  for (i in 1:N)
    Lmbd[i] = exp(LogLmbd[i]);
}

model {
  for (i in 1:N)
     Y[i] ~ poisson (Lmbd[i]);
}
"
parameters <- c("mu", 
                "a1", "a2", "a3", "a4",
                "b1", "b2", "b3",
                "a1b1", "a1b2", "a1b3",
                "a2b1", "a2b2", "a2b3",
                "a3b1", "a3b2", "a3b3",
                "a4b1", "a4b2", "a4b3")

library(rstan)

fit <- stan(model_code = stancode, 
            data   = data,
            pars   = parameters,
            iter   = 1000,
            chains = 4,
            thin   = 1,
            warmup = 500,
            seed   = 777
)
print(fit, digits_summary = 3)
## Inference for Stan model: cede7ce6929998090f76d95420ea7c5f.
## 4 chains, each with iter=1000; warmup=500; thin=1; 
## post-warmup draws per chain=500, total post-warmup draws=2000.
## 
##          mean se_mean    sd     2.5%      25%      50%      75%    97.5%
## mu      3.035   0.003 0.087    2.849    2.978    3.040    3.096    3.192
## a1     -1.099   0.008 0.213   -1.556   -1.236   -1.081   -0.951   -0.721
## a2     -0.211   0.004 0.135   -0.479   -0.301   -0.213   -0.123    0.057
## a3      0.524   0.003 0.114    0.306    0.447    0.524    0.598    0.748
## a4      0.786   0.003 0.109    0.570    0.712    0.783    0.857    1.006
## b1      0.612   0.003 0.107    0.413    0.539    0.608    0.683    0.823
## b2     -0.265   0.004 0.113   -0.493   -0.338   -0.264   -0.191   -0.043
## b3     -0.346   0.005 0.148   -0.648   -0.443   -0.340   -0.242   -0.093
## a1b1   -0.303   0.008 0.260   -0.774   -0.486   -0.313   -0.133    0.220
## a1b2    1.399   0.009 0.241    0.959    1.239    1.389    1.550    1.925
## a1b3   -1.096   0.014 0.395   -1.888   -1.354   -1.059   -0.823   -0.397
## a2b1   -0.118   0.005 0.169   -0.458   -0.230   -0.123   -0.006    0.226
## a2b2   -0.202   0.007 0.201   -0.595   -0.343   -0.197   -0.063    0.182
## a2b3    0.320   0.006 0.214   -0.117    0.181    0.327    0.459    0.724
## a3b1    0.115   0.004 0.133   -0.153    0.027    0.115    0.207    0.367
## a3b2   -0.381   0.005 0.164   -0.706   -0.492   -0.377   -0.269   -0.071
## a3b3    0.266   0.005 0.180   -0.071    0.141    0.264    0.384    0.618
## a4b1    0.306   0.003 0.133    0.059    0.213    0.306    0.398    0.565
## a4b2   -0.815   0.004 0.170   -1.170   -0.926   -0.809   -0.699   -0.508
## a4b3    0.510   0.006 0.176    0.185    0.388    0.505    0.618    0.877
## lp__ 1144.084   0.087 2.492 1138.538 1142.555 1144.436 1145.945 1147.949
##      n_eff  Rhat
## mu     976 1.004
## a1     649 1.007
## a2    1199 1.001
## a3    1067 1.007
## a4    1574 1.000
## b1    1196 1.000
## b2    1013 1.003
## b3     863 1.004
## a1b1   970 1.002
## a1b2   710 1.005
## a1b3   755 1.006
## a2b1  1238 1.001
## a2b2   929 1.002
## a2b3  1349 1.002
## a3b1  1305 1.001
## a3b2  1316 0.999
## a3b3  1229 1.001
## a4b1  2000 1.000
## a4b2  2000 1.001
## a4b3   940 1.003
## lp__   824 1.001
## 
## Samples were drawn using NUTS(diag_e) at Tue Apr 23 09:41:26 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).