##--------------------------------------------------------------------------------------------------------
## SCRIPT : TD Teneurs en protéines avec Stan
## Authors : ## Authors : Collectif BioBayes
## Last update : 2019-03-06
## R version 3.5.1 (2019-03-06) -- "Feather Spray"
## Copyright (C) 2018 The R Foundation for Statistical Computing
## Platform: x86_64-w64-mingw32/x64 (64-bit)
##--------------------------------------------------------------------------------------------------------

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

# rstan (Version 2.14.1, packaged: 2017-04-19 05:03:57 UTC, GitRev: 2e1f913d3ca3)
# For execution on a local, multicore CPU with excess RAM we recommend calling
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

rm(list = ls())

###Paramétrage graphique pour ggplot
theme_set(theme_bw(base_size = 20))


###Les données 
Y=c(11.9,10.3,12.4,12.1,12.4,12.7,11.5,13.6,9.9,11.1,10.6,12,12.8,11.6,11.2,10.3, 
    12.9,10.3,11.5,14.1,12.6,12.7,11.9,9.2,12.5,11.9,11.3,11,10,15.9,10.6,11.6, 
    11.4,10.9,14.6,10.7,9.4,13,12.5,12.9,11.1,12.9,13.2)
INN=c(1,0.8,0.8,1,0.7,0.9,1,0.9,0.6,0.7,0.7,0.7,0.9,0.8,0.8,0.6,0.9,0.69,0.72,1.1,0.97,0.98,0.87,
      0.62,0.73,1.09,0.76,0.86,0.75,1.03,0.75,0.86,0.82,0.77,0.88,0.84,0.72,0.81,0.74,0.94,0.78,0.92,1.07)
SPAD=c(45.5,42.6,47.8,48.6,48.2,48.8,48.9,49.9,48.6,44.9,50.2,46.1,51.6,50.9,40.4,
       47.6,49.9,37.2,44.87,48,46.03,45.8,46.87,38.67,42.67,44.23,39.67,39.4,39.43,44.73,38.27,44.6,45.5,46.1,
       38,43.6,43.5,44,49.9,51,48.8,48.7,49.6)
###Covariables centrées-réduites
INN_cr= (INN-mean(INN))/sd(INN)
SPAD_cr= (SPAD-mean(SPAD))/sd(SPAD)




################################
#Modèle à 1 variable explicative 
################################

### code du modele
teneur_univarie <- "
data {
  int<lower=0> N;       //Nombre de données
  vector[N] covar_cr;   //Vecteur du prédicteur centré-réduit
  vector[N] Y;          //Variable réponse
  real m;               //Moyenne empirique du prédicteur (quand non centré réduit)
  real<lower=0> s;      //Ecart-type empirique du prédicteur (quand non centrée réduit)
  vector[2] mu;         //Moyennes a priori des coefficients alpha
  matrix[2,2] W;        //Matrice de variance-covariance a priori des coefficients alpha
}

parameters {
  vector[2] alpha;                  
  real<lower=0, upper=100> sigma;   //Loi a priori Uniforme sur [0,100]
}

model{
  Y ~ normal(alpha[1]+alpha[2]*covar_cr,sigma);
  alpha ~ multi_normal(mu,W);     //Loi a priori normales 
}

generated quantities { //Quantités inconnues d'intérêt non définies dans le modèle
    vector[2] alpha_nonc;
    vector[N] loglike;

    alpha_nonc[1] = alpha[1]-alpha[2]*(m/s);
    alpha_nonc[2] = alpha[2]/s;
    
    for(i in 1:N){ 
      loglike[i] = normal_lpdf(Y[i] | alpha[1]+alpha[2]*covar_cr[i], sigma); //Log vraisemblance associée à chaque observation
    }
}

"

###Ajustement modèle à 1 variable explicative = INN 
stanfit_Inn <- stan(model_name = "TeneurUnivarie_Inn", 
                 model_code = teneur_univarie,
                 data =list(Y=Y,covar_cr=INN_cr,N=as.integer(43),m=mean(INN),s=sd(INN),mu=c(0,0),W=matrix(c(1000,0,0,1000),nr=2,ncol=2)),
                 pars = c("alpha_nonc", "sigma", "loglike"),
                 chains = 3, 
                 iter = 4200,  #warmup is included
                 warmup = 200,
                 thin = 1,
                 control = list(max_treedepth=15),
                 init = list(list(alpha=c(0,0), sigma=1),
                             list(alpha = c(5,5), sigma = 0.1),
                             list(alpha = c(-5,-5), sigma = 10)))


