R/describeISE.R

In ISEtools: Ion Selective Electrodes Analysis Methods

describeISE = function(data, model.path=NA, model.name = NA, Z=NA, temperature = 21, 
                       burnin=25000, iters = 50000, chains=4, thin = 1,
                       a.init= NA, b.init=NA, cstar.init=NA, logc.limits = c(-8.9, -1.9), sigma.upper = 5, diagnostic.print=FALSE, offset = 1, 
                       alpha = 0.05, beta = 0.05, SN = NA, 
                       keep.coda=TRUE, coda.n=1000, program="OpenBUGS") {
  
  ###################################################################
  # 				describeISE::                             #
  #   			estimate ISE parameters                         #
  ###################################################################
  
  ###
  # Calculates calibration curve to characterise an ISE
  # Required inputs are:
  #  data:        Data formatted using loadISEdata format
  #  model.path:  The directory where the BUGS model is located (e.g. "c:/myisefolder")
  #  model.name:  The name of the BUGS model (e.g. "Single ISE model Calibration only.txt")
  #
  # Technical inputs that do not require specification: 
  #  burnin:      The number of initial parameter estimates to discard while tuning the system
  #  iters:       The total number of iterations
  #  chains:      The number of parallel chains to run (useful for diagnostics)
  #  ### Note:    The final estimates are based on (iters - burnin)*chains number of simulations
  #  a.init:      Initial value passed to gen.inits.single (scalar) or gen.inits.multiple (vector)
  #  b.init:      Similar
  #  cstar.init:  Similar
  #  logc.limits: Upper and lower climits for log(c) initial values (defaults to -8.9 and -1.9)
  #  sigma.upper: Upper limit for sigma.init (use a smaller value if you receive an out of bounds error)
  #  diagnostic.print:  Prints information used to determine convergence (OpenBUGS only)
  #  offset:	The initial value for the slope uses the last and the last - offset point (defaults to 1).
  #  ### Note:    Default values may be altered if there are convergence problems in the 
  #               non-linear regression
  #  alpha:       False positive rate used for detection threshold (not output) to calculate LOD(alpha, beta) {only returned if SN = NA}
  #  beta:        False negative rate used to calculate LOD(alpha, beta) {only returned if SN = NA}
  #  SN:          Desired signal-to-noise ratio for LOD(S/N) calculations (default is to calculate the S/N equivalent based on alpha, beta)
  #  keep.coda    Logical flag indicating whether the MCMC simulations should be returned (keep.coda = TRUE) or not (keep.coda = FALSE)
  #  coda.n       Indicates how mnay simulations to return (sampled with replacement).  If coda.n > the total, all are returned
  #  program      Choice of "OpenBUGS" or "jags" for MCMC sampler, defaults to "OpenBUGS"
  ###
  
  # Error checking
  if (is.na(Z)) {
    stop("Please specify the ionic valence using Z=<value>, e.g. describeISE(..., Z = 2,...)", call.=FALSE)
  }
  
  # Check if required packages are installed and loaded
  if (!(program %in% c("OpenBUGS", "jags"))) {
    stop("Exiting: program must equal OpenBUGS or jags.", call.=FALSE)
  }
  
  if (program == "OpenBUGS") {
    OB.x = !identical(find.package("R2OpenBUGS", quiet=TRUE), character(0))
    if (!OB.x) {
      stop("The R2OpenBUGS package must be installed to use OpenBUGS.", call.=FALSE)
    }
  }
  
  if (program == "jags") {
    j.x = !identical(find.package("rjags", quiet=TRUE), character(0))
    if (!j.x) {
      stop("The rjags package must be installed to use jags.", call.=FALSE)
    }
  }
  
  # Flag for 1 or more than 1 ISEs
  single.ISE = FALSE
  if (data$R == 1) { single.ISE = TRUE }
  
  # If model.path is not specified, revert to appropriate defauls
  if (is.na(model.path)) { model.path = paste(path.package('ISEtools'), "/models", sep="") }
  
  # If model.name is not specified, revert to appropriate defaults
  if (is.na(model.name)) {
    if (single.ISE) {
      model.name = "Single_ISE_calibration_model.txt"
    }
    if (!single.ISE) {
      model.name = "Multiple_ISE_calibration_model.txt"
    }
  }
  
  # Location of the BUGs model
  ISE.model <- file.path(model.path, model.name)
  
  # Reduce the dataset to those variables expected by BUGS
  alldata = data
  data = list()
  data$N = alldata$N; data$log10x = alldata$log10x; data$emf = alldata$emf
  if (!single.ISE) { data$R = alldata$R; data$ISEID = alldata$ISEID }
  
  # Initial values
  message("Generating initial values...")
  if (single.ISE==TRUE) {
    # Number of ISEs = 1
    Num.ISEs = 1
    # Generate initial values
    inits <- rep(list(NA), chains)
    for (i in 1:chains) {
      inits[[i]] <- gen.inits.single(data, a.init=a.init, b.init=b.init,
                                     cstar.init= cstar.init, logc.limits=logc.limits, sigma.upper=sigma.upper, stdadd = F, offset=offset, calibration.only=T)
    }
  }
  if (single.ISE==FALSE) {
    # Generate initial values
    inits <- rep(list(NA), chains)
    # Number of ISEs given in data
    Num.ISEs = data$R
    for (i in 1:chains) {
      inits[[i]] <- gen.inits.multiple(data, a.init=a.init, b.init=b.init, 
                                       cstar.init= cstar.init, logc.limits=logc.limits, sigma.upper= sigma.upper, stdadd = F, offset=offset, calibration.only=T)
    }
  }
  # Parameters to monitor
  parameters <- c("a", "b", "c", "cstar", "sigma")
  
  # Calculate theoretical slope at mV level
  data$mu.b = 1000*log(10)*8.31433*(temperature + 273.15)/(96487*Z)
  
  if(program=="OpenBUGS") {
    # Call BUGS
    message("Running Bayesian model using OpenBUGS...")
    # Using R2OpenBUGS
    ISE.Bayes <- R2OpenBUGS::bugs(data, inits, parameters, n.iter = iters, 
			model.file=ISE.model, n.chains = chains, n.burnin=burnin, 
			n.thin=thin, digits=9)
 
       
    # Print and plot diagnostics
    if(diagnostic.print==TRUE) {
      print(ISE.Bayes, digits.summary=9)
      plot(ISE.Bayes)
    }
  }
  
  if(program=="jags") {
    # Call jags
    message("Running Bayesian model using jags...")
    ISE.Bayes.mod = rjags::jags.model(file = ISE.model, data=data, n.chains=chains, inits=inits, n.adapt=1000, quiet=TRUE)
    update(ISE.Bayes.mod, n.iter = burnin, progress.bar="none")
    my.coda = rjags::coda.samples(ISE.Bayes.mod, parameters, n.iter = (iters - burnin), thin = thin)
    
    # Align with OpenBUGS output
    ISE.Bayes = list()
    
    ISE.Bayes$sims.array = array(NA, c(chains, (iters - burnin)/thin, Num.ISEs*5) )
    for (i in 1:chains) {
      ISE.Bayes$sims.array[i,,] = my.coda[[i]]
    }
  }
  
  if(single.ISE==TRUE) {
    # Other vars
    ahat.coda = as.vector(ISE.Bayes$sims.array[,,1])
    coda.length = length(ahat.coda)
    bhat.coda = as.vector(ISE.Bayes$sims.array[,,2])
    cstarhat.coda = as.vector(ISE.Bayes$sims.array[,,4])
    chat.coda = cstarhat.coda^10
    sigmahat.coda = as.vector(ISE.Bayes$sims.array[,,5])
    ahat = median(ahat.coda)
    bhat = median(bhat.coda)
    cstarhat = median(cstarhat.coda)
    chat = cstarhat^10
    sigmahat = median(sigmahat.coda)
    
    ahat.lcl = quantile(ahat.coda, 0.025)
    bhat.lcl = quantile(bhat.coda, 0.025)
    cstarhat.lcl = quantile(cstarhat.coda, 0.025)
    chat.lcl = cstarhat.lcl^10
    ahat.ucl = quantile(ahat.coda, 0.975)
    bhat.ucl = quantile(bhat.coda, 0.975)
    cstarhat.ucl = quantile(cstarhat.coda, 0.975)
    chat.ucl = cstarhat.ucl^10
    sigmahat.lcl = quantile(sigmahat.coda, 0.025)
    sigmahat.ucl = quantile(sigmahat.coda, 0.975)
    
    if(is.na(SN)){
      if(is.na(alpha)) alpha = 0.05 # Reset to default if needed
      if(is.na(beta)) beta = 0.05
      zscore = qnorm(alpha, lower.tail=FALSE) + qnorm(beta, lower.tail=FALSE)
      LOD.info = list(type = "alpha, beta", alpha = alpha, beta = beta, SN = zscore)
    }else{
      zscore = SN
      LOD.info = list(type = "S/N", alpha = NA, beta = NA, SN = zscore)
    }
    LOD.coda = 10*log10(cstarhat.coda) + log10(10^(zscore*sigmahat.coda /abs( bhat.coda ) ) - 1)
    LOD.hat = median(LOD.coda)
    LOD.Q1 = quantile(LOD.coda, 0.25)
    LOD.Q3 = quantile(LOD.coda, 0.75)
    LOD.lcl = quantile(LOD.coda, 0.025)
    LOD.ucl = quantile(LOD.coda, 0.975)
  }
  
  if (single.ISE == FALSE) {	
    # Other vars
    coda.length = length(as.vector(ISE.Bayes$sims.array[,,1]))
    ahat.coda = matrix(NA, nrow = coda.length, ncol = data$R)
    bhat.coda <- cstarhat.coda <- chat.coda <- sigmahat.coda <- LOD.coda <- ahat.coda
    
    ahat = rep(NA, data$R)
    ahat.lcl = rep(NA, data$R)
    ahat.ucl = rep(NA, data$R)
    bhat = rep(NA, data$R)
    bhat.lcl = rep(NA, data$R)
    bhat.ucl = rep(NA, data$R)
    chat = rep(NA, data$R)
    chat.lcl = rep(NA, data$R)
    chat.ucl = rep(NA, data$R)
    cstarhat = rep(NA, data$R)
    cstarhat.lcl = rep(NA, data$R)
    cstarhat.ucl = rep(NA, data$R)
    sigmahat = rep(NA, data$R)
    sigmahat.lcl = rep(NA, data$R)
    sigmahat.ucl = rep(NA, data$R)
    
    LOD.hat = rep(NA, data$R)
    LOD.lcl = rep(NA, data$R)
    LOD.ucl = rep(NA, data$R)
    LOD.Q1 = rep(NA, data$R)
    LOD.Q3 = rep(NA, data$R)
    if(is.na(SN)){
      if(is.na(alpha)) alpha = 0.05 # Reset to default if needed
      if(is.na(beta)) beta = 0.05
      zscore = qnorm(alpha, lower.tail=FALSE) + qnorm(beta, lower.tail=FALSE)
      LOD.info = list(type = "alpha, beta", alpha = alpha, beta = beta, SN = zscore)
    }else{
      zscore = SN
      LOD.info = list(type = "S/N", alpha = NA, beta = NA, SN = zscore)
    }
    
    for (j in 1:data$R) {
      ahat.coda[,j] = as.vector(ISE.Bayes$sims.array[,,j])
      bhat.coda[,j] = as.vector(ISE.Bayes$sims.array[,,data$R + j])
      cstarhat.coda[,j] = as.vector(ISE.Bayes$sims.array[,,3*data$R + j])
      sigmahat.coda[,j] = as.vector(ISE.Bayes$sims.array[,,4*data$R + j])
      chat.coda[,j] = cstarhat.coda[,j]^10
      
      ahat.q = quantile(ahat.coda[,j], c(0.025, 0.5, 0.975))
      ahat[j] = ahat.q[2]
      ahat.lcl[j] = ahat.q[1]
      ahat.ucl[j] = ahat.q[3]
      bhat.q = quantile(bhat.coda[,j], c(0.025, 0.5, 0.975))
      bhat[j] = bhat.q[2]
      bhat.lcl[j] = bhat.q[1]
      bhat.ucl[j] = bhat.q[3]
      chat.q = quantile(chat.coda[,j], c(0.025, 0.5, 0.975))
      chat[j] = chat.q[2]
      chat.lcl[j] = chat.q[1]
      chat.ucl[j] = chat.q[3]
      cstarhat.q = quantile(cstarhat.coda[,j], c(0.025, 0.5, 0.975))
      cstarhat[j] = cstarhat.q[2]
      cstarhat.lcl[j] = cstarhat.q[1]
      cstarhat.ucl[j] = cstarhat.q[3]
      sigmahat.q = quantile(sigmahat.coda[,j], c(0.025, 0.5, 0.975))
      sigmahat[j] = sigmahat.q[2]
      sigmahat.lcl[j] = sigmahat.q[1]
      sigmahat.ucl[j] = sigmahat.q[3]
      
      LOD.coda[,j] = 10*log10(cstarhat.coda[,j]) + log10(10^(zscore*sigmahat.coda[,j] /abs( bhat.coda[,j] ) ) - 1)
      LOD.q = quantile(LOD.coda[,j], c(0.025, 0.25, 0.5, 0.75, 0.975))
      LOD.hat[j] = LOD.q[3]
      LOD.lcl[j] = LOD.q[1]
      LOD.ucl[j] = LOD.q[5]
      LOD.Q1[j] = LOD.q[2]
      LOD.Q3[j] = LOD.q[3]
    }	
  }
  
  ###
  # Relevant output
  ###
  if (keep.coda) {
    if (!is.na(coda.n)) {
      if (coda.n > coda.length) {
        message("Warning: requested coda sub-sample is longer than coda. Full coda returned.")
        coda.n = coda.length
      }
      coda.sample = sample(seq(coda.length), coda.n, replace=TRUE) 
      if (single.ISE == TRUE) {
        ahat.coda = ahat.coda[coda.sample]
        bhat.coda = bhat.coda[coda.sample]
        chat.coda = chat.coda[coda.sample]
        cstarhat.coda = cstarhat.coda[coda.sample]
        sigmahat.coda = sigmahat.coda[coda.sample]
        LOD.coda = LOD.coda[coda.sample]	
      }
      if (single.ISE == FALSE) {
        ahat.coda = ahat.coda[coda.sample,]
        bhat.coda = bhat.coda[coda.sample,]
        chat.coda = chat.coda[coda.sample,]
        cstarhat.coda = cstarhat.coda[coda.sample,]
        sigmahat.coda = sigmahat.coda[coda.sample,]
        LOD.coda = LOD.coda[coda.sample,]
      }
    }
    ISEdescription = list(ahat.coda=ahat.coda, bhat.coda=bhat.coda, chat.coda = chat.coda, sigmahat.coda=sigmahat.coda, cstarhat.coda=cstarhat.coda, 
                          LOD.coda=LOD.coda, 
                          ahat=ahat, bhat=bhat, chat=chat, cstarhat=cstarhat, sigmahat=sigmahat, 
                          ahat.lcl=ahat.lcl, ahat.ucl=ahat.ucl, bhat.lcl=bhat.lcl, bhat.ucl=bhat.ucl, 
                          chat.lcl=chat.lcl, chat.ucl=chat.ucl, cstarhat.lcl=cstarhat.lcl, cstarhat.ucl=cstarhat.ucl, sigmahat.lcl=sigmahat.lcl, sigmahat.ucl=sigmahat.ucl,
                          LOD.info=LOD.info, LOD.hat = LOD.hat, LOD.lcl = LOD.lcl, LOD.ucl = LOD.ucl,
                          LOD.Q1 = LOD.Q1, LOD.Q3 = LOD.Q3)
  }
  if (!keep.coda) {
    ISEdescription = list(
      ahat=ahat, bhat=bhat, chat=chat, cstarhat=cstarhat, sigmahat=sigmahat, 
      ahat.lcl=ahat.lcl, ahat.ucl=ahat.ucl, bhat.lcl=bhat.lcl, bhat.ucl=bhat.ucl, 
      chat.lcl=chat.lcl, chat.ucl=chat.ucl, cstarhat.lcl=cstarhat.lcl, cstarhat.ucl=cstarhat.ucl, sigmahat.lcl=sigmahat.lcl, sigmahat.ucl=sigmahat.ucl,
      LOD.info=LOD.info, LOD.hat = LOD.hat, LOD.lcl = LOD.lcl, LOD.ucl = LOD.ucl,
      LOD.Q1 = LOD.Q1, LOD.Q3 = LOD.Q3)	
  }
  class(ISEdescription) = "ISEdescription"
  return(ISEdescription)
}