model1 = "
functions{
  real frank_lpdf(matrix p_i, vector alpha, vector beta, real omega){
    real f1;
    real f2;
    int r;
    vector[rows(p_i)] f3;
    vector[rows(p_i)] f4;
    
    r = rows(p_i);
    
    f1 = beta_lpdf(col(p_i, 1)| alpha[1], beta[1]);
    f2 = beta_lpdf(col(p_i, 2)| alpha[2], beta[2]);
    for (i in 1:r){
      f3[i] = (omega*(1 - exp(-omega))*exp(-omega*(beta_cdf(p_i[i,1], alpha[1], beta[1]) + beta_cdf(p_i[i,2], alpha[2], beta[2]))));
      f4[i] = ((1 - exp(-omega)) - (1 - exp(-omega*beta_cdf(p_i[i,1], alpha[1], beta[1])))*(1 - exp(-omega*beta_cdf(p_i[i,2], alpha[2], beta[2]))));
    }
    return (f1 + f2 + sum(log((f3)./((f4).*(f4))))); 
  }
}
data{
  int N;
  int Nt;
  int Ns;
  
  int TP[N];
  int Dis[N];
  int TN[N];
  int NDis[N];
  int Study[N];
  int Test[N];
}
parameters{
  vector<lower=0, upper=1>[2] MU[Nt]; //mu_jk
  vector<lower=0, upper=1>[2] delta[Nt]; //delta_jk
  vector<lower=0, upper=1>[2] theta; //theta_j
  matrix<lower=0, upper=1>[Ns, 2] p_i[Nt]; //pi_ijk
  vector[Nt] omega; //omega_k
}
transformed parameters{
  vector<lower=0>[2] alpha[Nt]; //alpha_jk
  vector<lower=0>[2] beta[Nt]; //beta_jk
  
  
  for(k in 1:Nt){
    alpha[k] = (MU[k]).*((1 - (theta).*(delta[k]))./((theta).*(delta[k])));
    beta[k] = (1 - MU[k]).*((1 - (theta).*(delta[k]))./((theta).*(delta[k])));
  }
  
}
model{
  //Priors
  for (k in 1:Nt){
    MU[k] ~ uniform(0, 1);
    delta[k] ~ uniform(0, 1);
  }
  theta ~ uniform(0, 1);
  omega ~ normal(0, 5);
  
  for (k in 1:Nt)
    p_i[k] ~ frank(alpha[k], beta[k], omega[k]);
  
  for (n in 1:N){
    TP[n] ~ binomial(Dis[n], p_i[Test[n]][Study[n], 1]);
    TN[n] ~ binomial(NDis[n], p_i[Test[n]][Study[n], 2]);
  }
  
}
generated quantities{
  vector[2*N] loglik;
  
  for (n in 1:N)
    loglik[n] = binomial_lpmf(TN[n]| NDis[n], p_i[Test[n]][Study[n], 1]);
  
  for (n in (N+1):(2*N))
    loglik[n] = binomial_lpmf(TN[n-N]| NDis[n-N], p_i[Test[n-N]][Study[n-N], 2]);
  
}
"

N <- nrow(data)
Ns <- max(data$Study)
Nt <- max(data$Test)


datalist1 <- list(
    N = N,
    Ns = Ns,
    Nt = Nt,
    TP = data$TP,
    Dis = data$Dis,
    TN = data$TN,
    NDis = data$NDis,
    Test = data$Test,
    Study = data$Study)


fit1 <- stan(model_code = model1, 
             data=datalist1, 
             chains = 3,
             iter = 4000, 
             warmup = 1000, 
             thin = 5, 
             seed = 2,
             cores=3)