#Représentations graphiques des chaînes 
stan_trace(stanfit_Inn,pars= c("alpha_nonc","sigma"), inc_warmup=TRUE)
stan_trace(stanfit_Inn,pars= c("alpha_nonc","sigma"), inc_warmup=FALSE)
#Utiliser le diagnostic de convergence de Gelman et Rubin
gelman.plot(As.mcmc.list(stanfit_Inn,pars= c("alpha_nonc","sigma"), inc_warmup=TRUE))
#Analyser les auto-corrélations des chaînes
autocorr.plot(As.mcmc.list(stanfit_Inn,pars= c("alpha_nonc","sigma")))
#Densités a posteriori
stan_dens(stanfit_Inn,pars= c("alpha_nonc","sigma"))
#Statistiques a posteriori
print(stanfit_Inn, digits_summary = 2)
#Temps de calcul
get_elapsed_time(stanfit_Inn)

#Calcul du critère WAIC
loo::waic(rstan::extract(stanfit_Inn, 'loglike')$loglike)

#Calcul du critère DIC
DIC<- function(stanfit,chain,covar1,covar2=NULL,modele){
  M<- as.matrix(chain$loglike)
  D<- -2*apply(M,FUN=sum,MARGIN=1)
  Dbar<- mean(D)

  if(modele=="univarie"){
    Y_postmean = get_posterior_mean(stanfit,pars= c("alpha_nonc[1]"))[,4]+get_posterior_mean(stanfit,pars= c("alpha_nonc[2]"))[,4]*covar1
  }
  if(modele=="bivarie"){
    Y_postmean = get_posterior_mean(stanfit,pars= c("alpha_nonc[1]"))[,4]+get_posterior_mean(stanfit,pars= c("alpha_nonc[2]"))[,4]*covar1 + get_posterior_mean(stanfit,pars= c("alpha_nonc[3]"))[,4]*covar2 
  }
  if(modele=="sanscov"){
    Y_postmean = get_posterior_mean(stanfit,pars= c("mu"))[,4]
  }
  loglik_postmean = dnorm(Y,Y_postmean, get_posterior_mean(stanfit,pars= c("sigma"))[,4], log=TRUE)
  Dhat = -2*sum(loglik_postmean)
  pD<- Dbar - Dhat 
  return(list("Dbar"=Dbar, "pD"=pD, "DIC"=Dbar + pD))
}

chain_INN<- extract(stanfit_Inn)
DIC(stanfit_Inn,chain_INN, covar1=INN,modele="univarie")




###Ajustement modèle à 1 variable explicative = SPAD
stanfit_SPAD <- stan(model_name = "TeneurUnivarie_SPAD", 
                    model_code = teneur_univarie,
                    data =list(Y=Y,covar_cr=SPAD_cr,N=as.integer(43),m=mean(SPAD),s=sd(SPAD),mu=c(0,0),W=matrix(c(1000,0,0,1000),nr=2,ncol=2)),
                    pars = c("alpha_nonc", "sigma", "loglike"),
                    chains = 3, 
                    iter = 4200,  #warmup is included
                    warmup = 200,
                    thin = 1,
                    control = list(max_treedepth=15),
                    init = list(list(alpha=c(0,0), sigma=1),
                                list(alpha = c(5,5), sigma = 0.1),
                                list(alpha = c(-5,-5), sigma = 10)))



#Représentations graphiques des chaînes 
stan_trace(stanfit_SPAD,pars= c("alpha_nonc","sigma"), inc_warmup=TRUE)
stan_trace(stanfit_SPAD,pars= c("alpha_nonc","sigma"), inc_warmup=FALSE)
#Utiliser le diagnostic de convergence de Gelman et Rubin
gelman.plot(As.mcmc.list(stanfit_SPAD,pars= c("alpha_nonc","sigma"), inc_warmup=TRUE))
#Analyser les auto-corrélations des chaînes
autocorr.plot(As.mcmc.list(stanfit_SPAD,pars= c("alpha_nonc","sigma")))
#Densités a posteriori
stan_dens(stanfit_SPAD,pars= c("alpha_nonc","sigma"))
#Statistiques a posteriori
print(stanfit_SPAD, digits_summary = 2)
#Temps de calcul
get_elapsed_time(stanfit_SPAD)


#Calcul du critère WAIC
loo::waic(rstan::extract(stanfit_SPAD, 'loglike')$loglike)

