##--------------------------------------------------------------------------------------------------------
## SCRIPT : Meta-analyse
##
## Authors : Collectif BioBayes
## Last update : 2019-03-07
## R version 3.5.2 (2018-12-20) -- "Eggshell Igloo"
## Copyright (C) 2018 The R Foundation for Statistical Computing
## Platform: x86_64-w64-mingw32/x64 (64-bit)
##--------------------------------------------------------------------------------------------------------

lapply(c("ggplot2", "ggthemes", "rstan"), library, character.only=TRUE)

theme_set(theme_bw(base_size = 20))

### clean up
rm(list = ls())

setwd(WorkDir <- getwd())
DataDir <- paste(WorkDir, "data", sep = "/")
OutDir <- paste(WorkDir, "output", sep = "/")

### charger les données
don <- list(K = 14,
            rm = c(1, 9, 2, 1, 10, 1, 1, 90, 0, 6, 1, 5, 4, 4),
            rc = c(2, 23, 7, 1, 8, 9, 3, 118, 1, 11, 7, 13, 8, 17),
            nm = c(40, 135, 200, 48, 150, 59, 25, 1159, 22, 76, 27, 23, 130, 107),
            nc = c(36, 135, 200, 46, 148, 56, 23, 1157, 21, 75, 27, 33, 122, 108)
            )


### Modélisation avec Stan
# rstan (Version 2.18.2, GitRev: 2e1f913d3ca3)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

meta <- '
 data {
  int<lower = 1> K;
  int<lower = 1> nm[K];
  int<lower = 0> rm[K];
  int<lower = 1> nc[K];
  int<lower = 0> rc[K];
  real<lower = 0.0> hyperprior_scale_mu;
  real hyperprior_location_mu;
 }

 parameters {
  real unscaled_mu;
  real<lower = 0.0, upper = 10.0> sigma;
  vector[K] eps;
  vector<lower = 0.0, upper = 1.0>[K] pc;
 }

 transformed parameters {
  real mu;
  real odds_ratio_pop;
  vector[K] delta;
  vector<lower = 0.0, upper = 1.0>[K] pm;
  // equivalent to mu ~ normal(hyperprior_location_mu, hyperprior_scale_mu)
  mu = hyperprior_location_mu + hyperprior_scale_mu * unscaled_mu;
  // equivalent to delta ~ normal(mu, sigma)
  delta = rep_vector(mu, K) + sigma * eps;
  pm = inv_logit(logit(pc) + delta);
  odds_ratio_pop = exp(mu);
 }
 
 model {
  // no explicit prior on pc or sigma: 
  // Stan will take use independent uniform priors on the range defined in the parameter block
  pc ~ beta(1.0, 1.0);
  eps ~ normal(0.0, 1.0);
  unscaled_mu ~ normal(0.0, 1.0);
  rc ~ binomial(nc, pc);
  rm ~ binomial(nm, pm);
 }

 generated quantities {
  real log_new_odds_ratio;
  real new_odds_ratio;
  vector[2 * K] log_lik;
  log_new_odds_ratio = normal_rng(mu, sigma);
  new_odds_ratio = exp(log_new_odds_ratio);
  for(i in 1:K) {
   log_lik[i] = binomial_lpmf(rc[i]| nc[i], pc[i]);
   log_lik[i + K] = binomial_lpmf(rm[i]| nm[i], pm[i]);
  }
 }

'

## Compilation du modèle 
model_meta <- stan_model(model_code = meta,
                         model_name = "Meta-analyse"
                         )

## Echantillonnage
n_chains <- 3
n_iter <- 500
n_warm <- 200
n_thin <- 1

standata <- don
standata$hyperprior_location_mu <- 0
standata$hyperprior_scale_mu <- log(2) / 2 # apriori mg treatment can halve or double mortality

### 
fitted_model_1 <- sampling(object = model_meta,
                           data = standata, 
                           pars = c("pc", "pm", "mu", "sigma", "odds_ratio_pop", "new_odds_ratio", "log_lik"),
                           chains = n_chains, 
                           iter = n_iter, 
                           warmup = n_warm, 
                           thin = n_thin
                           )

print(fitted_model_1, digits_summary = 3)
stan_rhat(fitted_model_1)
get_elapsed_time(fitted_model_1)

### regarder rapidement quelques paramètres
pairs(x = fitted_model_1, 
      pars = c("sigma", "mu", "odds_ratio_pop", "lp__")
      )

### regarder la trace des parametres
plot(fitted_model_1, 
     plotfun = "trace", 
     pars = c("sigma", "mu", "odds_ratio_pop", "lp__"), 
     inc_warmup = TRUE
     )

plot(fitted_model_1, 
     plotfun = "trace", 
     pars = c("sigma", "mu", "odds_ratio_pop", "lp__"), 
     inc_warmup = TRUE
     )

### regarder la distribution a posteriori
plot(fitted_model_1, 
     plotfun = "hist", 
     pars = c("sigma", "mu", "odds_ratio_pop", "lp__"), 
     inc_warmup = FALSE
     )

### calcul du WAIC et du critere Leave-One-Out
loo::waic(loo::extract_log_lik(fitted_model_1))
loo::loo(loo::extract_log_lik(fitted_model_1))

### calculer la probabilite que 
## le parametre odds_ratio_pop > 1 ?
## le parametre new_odds_ratio > 1 ?
# que conclure?

### refaire l'analyse en utilisant le prior informatif de Higgins & Spiegelhalter 2002
## que faut-il changer?

### refaire l'analyse en intégrant les données du mégatrial de 1994
## pour cette 15eme etude, motivee par les resultats de la meta-analyse de nombreux patients
## ont ete enrolles
#  control : rc = 2103 / nc = 29039
#  magnesium : rm = 2216 / nm = 29011