# TD Magnésium - infarctus -WINBUGS

# à passer en JAGS

setwd("C:/Users/ewalker/Documents/docs_emily/FORMATIONS_RESEAUX/EC_BIOBAYES2019/TD4_metaanalyse")

library(rjags)
library(coda)


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))



cat("
    model{
      for(i in 1:K) {

      #Binomial structure      
      rc[i] ~ dbin(pc[i],nc[i]);
      rm[i] ~ dbin(pm[i],nm[i]);

      phi[i] <- logit(pc[i]);
      logit(pm[i]) <- phi[i] + delta[i];
      
      #Define log(odds ratio)
      delta[i] ~ dnorm(mu, precision);
      #Odds ratio
      Odds[i] <-exp(delta[i]); 
      
      #Random effects
      pc[i] ~ dunif(0,1)
      #Prior for pc
      }

      delta.new ~ dnorm(mu, precision);
      Odds.new<-exp(delta.new);
      Odds.pop<-exp(mu);

      #Predicted effect
      mu ~ dnorm(0.0, 0.0001);

      #prior for mean log Odds ratio
      tau ~ dunif(0, 100);

      #Prior for tau
      precision <- 1/(tau*tau);
      tau.sq <- 1/precision;

      #Calculate probability less 0
      less0 <- min(delta.new, 0);
      less0mean<-min(mu, 0);
      prob0 <- 1 - equals(less0, 0);
      prob0mean<- 1 - equals(less0mean, 0);
    }
", file="modele_magnesium.jags")

# n.chains<-3
# n.iter<-1500000
# n.burn<-1500000
# n.thin<- 1500

n.chains<-1
n.iter<-3000
n.burn<-3000
n.thin<- 100

modele <- jags.model(file="modele_magnesium.jags", data=don, 
                     #inits=list(X=matrix(1,ncol=Ntemps,nrow=Nind),beta1=), 
                     n.chains=n.chains, quiet=FALSE)
update(modele, n.iter=n.burn)
RES <- coda.samples(model=modele, variable.names=c("rc","rm","pc","pm","delta.new","delta"), 
                    n.iter=n.iter,thin=n.thin)

#save(RES,file="RES.rda")
#load("RES.rda")

plot(exp(RES[,"delta[1]"][[1]]))

exp(mean(RES[,"delta.new"][[1]]))

exp(mean(RES[,"delta[1]"][[1]]))
exp(mean(RES[,"delta[2]"][[1]]))
exp(mean(RES[,"delta[3]"][[1]]))
exp(mean(RES[,"delta[4]"][[1]]))
exp(mean(RES[,"delta[5]"][[1]]))
exp(mean(RES[,"delta[6]"][[1]]))
exp(mean(RES[,"delta[7]"][[1]]))
exp(mean(RES[,"delta[8]"][[1]]))
exp(mean(RES[,"delta[9]"][[1]]))
exp(mean(RES[,"delta[10]"][[1]]))
exp(mean(RES[,"delta[11]"][[1]]))
exp(mean(RES[,"delta[12]"][[1]]))
exp(mean(RES[,"delta[13]"][[1]]))
exp(mean(RES[,"delta[14]"][[1]]))