#Calcul du critère DIC 
chain_SPAD<- extract(stanfit_SPAD)
DIC(stanfit_SPAD, chain_SPAD, covar1=SPAD,modele="univarie")

  


##################################
#Modèle à 2 variables explicatives 
##################################

### code du modèle
teneur_bivarie <- "

data {
  int<lower=0> N;    //Nombre de données
  vector[N] INN_cr;   //Vecteur du prédicteur INN centré réduit
  vector[N] SPAD_cr;   //Vecteur du prédicteur SPAD centré réduit
  vector[N] Y;       //Variable réponse
  real m_INN;            //Moyenne empirique de INN
  real<lower=0> s_INN;   //Ecart-type empirique de INN
  real m_SPAD;            //Moyenne empirique de SPAD
  real<lower=0> s_SPAD;   //Ecart-type empirique de SPAD
  vector[3] mu;      //Moyennes a priori des coefficients alpha
  matrix[3,3] W;     //Matrice de variance-covariance a priori des coefficients alpha
}


parameters {
  vector[3] alpha;                 
  real<lower=0, upper=100> sigma;   //Loi a priori Uniforme sur [0,100]
}

model{
  Y ~ normal(alpha[1]+alpha[2]*SPAD_cr+alpha[3]*INN_cr,sigma);
  alpha ~ multi_normal(mu,W);     //Loi a priori normales 
}

generated quantities { //Quantités inconnues d''intérêt non définies dans le modèle
vector[3] alpha_nonc;
vector[N] loglike;

alpha_nonc[1] = alpha[1]-alpha[2]*(m_SPAD/s_SPAD)-alpha[3]*(m_INN/s_INN);
alpha_nonc[2] = alpha[2]/s_SPAD;
alpha_nonc[3] = alpha[3]/s_INN;

for(i in 1:N){ 
  loglike[i] = normal_lpdf(Y[i] | alpha[1]+alpha[2]*SPAD_cr[i]+alpha[3]*INN_cr[i], sigma);
}
}

"

###Ajustement modèle bivarié
stanfit_bivarie <- stan(model_name = "TeneurBivarie", 
                     model_code = teneur_bivarie,
                     data =list(Y=Y,SPAD_cr=SPAD_cr,INN_cr=INN_cr,N=as.integer(43),m_SPAD=mean(SPAD),s_SPAD=sd(SPAD),m_INN=mean(INN),s_INN=sd(INN),mu=c(0,0,0),W=matrix(c(1000,0,0,0,1000,0,0,0,1000),nr=3,ncol=3,byrow=TRUE)),
                     pars = c("alpha_nonc", "sigma", "loglike"),
                     chains = 3, 
                     iter = 4200,  #warmup is included
                     warmup = 200,
                     thin = 1,
                     control = list(max_treedepth=15),
                     init = list(list(alpha=c(0,0,0), sigma=1),
                                 list(alpha = c(5,5,5), sigma = 0.1),
                                 list(alpha = c(-5,-5,-5), sigma = 10)))



#Représentations graphiques des chaînes 
stan_trace(stanfit_bivarie,pars= c("alpha_nonc","sigma"), inc_warmup=TRUE)
stan_trace(stanfit_bivarie,pars= c("alpha_nonc","sigma"), inc_warmup=FALSE)
#Utiliser le diagnostic de convergence de Gelman et Rubin
gelman.plot(As.mcmc.list(stanfit_bivarie,pars= c("alpha_nonc","sigma"), inc_warmup=TRUE))
#Analyser les auto-corrélations des chaînes
autocorr.plot(As.mcmc.list(stanfit_bivarie,pars= c("alpha_nonc","sigma")))
#Densités a posteriori
stan_dens(stanfit_bivarie,pars= c("alpha_nonc","sigma"))
print(stanfit_bivarie, digits_summary = 2)
#Temps de calcul
get_elapsed_time(stanfit_bivarie)


#Calcul du critère WAIC
loo::waic(rstan::extract(stanfit_bivarie, 'loglike')$loglike)

#Calcul du critère DIC 
chain_bivarie<- extract(stanfit_bivarie)
DIC(stanfit_bivarie, chain_bivarie, covar1=SPAD,covar2= INN, modele="bivarie")




################################
#Modèle sans variable explicative
################################

### code du modele
teneur_sanscov <- "

data {
  int<lower=0> N;    //Nombre de données
  vector[N] Y;       //Variable réponse
  real a;           //Moyennes a priori du coefficient mu
  real b;           //Ecart-type a priori du coefficient mu
}

