#####################################################################
## This program is distributed in the hope that it will be useful, ##
## but WITHOUT ANY WARRANTY; without even the implied warranty of ##
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ##
## GNU General Public License for more details. ##
#####################################################################
#-------------------------------------------------------------------------------
# plotFit: Create plot of the dose-response with associated models
#-------------------------------------------------------------------------------
#' @keywords internal
#'
#' @importFrom graphics par layout plot rect abline curve axis axTicks points
#' @importFrom graphics plot.window text
.plotFit <- function(resp, logc, pars, bline="bmad") {
###--------------------------- Draw Left Panel ----------------------------###
layout(mat=matrix(seq_len(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(vapply(pars, is.numeric, logical(1)))
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)
}
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)
}
}
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=expression(bold(paste("Concentration (",mu,"M)"))),
ylab=ylab,
main="",
bty="n",
xaxt="n",
yaxt="n",
type="n"
)
do.call(what=plot, args=c(resp[md] ~ logc[md], p), quote=TRUE)
useBmad <- bline != "coff" | is.null(pars$coff)
rect(
xleft=par()$usr[1],
xright=par()$usr[2],
ybottom=-1*ifelse(useBmad, 3*pars$bmad, pars$coff),
ytop=ifelse(useBmad, 3*pars$bmad, pars$coff),
border=NA,
col="gray70",
density=15,
angle=ifelse(useBmad, 45, -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 (!is.na(pars$hill) & pars$hill) {
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=TRUE,
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) {
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=TRUE,
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"
)
}
axis(
side=1,
at=axTicks(side=1),
labels=signif(10^axTicks(side=1), digits=2),
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[seq_len(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[seq_len(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
)
)
} else {
ftxt <- NULL
}
plot_txt1 <- paste0(itxt, htxt, gtxt, atxt, ntxt, ctxt, ftxt)
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"
)
)
}
#-------------------------------------------------------------------------------
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.