Nothing
#Copyright 2022 NIEHS <matt.wheeler@nih.gov>
#
#
#Permission is hereby granted, free of charge, to any person obtaining a copy of this software
#and associated documentation files (the "Software"), to deal in the Software without restriction,
#including without limitation the rights to use, copy, modify, merge, publish, distribute,
#sublicense, and/or sell copies of the Software, and to permit persons to whom the Software
#is furnished to do so, subject to the following conditions:
#
#The above copyright notice and this permission notice shall be included in all copies
#or substantial portions of the Software.
#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
#INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
#PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
#HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
#CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
#OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
.evaluate_alpha <- function(...){
args <- list(...)
cl <- match.call()
ev <- match("alpha",names(cl))
if (is.na(ev)){
alpha = 0.05
}else{
alpha = args$alpha
if ( (alpha <= 0) || alpha > 0.5){
stop("alpha must be in (0,0.5]")
}
}
return(alpha)
}
.summary_continuous_max<-function(object,...){
model = object
returnV <- list()
alpha = .evaluate_alpha(...)
if (is.null(model$prior)){
returnV$fit_method <- "MLE"
returnV$prior <- NA
}else{
returnV$fit_method <- "Bayesian:MAP"
returnV$prior <- model$prior
}
returnV$fit <- model$full_model
temp_function <- tryCatch(splinefun(model$bmd_dist[,2],model$bmd_dist[,1],method="monoH.FC"),error = function(e) {NULL})
if (is.null(temp_function)){
returnV$BMD <- c(NA,NA,NA)
}else{
returnV$BMD <- temp_function(1-c(1-alpha,0.5,alpha))
}
names(returnV$BMD) <- c("BMDL","BMD","BMDU")
returnV$alpha <- alpha
temp <- model$Deviance - matrix(rep(model$Deviance[5,],5),nrow=5,ncol=2,byrow=T)
if (grepl("NCV",model$full_model)){
temp = rbind(temp[3,],temp[2,])
}else{
temp = rbind(temp[1,],temp[2,])
}
temp[,1] <- -1*temp[,1]
returnV$GOF <- cbind(temp,pchisq(temp[,1],temp[,2],lower.tail = F))
colnames(returnV$GOF) <- c("-2LL","DF","P-Value")
class(returnV) <- "summary_continuous_max"
return(returnV)
}
.print_summary_continuous_max<-function(x, ...){ # nolint
s_fit <- x
if (grepl("MLE",s_fit$fit_method)){
cat(sprintf("Summary of single model fit (%s) using ToxicR\n","MLE"))
cat(s_fit$fit,"\n")
}else{
cat(sprintf("Summary of single model fit (%s) using ToxicR\n\n","Bayesian-MAP"))
s_fit$GOF[,2] = round(s_fit$GOF[,2],2)
}
cat("\n")
cat("BMD: ")
cat(sprintf("%1.2f (%1.2f, %1.2f) %1.1f%% CI\n",s_fit$BMD[2],s_fit$BMD[1],s_fit$BMD[3],100*(1-2*s_fit$alpha)))
cat("\n")
cat("Model GOF\n")
cat("--------------------------------------------------\n")
s_fit$GOF[,1] <- round(s_fit$GOF[,1],2)
s_fit$GOF[,3] <- round(s_fit$GOF[,3],3)
rownames(s_fit$GOF) <- c("Test: Mean Adequate","Test: Mean/Variance Adequate")
print(s_fit$GOF)
}
.summary_continuous_mcmc<-function(object,...){
model = object
returnV <- list()
alpha = .evaluate_alpha(...)
returnV$fit_method <- "Bayesian:MCMC"
returnV$prior <- model$prior
returnV$fit <- model$full_model
temp_function <- tryCatch(splinefun(model$bmd_dist[,2],model$bmd_dist[,1],method="monoH.FC"),error = function(e) {NULL})
if (is.null(temp_function)){
returnV$BMD <- c(NA,NA,NA)
}else{
returnV$BMD <- temp_function(1-c(1-alpha,0.5,alpha))
}
names(returnV$BMD) <- c("BMDL","BMD","BMDU")
returnV$alpha <- alpha
returnV$eff_size <- coda::effectiveSize(model$mcmc_result$BMD_samples)
returnV$geweke_z <- coda::geweke.diag(coda::as.mcmc(model$mcmc_result$BMD_samples),frac1 = 0.3, frac2 = 0.4)$z
class(returnV) <- "summary_mcmc"
return(returnV)
}
.print_summary_continuous_mcmc<-function(x, ...){
s_fit = x
cat(sprintf("Summary of single model fit (%s) using ToxicR\n","MCMC"))
cat(s_fit$fit,"\n")
cat("\n")
cat("BMD: ")
cat(sprintf("%1.2f (%1.2f, %1.2f) %1.1f%% CI\n",s_fit$BMD[2],s_fit$BMD[1],s_fit$BMD[3],100*(1-2*s_fit$alpha)))
cat("\n")
cat("Convergence Diagnostics on BMD\n")
cat("--------------------------------------------------\n")
cat(sprintf("Effective Sample Size: %1.2f\n\n", s_fit$eff_size) )
cat(sprintf("Geweke Z-score that mean of first 30%% of \nMCMC chain is different from last 40%%\nZ-Score: %1.3f P-value %1.3f\n",
s_fit$geweke_z,2*pnorm(abs(s_fit$geweke_z),lower.tail=F)))
}
.summary_ma_max<-function(object, ...){
model = object
alpha = .evaluate_alpha(...)
returnV <- list()
returnV$fit_method <- "Bayesian:MAX"
returnV$fit_table <- data.frame(post_p = round(model$posterior_probs,3))
tmp_idx <- grep("Individual",names(model))
temp_mfit <- rep(" ",length(tmp_idx)) #model name
temp_BMD <- rep(" ",length(tmp_idx)) #bmd
temp_BMDL <- rep(" ",length(tmp_idx)) #bmdl
temp_BMDU <- rep(" ",length(tmp_idx)) #bmdu
for (ii in tmp_idx){
tmp_fit <- model[[ii]]
data_temp = tmp_fit$bmd_dist
dist = data_temp[!is.infinite(data_temp[,1]) & !is.na(data_temp[,1]),]
dist = data_temp[!is.nan(data_temp[,1])]
if (length(dist)>10 & !identical(dist, numeric(0))){
temp_function <- splinefun(data_temp[,2],data_temp[,1],method="monoH.FC")
temp_bmds <- temp_function(1-c(1-alpha,0.5,alpha))
temp_mfit[ii] <- sub("Model: ","",tmp_fit$full_model)
temp_BMD[ii] <- round(temp_bmds[2],3)
temp_BMDL[ii] <- round(temp_bmds[1],3)
temp_BMDU[ii] <- round(temp_bmds[3],3)
}else{
temp_mfit[ii] <- sub("Model: ","",tmp_fit$full_model)
temp_BMD[ii] <- NA
temp_BMDL[ii] <- NA
temp_BMDU[ii] <- NA
}
}
returnV$fit_table$model_names = temp_mfit
returnV$fit_table$BMD = temp_BMD
returnV$fit_table$BMDL = temp_BMDL
returnV$fit_table$BMDU = temp_BMDU
tmp_idx = order(returnV$fit_table$post_p,decreasing=T)
returnV$fit_table = returnV$fit_table[tmp_idx,c(2,3,4,5,1)]
warnFunc <- function(w){
return()
}
tryCatch({temp_function <- splinefun(model$ma_bmd[,2],model$ma_bmd[,1],method="monoH.FC")},
warning = warnFunc)
returnV$BMD <- temp_function(1-c(1-alpha,0.5,alpha))
names(returnV$BMD) <- c("BMDL","BMD","BMDU")
returnV$alpha <- alpha
class(returnV) <- "ma_summary_max"
return(returnV)
}
.summary_ma_mcmc<-function(object, ...){
model = object
alpha = .evaluate_alpha(...)
returnV <- list()
returnV$fit_method <- "Bayesian:MCMC"
returnV$fit_table <- data.frame(post_p = round(model$posterior_probs,3))
tmp_idx <- grep("Individual",names(model))
temp_mfit <- rep(" ",length(tmp_idx)) #model name
temp_BMD <- rep(" ",length(tmp_idx)) #bmd
temp_BMDL <- rep(" ",length(tmp_idx)) #bmdl
temp_BMDU <- rep(" ",length(tmp_idx)) #bmdu
for (ii in tmp_idx){
tmp_fit <- model[[ii]]
data_temp = tmp_fit$bmd_dist
dist = data_temp[!is.infinite(data_temp[,1]) & !is.na(data_temp[,1]),]
dist = data_temp[!is.nan(data_temp[,1])]
if (length(dist)>10 & !identical(dist, numeric(0))){
temp_function <- splinefun(data_temp[,2],data_temp[,1],method="monoH.FC")
temp_bmds <- temp_function(1-c(1-alpha,0.5,alpha))
temp_mfit[ii] <- sub("Model: ","",tmp_fit$full_model)
temp_BMD[ii] <- round(temp_bmds[2],3)
temp_BMDL[ii] <- round(temp_bmds[1],3)
temp_BMDU[ii] <- round(temp_bmds[3],3)
}else{
temp_mfit[ii] <- sub("Model: ","",tmp_fit$full_model)
temp_BMD[ii] <- NA
temp_BMDL[ii] <- NA
temp_BMDU[ii] <- NA
}
}
returnV$fit_table$model_names = temp_mfit
returnV$fit_table$BMD = temp_BMD
returnV$fit_table$BMDL = temp_BMDL
returnV$fit_table$BMDU = temp_BMDU
tmp_idx = order(returnV$fit_table$post_p,decreasing=T)
returnV$fit_table = returnV$fit_table[tmp_idx,c(2,3,4,5,1)]
warnFunc <- function(w){
return()
}
tryCatch({temp_function <- splinefun(model$ma_bmd[,2],model$ma_bmd[,1],method="monoH.FC")},
warning = warnFunc)
returnV$BMD <- temp_function(1-c(1-alpha,0.5,alpha))
names(returnV$BMD) <- c("BMDL","BMD","BMDU")
returnV$alpha <- alpha
class(returnV) <- "ma_summary_mcmc"
return(returnV)
}
.print_summary_ma<-function(x, ...){ # nolint
s_fit <- x
cat("Summary of single MA BMD\n\n")
cat("Individual Model BMDS\n")
cat(paste("Model",strrep(' ',34),sep=''),"\t\t BMD (BMDL, BMDU)\tPr(M|Data)\n")
cat("___________________________________________________________________________________________\n")
for (ii in 1:nrow(s_fit$fit_table)){
tmp_length = nchar(s_fit$fit_table[ii,1])
pad <- paste(substr(s_fit$fit_table[ii,1],1,38),strrep(" ",39-tmp_length),sep="")
cat(sprintf("%s\t\t\t%1.2f (%1.2f ,%1.2f) \t %1.3f\n",pad,as.numeric(s_fit$fit_table[ii,2]),
as.numeric(s_fit$fit_table[ii,3]),as.numeric(s_fit$fit_table[ii,4]),as.numeric(s_fit$fit_table[ii,5])))
}
cat("___________________________________________________________________________________________\n")
cat("Model Average BMD: ")
cat(sprintf("%1.2f (%1.2f, %1.2f) %1.1f%% CI\n",s_fit$BMD[2],s_fit$BMD[1],s_fit$BMD[3],
100*(1-2*s_fit$alpha)))
}
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.