parameters {
  real mu;                  
  real<lower=0, upper=100> sigma;   //Loi a priori Uniforme sur [0,100]
}

model{
  Y ~ normal(mu,sigma);
  mu ~ normal(a,b);     //Loi a priori normale 
}

generated quantities { //Quantités inconnues d''intérêt non définies dans le modèle
  vector[N] loglike;

  for(i in 1:N){ 
    loglike[i] = normal_lpdf(Y[i] | mu, sigma);
  }
}

"

###Ajustement modèle sans covariable
stanfit_sanscov <- stan(model_name = "TeneurUnivarie_sanscov", 
                    model_code = teneur_sanscov,
                    data =list(Y=Y,N=as.integer(43),a=0,b=1000),
                    pars = c("mu", "sigma", "loglike"),
                    chains = 3, 
                    iter = 4200,  #warmup is included
                    warmup = 200,
                    thin = 1,
                    control = list(max_treedepth=15),
                    init = list(list(mu=0, sigma=1),
                                list(mu=5, sigma = 0.1),
                                list(mu=-5, sigma = 10)))




#Représentations graphiques des chaînes 
stan_trace(stanfit_sanscov,pars= c("mu","sigma"), inc_warmup=FALSE)
#Utiliser le diagnostic de convergence de Gelman et Rubin
gelman.plot(As.mcmc.list(stanfit_sanscov,pars= c("mu","sigma"), inc_warmup=TRUE))
#Analyser les auto-corrélations des chaînes
autocorr.plot(As.mcmc.list(stanfit_sanscov,pars= c("mu","sigma")))
#Densités a posteriori
stan_dens(stanfit_sanscov,pars= c("mu","sigma"))
print(stanfit_sanscov, digits_summary = 2)
#Temps de calcul
get_elapsed_time(stanfit_sanscov)


#Calcul du critère WAIC
loo::waic(rstan::extract(stanfit_sanscov, 'loglike')$loglike)

#Calcul du critère DIC 
chain_sanscov<- extract(stanfit_sanscov)
DIC(stanfit_sanscov, chain_sanscov, covar1=NULL,covar2= NULL, modele="sanscov")




########################################################################
#Prediction a posteriori selon le modèle à 1 variable explicative = INN
###########################################################################
### code du modele
teneur_pred <- "
data {
  int<lower=0> N;       //Nombre de données
  vector[N] covar_cr;   //Vecteur du prédicteur centré-réduit
  vector[N] Y;          //Variable réponse
  real covar_pred;      //Valeur du prédicteur associé à la variable à prédire
  real m;               //Moyenne empirique du prédicteur (quand non centré réduit)
  real<lower=0> s;      //Ecart-type empirique du prédicteur (quand non centrée réduit)
  vector[2] mu;         //Moyennes a priori des coefficients alpha
  matrix[2,2] W;        //Matrice de variance-covariance a priori des coefficients alpha
}

parameters {
  vector[2] alpha;                  
  real<lower=0, upper=100> sigma;   //Loi a priori Uniforme sur [0,100]
}

model{
  Y ~ normal(alpha[1]+alpha[2]*covar_cr,sigma);
  alpha ~ multi_normal(mu,W);     //Loi a priori normales 
}

generated quantities { //Quantités inconnues d'intérêt non définies dans le modèle
  real Y_pred;
  Y_pred = normal_rng(alpha[1]+alpha[2]*((covar_pred -m)/s),sigma);
}

"

###Ajustement modèle à 1 variable explicative = INN 
#Vérifier tout d'abord que les programmes fonctionnent en réalisant 5000 itérations MCMC
#Stocker les valeurs échantillonnées pour trois chaînes de Markov
stanfit_pred<- stan(model_name = "Teneur_pred", 
                    model_code = teneur_pred,
                    data =list(Y=Y,covar_cr=INN_cr,N=as.integer(43),m=mean(INN),s=sd(INN),covar_pred=0.95,mu=c(0,0),W=matrix(c(1000,0,0,1000),nr=2,ncol=2)),
                    pars = c("Y_pred"),
                    chains = 3, 
                    iter = 4200,  #warmup is included
                    warmup = 200,
                    thin = 1,
                    control = list(max_treedepth=15),
                    init = list(list(alpha=c(0,0), sigma=1),
                                list(alpha = c(5,5), sigma = 0.1),
                                list(alpha = c(-5,-5), sigma = 10)))


#Statistiques a posteriori
print(stanfit_pred, digits_summary = 2)

#Densités a posteriori
stan_dens(stanfit_pred,pars=c("Y_pred"))