R/plotFit.R
In USEPA/CompTox-ToxCast-tcpl: ToxCast Data Analysis Pipeline
Defines functions
.plotFit
#-------------------------------------------------------------------------------
# plotFit: Create plot of the dose-response with associated models
#-------------------------------------------------------------------------------
#' @importFrom graphics par layout plot rect abline curve axis axTicks points
#' @importFrom graphics plot.window text arrows
.plotFit <- function(resp, logc, pars) {
###--------------------------- Draw Left Panel ----------------------------###
layout(mat = matrix(1:2, nrow = 1), widths = c(4, 5.5), heights = 3.5)
on.exit(layout(1))
opar <- par()[c("pty", "mar", "family")]
on.exit(par(opar), add = TRUE)
par(pty = "s",
mar = c(4, 4.5, 2, 3) + 0.1,
family = "mono")
## Round all numeric values in 'pars' to 99 digits
nind <- which(sapply(pars, is.numeric))
pars[nind] <- lapply(pars[nind], round, digits = 99)
ylab <- NULL
if (pars$resp_unit == "percent_activity") {
#y0 <- c(-50, 150)
ylab <- "Percent Activity"
}
if (pars$resp_unit == "log2_fold_induction") {
#y0 <- c(-1, 4)
ylab <- "Log2(Fold Induction)"
}
if (pars$resp_unit == "log10_fold_induction") {
#y0 <- c(-0.1, 2)
ylab <- "Log10(Fold Induction)"
}
if (is.null(ylab)) {
ylab <- pars$resp_unit
#y0 <- c(-50, 150)
}
#new labels for x axis
xlab <- NULL
if (is.na(pars$conc_unit)) {
xlab <- expression(bold(paste("Concentration (",mu,"M)")))
}else{
if (pars$conc_unit == "uM") {
xlab <- expression(bold(paste("Concentration (",mu,"M)")))
}
if (pars$conc_unit == "mg/l") {
xlab <- "Concentration (mg/l)"
}
if (pars$conc_unit == "CF") {
xlab <- "Concentration (CF)"
}
}
if (is.null(xlab)) {
xlab <- paste0("Concentration (",pars$conc_unit,")")
}
if(pars$bmad != 0){
y0 <- c(signif(-10*pars$bmad,2), signif(20*pars$bmad,2))
}else{
if ("coff" %in% names(pars)) {
y0 <- c(signif(-5*pars$coff,2), signif(10*pars$coff,2))
}else{
y0 <- c(-0.1, 2)
}
}
if (!("fitc" %in% names(pars))) pars$fitc = -10000L # Create placeholder fitc for lvl 4
fmax <- suppressWarnings(with(pars, 1.05*max(hill_tp, gnls_tp, na.rm = TRUE)))
if (is.infinite(fmax)) fmax <- NA_real_
view <- fmax/diff(range(resp))
hbrk <- pars$resp_max > y0[2]
lbrk <- pars$resp_min < y0[1]
brk <- with(pars, view < 0.5 & (hbrk | lbrk))
pad <- if (hbrk & lbrk) 0.1 else 0.2
if (!is.na(brk) & brk) {
yrng <- (fmax - y0[1])/(1 - hbrk*pad - lbrk*pad)
ylim <- c(y0[1] - pad*yrng*lbrk, fmax + pad*yrng*hbrk)
md <- resp < fmax & resp > ylim[1]
if (all(md)) {
brk <- FALSE
} else {
hi <- resp > fmax
if (!any(hi)) hbrk <- FALSE
lo <- resp < ylim[1]
if (!any(lo)) lbrk <- FALSE
}
if (!any(lbrk, hbrk)) brk <- FALSE
} else {
if (is.na(brk)) {
brk <- if(hbrk | lbrk) TRUE else FALSE
if (brk) {
yrng <- diff(y0)/(1 - hbrk*pad - lbrk*pad)
ylim <- c(y0[1] - pad*yrng*lbrk, y0[2] + pad*yrng*hbrk)
md <- resp < ylim[2] & resp > ylim[1]
if (all(md)) {
brk <- FALSE
} else {
hi <- resp > ylim[2]
if (!any(hi)) hbrk <- FALSE
lo <- resp < ylim[1]
if (!any(lo)) lbrk <- FALSE
}
} else {
ylim <- y0
md <- rep(TRUE, length(resp))
}
} else {
ylim <- with(pars, c(min(y0[1], 1.2*resp_min), max(y0[2], 1.2*resp_max)))
md <- rep(TRUE, length(resp))
}
}
p <- list(ylim = ylim,
xlim = range(logc),
cex.lab = 1.2,
cex.axis = 1.2,
font.lab = 2,
col = "black",
cex = 2,
xlab = xlab,
ylab = ylab,
main = "",
bty = "n",
xaxt = "n",
yaxt = "n",
type = "n")
do.call(what = plot, args = c(resp[md] ~ logc[md], p), quote = TRUE)
rect(xleft = par()$usr[1],
xright = par()$usr[2],
ybottom = -3 * pars$bmad,
ytop = 3 * pars$bmad,
border = NA,
col = "gray70",
density = 15,
angle = 45)
if ("coff" %in% names(pars)) abline(h = pars$coff, lwd = 1.5, col = "gray70")
if (is.null(pars$modl)) pars$modl <- "none"
if (is.na(pars$modl)) pars$modl <- "none"
if (!is.na(pars$cnst) & pars$cnst) {
abline(h = 0,
lwd = 4,
col = "darkorange",
lty = ifelse(pars$modl == "cnst", "solid", "dashed"))
if (pars$fitc == 100L) {
# loec analysis
abline(v = pars$modl_acc,
lwd = 4,
col = 'tomato3')
}
}
if (!is.na(pars$hill) & pars$hill & pars$fitc != 100L) {
hill.eq <- function(x) with(pars, hill_tp/(1 + 10^((hill_ga - x)*hill_gw)))
curve(hill.eq,
from = pars$logc_min,
to = pars$logc_max,
add = T,
n = 1e4,
lwd = 4,
col = "tomato3",
lty = ifelse(pars$modl == "hill", "solid", "dashed"))
abline(v = pars$hill_ga,
lwd = 2.5,
lty = ifelse(pars$modl == "hill", "solid", "dashed"),
col = "tomato3")
}
if (!is.na(pars$gnls) & pars$gnls & pars$fitc != 100L) {
gnls.eq <- function(x) {
with(pars, {
h1 <- (1/(1 + 10^((gnls_ga - x)*gnls_gw)))
h2 <- (1/(1 + 10^((x - gnls_la)*gnls_lw)))
gnls_tp*h1*h2
})
}
curve(gnls.eq,
from = pars$logc_min,
to = pars$logc_max,
add = T,
n = 1e4,
lwd = 4,
col = "dodgerblue2",
lty = ifelse(pars$modl == "gnls", "solid", "dashed"))
abline(v = pars$gnls_ga,
lwd = 2.5,
lty = ifelse(pars$modl == "gnls", "solid", "dashed"),
col = "dodgerblue2")
}
if (is.null(pars$toxboot)) pars$toxboot <- NA
#error bars plotting
if (!is.na(pars$toxboot) & pars$toxboot) {
y = pars$modl_tp/2 #50% resp
if(!is.na(pars$modl_ga_min) & !is.na(pars$modl_ga_max))
suppressWarnings(arrows(pars$modl_ga_min, y, pars$modl_ga_max, y, code=3, angle=90, length=0.1, lwd = 2))
}
axis(side = 1,
at = axTicks(side = 1),
labels = signif(10^axTicks(side = 1), digits = 1),
font = 1,
lwd = 2,
cex.axis = 1.2,
col = "gray35")
axis(side = 2,
at = axTicks(side = 2),
labels = axTicks(side = 2),
font = 1,
lwd = 2,
cex.axis = 1.2,
col = "gray35")
# points(x = pars$emax_conc,
# y = pars$emax,
# pch = 22,
# cex = 2,
# col = "gray35",
# lwd = 1,
# bg = "yellow2")
points(resp[md] ~ logc[md], cex = 1.5, lwd = 2.5, col = "gray30")
if (brk) {
if (hbrk) {
hrng <- unique(range(resp[hi]))
if (length(hrng) != 1) {
hlim <- with(pars, c(resp_max - diff(hrng)/pad, resp_max))
} else {
hlim <- with(pars, c(resp_max - (hrng - y0[2])/pad, resp_max))
}
par(new = TRUE)
plot.window(xlim = par()$usr[1:2], ylim = hlim)
points(resp[hi] ~ logc[hi], cex = 0.5, lwd = 2.5, col = "gray60")
axis(side = 4,
at = hrng,
labels = signif(hrng, 2),
font = 1,
lwd = 2,
cex.axis = 0.5,
col = "gray60")
}
if (lbrk) {
lrng <- unique(range(resp[lo]))
if (length(lrng) != 1) {
llim <- with(pars, c(resp_min, resp_min + diff(lrng)/pad))
} else {
llim <- with(pars, c(resp_min, resp_min + (y0[1] - lrng)/pad))
}
par(new = TRUE)
plot.window(xlim = par()$usr[1:2], ylim = llim)
points(resp[lo] ~ logc[lo], cex = 0.5, lwd = 2.5, col = "gray60")
axis(side = 4,
at = lrng,
labels = signif(lrng, 2),
font = 1,
lwd = 2,
cex.axis = 0.5,
col = "gray60")
}
}
###--------------------- Prepare Text for Right Panel ---------------------###
spaces <- function(x) paste(rep(" ", x), collapse = "")
itxt <- with(pars, {
paste0("ASSAY: ", aenm, "\n\n",
"NAME: ", chnm, "\n",
"CHID: ", chid, spaces(8 - ifelse(is.na(chid), 2, nchar(chid))),
"CASRN: ", casn, "\n",
"SPID(S): ", spid, "\n",
"M4ID: ", m4id, " ", ifelse(brk, "BRK", ""), "\n\n"
)
})
if (!is.na(pars$hill) & pars$hill) {
if (pars$hcov) {
hsds <- with(pars, signif(c(hill_tp_sd, hill_ga_sd, hill_gw_sd), 3))
hsds[is.na(hsds)] <- "NaN"
} else {
hsds <- rep("NA", 3)
}
hprs <- with(pars, signif(c(hill_tp, hill_ga, hill_gw), 3))
htxt1 <- paste("HILL MODEL (in red):\n tp",
"ga",
"gw\n",
sep = spaces(7))
htxt2 <- paste0(c("val: ", "sd: "),
c(paste(sapply(hprs,
function(x) {
paste0(x, spaces(9 - nchar(x)))
}),
collapse = ""),
paste(sapply(hsds,
function(x) {
paste0(x, spaces(9 - nchar(x)))
}),
collapse = "")),
collapse = "\n")
htxt <- paste0(htxt1, htxt2, "\n\n")
} else {
if (is.na(pars$hill)) {
htxt <- "HILL MODEL: Not applicable.\n\n"
} else {
htxt <- "HILL MODEL: Failed to converge.\n\n"
}
}
if (!is.na(pars$gnls) & pars$gnls) {
if (pars$gcov) {
gsds <- with(pars,
signif(c(gnls_tp_sd,
gnls_ga_sd,
gnls_gw_sd,
gnls_la_sd,
gnls_lw_sd),
3)
)
gsds[is.na(gsds)] <- "NaN"
} else {
gsds <- rep("NA", 5)
}
gprs <- with(pars,
signif(c(gnls_tp, gnls_ga, gnls_gw, gnls_la, gnls_lw), 3))
gtxt1 <- paste("GAIN-LOSS MODEL (in blue):\n tp",
"ga",
"gw",
"la",
"lw\n",
sep = spaces(7))
gtxt2 <- paste0(c("val: ", "sd: "),
c(paste(sapply(gprs,
function(x) {
paste0(x, spaces(9 - nchar(x)))
}),
collapse = ""),
paste(sapply(gsds,
function(x) {
paste0(x, spaces(9 - nchar(x)))
}),
collapse = "")),
collapse = "\n")
gtxt <- paste0(gtxt1, gtxt2, "\n\n")
} else {
if (is.na(pars$hill)) {
gtxt <- "GAIN-LOSS MODEL: Not applicable.\n\n"
} else {
gtxt <- "GAIN-LOSS MODEL: Failed to converge.\n\n"
}
}
cvals <- function(x) {
x <- as.character(x)
x[is.na(x)] <- "NA"
x
}
aics <- cvals(with(pars, round(c(cnst_aic, hill_aic, gnls_aic), 2)))
prob <- cvals(with(pars, round(c(cnst_prob, hill_prob, gnls_prob), 2)))
rmse <- cvals(with(pars, round(c(cnst_rmse, hill_rmse, gnls_rmse), 2)))
models <- c("CNST", "HILL", "GNLS")
atxt <- paste0(spaces(6),
models[1],
spaces(8),
models[2],
spaces(8),
models[3],
"\n",
paste0("AIC: ",
aics[1],
spaces(12 - nchar(aics[1])),
aics[2],
spaces(12 - nchar(aics[2])),
aics[3]),
"\n",
paste0("PROB: ",
prob[1],
spaces(12 - nchar(prob[1])),
prob[2],
spaces(12 - nchar(prob[2])),
prob[3]),
"\n",
paste0("RMSE: ",
rmse[1],
spaces(12 - nchar(rmse[1])),
rmse[2],
spaces(12 - nchar(rmse[2])),
rmse[3]),
"\n\n")
pars$max_mean <- signif(pars$max_mean, 3)
pars$max_med <- signif(pars$max_med, 3)
ntxt <- paste0("MAX_MEAN: ", pars$max_mean,
spaces(10 - nchar(pars$max_mean)),
"MAX_MED: ", pars$max_med,
spaces(10 - nchar(pars$max_med)),
"BMAD: ", signif(pars$bmad, 3),
"\n\n")
if (!is.null(pars$hitc)) {
pars$coff <- signif(pars$coff, 3)
ctxt <- paste0("COFF: ", pars$coff, spaces(7 - nchar(pars$coff)),
"HIT-CALL: ", pars$hitc, spaces(5 - nchar(pars$hitc)),
"FITC: ", pars$fitc, spaces(5 - nchar(pars$fitc)),
"ACTP: ", round(pars$actp, 2),
"\n\n")
} else {
ctxt <- NULL
}
if (!is.null(pars$flgo)) {
ftxt <- paste0("FLAGS: ", ifelse(is.na(pars$flgo), "", pars$flgo),"\n\n")
} else {
ftxt <- NULL
}
if (!is.na(pars$toxboot)) {
btxt <- paste0("HIT-PCT: ", pars$hit_pct, spaces(2),
"MED-GA: ", round(pars$modl_ga_med,4), spaces(2),
"GA-CI: ", round(pars$modl_ga_delta,4), spaces(2),
"\n\n")
} else {
btxt <- NULL
}
plot_txt1 <- paste0(itxt, htxt, gtxt, atxt, ntxt, ctxt, ftxt, btxt)
if (pars$modl != "none") {
nlines <- sum(7,
length(gregexpr("\n", htxt)[[1]]),
length(gregexpr("\n", gtxt)[[1]]))
winner <- with(pars, which(c("cnst", "hill", "gnls") == pars$modl))
if (length(winner) > 1) winner <- winner[1]
plot_txt2 <- paste0(paste(rep("\n", nlines), collapse = ""),
spaces(6 + 12*(winner - 1)),
models[winner],
"\n",
spaces(6 + 12*(winner - 1)),
aics[winner],
"\n",
spaces(6 + 12*(winner - 1)),
prob[winner],
"\n",
spaces(6 + 12*(winner - 1)),
rmse[winner])
} else {
plot_txt2 <- NULL
}
###--------------------------- Draw Right Panel ---------------------------###
par(pty = "m",
family = "mono",
mar = rep(2,4) + 0.1)
plot(0,
type = "n",
bty = "n",
xaxt = "n",
yaxt = "n",
ylab = "",
xlab = "",
xlim = c(0, 16),
ylim = c(0, 16))
suppressWarnings(
text(y = 15,
x = 1,
labels = plot_txt1,
adj = c(0, 1),
font = 2,
cex = 1)
)
suppressWarnings(
text(y = 15,
x = 1,
labels = plot_txt2,
adj = c(0, 1),
font = 2,
cex = 1,
col = "red")
)
}
#-------------------------------------------------------------------------------
USEPA/CompTox-ToxCast-tcpl documentation built on April 22, 2024, 12:37 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions .plotFit
#-------------------------------------------------------------------------------
# plotFit: Create plot of the dose-response with associated models
#-------------------------------------------------------------------------------
#' @importFrom graphics par layout plot rect abline curve axis axTicks points
#' @importFrom graphics plot.window text arrows
.plotFit <- function(resp, logc, pars) {
###--------------------------- Draw Left Panel ----------------------------###
layout(mat = matrix(1:2, nrow = 1), widths = c(4, 5.5), heights = 3.5)
on.exit(layout(1))
opar <- par()[c("pty", "mar", "family")]
on.exit(par(opar), add = TRUE)
par(pty = "s",
mar = c(4, 4.5, 2, 3) + 0.1,
family = "mono")
## Round all numeric values in 'pars' to 99 digits
nind <- which(sapply(pars, is.numeric))
pars[nind] <- lapply(pars[nind], round, digits = 99)
ylab <- NULL
if (pars$resp_unit == "percent_activity") {
#y0 <- c(-50, 150)
ylab <- "Percent Activity"
}
if (pars$resp_unit == "log2_fold_induction") {
#y0 <- c(-1, 4)
ylab <- "Log2(Fold Induction)"
}
if (pars$resp_unit == "log10_fold_induction") {
#y0 <- c(-0.1, 2)
ylab <- "Log10(Fold Induction)"
}
if (is.null(ylab)) {
ylab <- pars$resp_unit
#y0 <- c(-50, 150)
}
#new labels for x axis
xlab <- NULL
if (is.na(pars$conc_unit)) {
xlab <- expression(bold(paste("Concentration (",mu,"M)")))
}else{
if (pars$conc_unit == "uM") {
xlab <- expression(bold(paste("Concentration (",mu,"M)")))
}
if (pars$conc_unit == "mg/l") {
xlab <- "Concentration (mg/l)"
}
if (pars$conc_unit == "CF") {
xlab <- "Concentration (CF)"
}
}
if (is.null(xlab)) {
xlab <- paste0("Concentration (",pars$conc_unit,")")
}
if(pars$bmad != 0){
y0 <- c(signif(-10*pars$bmad,2), signif(20*pars$bmad,2))
}else{
if ("coff" %in% names(pars)) {
y0 <- c(signif(-5*pars$coff,2), signif(10*pars$coff,2))
}else{
y0 <- c(-0.1, 2)
}
}
if (!("fitc" %in% names(pars))) pars$fitc = -10000L # Create placeholder fitc for lvl 4
fmax <- suppressWarnings(with(pars, 1.05*max(hill_tp, gnls_tp, na.rm = TRUE)))
if (is.infinite(fmax)) fmax <- NA_real_
view <- fmax/diff(range(resp))
hbrk <- pars$resp_max > y0[2]
lbrk <- pars$resp_min < y0[1]
brk <- with(pars, view < 0.5 & (hbrk | lbrk))
pad <- if (hbrk & lbrk) 0.1 else 0.2
if (!is.na(brk) & brk) {
yrng <- (fmax - y0[1])/(1 - hbrk*pad - lbrk*pad)
ylim <- c(y0[1] - pad*yrng*lbrk, fmax + pad*yrng*hbrk)
md <- resp < fmax & resp > ylim[1]
if (all(md)) {
brk <- FALSE
} else {
hi <- resp > fmax
if (!any(hi)) hbrk <- FALSE
lo <- resp < ylim[1]
if (!any(lo)) lbrk <- FALSE
}
if (!any(lbrk, hbrk)) brk <- FALSE
} else {
if (is.na(brk)) {
brk <- if(hbrk | lbrk) TRUE else FALSE
if (brk) {
yrng <- diff(y0)/(1 - hbrk*pad - lbrk*pad)
ylim <- c(y0[1] - pad*yrng*lbrk, y0[2] + pad*yrng*hbrk)
md <- resp < ylim[2] & resp > ylim[1]
if (all(md)) {
brk <- FALSE
} else {
hi <- resp > ylim[2]
if (!any(hi)) hbrk <- FALSE
lo <- resp < ylim[1]
if (!any(lo)) lbrk <- FALSE
}
} else {
ylim <- y0
md <- rep(TRUE, length(resp))
}
} else {
ylim <- with(pars, c(min(y0[1], 1.2*resp_min), max(y0[2], 1.2*resp_max)))
md <- rep(TRUE, length(resp))
}
}
p <- list(ylim = ylim,
xlim = range(logc),
cex.lab = 1.2,
cex.axis = 1.2,
font.lab = 2,
col = "black",
cex = 2,
xlab = xlab,
ylab = ylab,
main = "",
bty = "n",
xaxt = "n",
yaxt = "n",
type = "n")
do.call(what = plot, args = c(resp[md] ~ logc[md], p), quote = TRUE)
rect(xleft = par()$usr[1],
xright = par()$usr[2],
ybottom = -3 * pars$bmad,
ytop = 3 * pars$bmad,
border = NA,
col = "gray70",
density = 15,
angle = 45)
if ("coff" %in% names(pars)) abline(h = pars$coff, lwd = 1.5, col = "gray70")
if (is.null(pars$modl)) pars$modl <- "none"
if (is.na(pars$modl)) pars$modl <- "none"
if (!is.na(pars$cnst) & pars$cnst) {
abline(h = 0,
lwd = 4,
col = "darkorange",
lty = ifelse(pars$modl == "cnst", "solid", "dashed"))
if (pars$fitc == 100L) {
# loec analysis
abline(v = pars$modl_acc,
lwd = 4,
col = 'tomato3')
}
}
if (!is.na(pars$hill) & pars$hill & pars$fitc != 100L) {
hill.eq <- function(x) with(pars, hill_tp/(1 + 10^((hill_ga - x)*hill_gw)))
curve(hill.eq,
from = pars$logc_min,
to = pars$logc_max,
add = T,
n = 1e4,
lwd = 4,
col = "tomato3",
lty = ifelse(pars$modl == "hill", "solid", "dashed"))
abline(v = pars$hill_ga,
lwd = 2.5,
lty = ifelse(pars$modl == "hill", "solid", "dashed"),
col = "tomato3")
}
if (!is.na(pars$gnls) & pars$gnls & pars$fitc != 100L) {
gnls.eq <- function(x) {
with(pars, {
h1 <- (1/(1 + 10^((gnls_ga - x)*gnls_gw)))
h2 <- (1/(1 + 10^((x - gnls_la)*gnls_lw)))
gnls_tp*h1*h2
})
}
curve(gnls.eq,
from = pars$logc_min,
to = pars$logc_max,
add = T,
n = 1e4,
lwd = 4,
col = "dodgerblue2",
lty = ifelse(pars$modl == "gnls", "solid", "dashed"))
abline(v = pars$gnls_ga,
lwd = 2.5,
lty = ifelse(pars$modl == "gnls", "solid", "dashed"),
col = "dodgerblue2")
}
if (is.null(pars$toxboot)) pars$toxboot <- NA
#error bars plotting
if (!is.na(pars$toxboot) & pars$toxboot) {
y = pars$modl_tp/2 #50% resp
if(!is.na(pars$modl_ga_min) & !is.na(pars$modl_ga_max))
suppressWarnings(arrows(pars$modl_ga_min, y, pars$modl_ga_max, y, code=3, angle=90, length=0.1, lwd = 2))
}
axis(side = 1,
at = axTicks(side = 1),
labels = signif(10^axTicks(side = 1), digits = 1),
font = 1,
lwd = 2,
cex.axis = 1.2,
col = "gray35")
axis(side = 2,
at = axTicks(side = 2),
labels = axTicks(side = 2),
font = 1,
lwd = 2,
cex.axis = 1.2,
col = "gray35")
# points(x = pars$emax_conc,
# y = pars$emax,
# pch = 22,
# cex = 2,
# col = "gray35",
# lwd = 1,
# bg = "yellow2")
points(resp[md] ~ logc[md], cex = 1.5, lwd = 2.5, col = "gray30")
if (brk) {
if (hbrk) {
hrng <- unique(range(resp[hi]))
if (length(hrng) != 1) {
hlim <- with(pars, c(resp_max - diff(hrng)/pad, resp_max))
} else {
hlim <- with(pars, c(resp_max - (hrng - y0[2])/pad, resp_max))
}
par(new = TRUE)
plot.window(xlim = par()$usr[1:2], ylim = hlim)
points(resp[hi] ~ logc[hi], cex = 0.5, lwd = 2.5, col = "gray60")
axis(side = 4,
at = hrng,
labels = signif(hrng, 2),
font = 1,
lwd = 2,
cex.axis = 0.5,
col = "gray60")
}
if (lbrk) {
lrng <- unique(range(resp[lo]))
if (length(lrng) != 1) {
llim <- with(pars, c(resp_min, resp_min + diff(lrng)/pad))
} else {
llim <- with(pars, c(resp_min, resp_min + (y0[1] - lrng)/pad))
}
par(new = TRUE)
plot.window(xlim = par()$usr[1:2], ylim = llim)
points(resp[lo] ~ logc[lo], cex = 0.5, lwd = 2.5, col = "gray60")
axis(side = 4,
at = lrng,
labels = signif(lrng, 2),
font = 1,
lwd = 2,
cex.axis = 0.5,
col = "gray60")
}
}
###--------------------- Prepare Text for Right Panel ---------------------###
spaces <- function(x) paste(rep(" ", x), collapse = "")
itxt <- with(pars, {
paste0("ASSAY: ", aenm, "\n\n",
"NAME: ", chnm, "\n",
"CHID: ", chid, spaces(8 - ifelse(is.na(chid), 2, nchar(chid))),
"CASRN: ", casn, "\n",
"SPID(S): ", spid, "\n",
"M4ID: ", m4id, " ", ifelse(brk, "BRK", ""), "\n\n"
)
})
if (!is.na(pars$hill) & pars$hill) {
if (pars$hcov) {
hsds <- with(pars, signif(c(hill_tp_sd, hill_ga_sd, hill_gw_sd), 3))
hsds[is.na(hsds)] <- "NaN"
} else {
hsds <- rep("NA", 3)
}
hprs <- with(pars, signif(c(hill_tp, hill_ga, hill_gw), 3))
htxt1 <- paste("HILL MODEL (in red):\n tp",
"ga",
"gw\n",
sep = spaces(7))
htxt2 <- paste0(c("val: ", "sd: "),
c(paste(sapply(hprs,
function(x) {
paste0(x, spaces(9 - nchar(x)))
}),
collapse = ""),
paste(sapply(hsds,
function(x) {
paste0(x, spaces(9 - nchar(x)))
}),
collapse = "")),
collapse = "\n")
htxt <- paste0(htxt1, htxt2, "\n\n")
} else {
if (is.na(pars$hill)) {
htxt <- "HILL MODEL: Not applicable.\n\n"
} else {
htxt <- "HILL MODEL: Failed to converge.\n\n"
}
}
if (!is.na(pars$gnls) & pars$gnls) {
if (pars$gcov) {
gsds <- with(pars,
signif(c(gnls_tp_sd,
gnls_ga_sd,
gnls_gw_sd,
gnls_la_sd,
gnls_lw_sd),
3)
)
gsds[is.na(gsds)] <- "NaN"
} else {
gsds <- rep("NA", 5)
}
gprs <- with(pars,
signif(c(gnls_tp, gnls_ga, gnls_gw, gnls_la, gnls_lw), 3))
gtxt1 <- paste("GAIN-LOSS MODEL (in blue):\n tp",
"ga",
"gw",
"la",
"lw\n",
sep = spaces(7))
gtxt2 <- paste0(c("val: ", "sd: "),
c(paste(sapply(gprs,
function(x) {
paste0(x, spaces(9 - nchar(x)))
}),
collapse = ""),
paste(sapply(gsds,
function(x) {
paste0(x, spaces(9 - nchar(x)))
}),
collapse = "")),
collapse = "\n")
gtxt <- paste0(gtxt1, gtxt2, "\n\n")
} else {
if (is.na(pars$hill)) {
gtxt <- "GAIN-LOSS MODEL: Not applicable.\n\n"
} else {
gtxt <- "GAIN-LOSS MODEL: Failed to converge.\n\n"
}
}
cvals <- function(x) {
x <- as.character(x)
x[is.na(x)] <- "NA"
x
}
aics <- cvals(with(pars, round(c(cnst_aic, hill_aic, gnls_aic), 2)))
prob <- cvals(with(pars, round(c(cnst_prob, hill_prob, gnls_prob), 2)))
rmse <- cvals(with(pars, round(c(cnst_rmse, hill_rmse, gnls_rmse), 2)))
models <- c("CNST", "HILL", "GNLS")
atxt <- paste0(spaces(6),
models[1],
spaces(8),
models[2],
spaces(8),
models[3],
"\n",
paste0("AIC: ",
aics[1],
spaces(12 - nchar(aics[1])),
aics[2],
spaces(12 - nchar(aics[2])),
aics[3]),
"\n",
paste0("PROB: ",
prob[1],
spaces(12 - nchar(prob[1])),
prob[2],
spaces(12 - nchar(prob[2])),
prob[3]),
"\n",
paste0("RMSE: ",
rmse[1],
spaces(12 - nchar(rmse[1])),
rmse[2],
spaces(12 - nchar(rmse[2])),
rmse[3]),
"\n\n")
pars$max_mean <- signif(pars$max_mean, 3)
pars$max_med <- signif(pars$max_med, 3)
ntxt <- paste0("MAX_MEAN: ", pars$max_mean,
spaces(10 - nchar(pars$max_mean)),
"MAX_MED: ", pars$max_med,
spaces(10 - nchar(pars$max_med)),
"BMAD: ", signif(pars$bmad, 3),
"\n\n")
if (!is.null(pars$hitc)) {
pars$coff <- signif(pars$coff, 3)
ctxt <- paste0("COFF: ", pars$coff, spaces(7 - nchar(pars$coff)),
"HIT-CALL: ", pars$hitc, spaces(5 - nchar(pars$hitc)),
"FITC: ", pars$fitc, spaces(5 - nchar(pars$fitc)),
"ACTP: ", round(pars$actp, 2),
"\n\n")
} else {
ctxt <- NULL
}
if (!is.null(pars$flgo)) {
ftxt <- paste0("FLAGS: ", ifelse(is.na(pars$flgo), "", pars$flgo),"\n\n")
} else {
ftxt <- NULL
}
if (!is.na(pars$toxboot)) {
btxt <- paste0("HIT-PCT: ", pars$hit_pct, spaces(2),
"MED-GA: ", round(pars$modl_ga_med,4), spaces(2),
"GA-CI: ", round(pars$modl_ga_delta,4), spaces(2),
"\n\n")
} else {
btxt <- NULL
}
plot_txt1 <- paste0(itxt, htxt, gtxt, atxt, ntxt, ctxt, ftxt, btxt)
if (pars$modl != "none") {
nlines <- sum(7,
length(gregexpr("\n", htxt)[[1]]),
length(gregexpr("\n", gtxt)[[1]]))
winner <- with(pars, which(c("cnst", "hill", "gnls") == pars$modl))
if (length(winner) > 1) winner <- winner[1]
plot_txt2 <- paste0(paste(rep("\n", nlines), collapse = ""),
spaces(6 + 12*(winner - 1)),
models[winner],
"\n",
spaces(6 + 12*(winner - 1)),
aics[winner],
"\n",
spaces(6 + 12*(winner - 1)),
prob[winner],
"\n",
spaces(6 + 12*(winner - 1)),
rmse[winner])
} else {
plot_txt2 <- NULL
}
###--------------------------- Draw Right Panel ---------------------------###
par(pty = "m",
family = "mono",
mar = rep(2,4) + 0.1)
plot(0,
type = "n",
bty = "n",
xaxt = "n",
yaxt = "n",
ylab = "",
xlab = "",
xlim = c(0, 16),
ylim = c(0, 16))
suppressWarnings(
text(y = 15,
x = 1,
labels = plot_txt1,
adj = c(0, 1),
font = 2,
cex = 1)
)
suppressWarnings(
text(y = 15,
x = 1,
labels = plot_txt2,
adj = c(0, 1),
font = 2,
cex = 1,
col = "red")
)
}
#-------------------------------------------------------------------------------
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.