Nothing
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)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.