##########################################################
#                                                        #
# Functions for running Zero inflated Poisson Log-Normal #
#                                                        #
# By: Viktor Jonsson,   2016-10-15                       #
#                                                        #
##########################################################



ZiMGLMMLoan <- function(Y,group,N,nR=dim(Y)[1],nrAdapt=1000,
                             nrBurnIn=1000, nrSamples=10000, nrThin = 1,
                             modelName = "ZiGLMM.jags",
                             returnBeta1=FALSE,nRchains=1) {
  
  #Takes a matrix as input a long with group definitions and normalizing constants
  #runs MetaBayes and returns a vector of pvalues.
  
  
  nC <- sum(groups)
  
  groups <- c(sum(group==0),sum(group==1))
  
  
  
  
  
  
  nrColumns <- dim(Y)[2]
  nrRows <- dim(Y)[1]
  
  
  
  ###Obs out of habit the median of the total sums are used.
  normalizingConstants <- N/median(N)
  
  #DATA FOR MCMC OBJECT
  
  dataL <- list(Y=Y,N = normalizingConstants,group=groups)
  
  #INITS FOR MCMC OBJECT
  
  initsL1 <- ZiGLMMloaninits(Y=Y,n.chain = nRchains)
  
  
  ## GAMMA MODEL
  
  
  
  model1 <- jags.model(modelName,data = dataL,inits = initsL1, n.chains = 1,n.adapt = nrAdapt)
  bf <- proc.time()
  update(model1,n.iter=nrBurnIn)
  
  model1$rawData <- Y
  
  
  
  model1$samples <- coda.samples(model1,c("beta","p"),n.iter = nrSamples,thin=nrThin)
  
  
  
  
  
  #significance score calculation for ZiGLMM
  
  tempMatrix <- as.matrix(model1$samples[,1:nrRows])
  
  
  
  betaMeans <- apply(tempMatrix,2,mean)
  betaSds <- apply(tempMatrix,2,sd)
  
  betaProbUpper <- pnorm(0,betaMeans,betaSds)
  betaProbLower <- pnorm(0,betaMeans,betaSds,lower.tail=FALSE)
  
  betaProbMatrix <- cbind(betaProbUpper,betaProbLower)
  
  betaProb <- apply(betaProbMatrix,1,min)*2
  
  tempMatrix2 <- as.matrix(model1$samples[,(nrRows+1):(2*nrRows)])
  pMeans <- apply(tempMatrix2,2,mean)
  

  
  if(returnBeta1) {
    return(list(betaProb,betaMeans,betaSds,pMeans))
  } else {
    return(betaProb)
  }
  
  
}

ZiGLMMloaninits <- function(Y,n.chain = 1){
  nR = dim(Y)[1]
  Y <- as.matrix(Y)
  if (n.chain == 1) {
    
    alpha <- log(apply(Y,1,mean))
    beta <- runif(nR,-0.2,0.2)

    kappa <- 1
    eta <- 1
    binVar <- runif(nR,0.05,0.1)
    phi <- alpha + (log(Y+0.1)-alpha)
    p <- rep(0,nR)
    pi <- Y*0
    

    
    return(list(alpha = alpha,beta=beta,phi=phi,binVar = binVar,
                kappa=kappa, eta=eta,p=p,pi=pi))
  } else {
    initList <- list()
    for (i in 1:n.chain) {
      alpha <- log(apply(Y,1,mean))
      beta <- runif(nR,-0.2,0.2)
      
      
      kappa <- 1
      eta <- 1
      binVar <- runif(nR,0.05,0.1)
      phi <- alpha + (log(Y+0.1)-alpha)
      p <- rep(0,nR)
      pi <- Y*0
      

      
      initList[[i]] <- list(alpha = alpha,beta=beta,phi=phi,binVar = binVar,
                            kappa=kappa, eta=eta,p=p,pi=pi)
    }
    return(initList)
    
  }
  
}