Nothing
# Plot of BMD values as an ECDF plot, with an horizontal color gradient coding
# for the signal level at each dose (or concentration)for each item,
# form an extended results dataframe (e.g. with annotation of items)
# with optionnal use of columns for shape and or facet
bmdplotwithgradient <- function(extendedres, BMDtype = c("zSD", "xfold"),
xmin, xmax, y0shift = TRUE, scaling = FALSE,
facetby, facetby2, shapeby, npoints = 50,
line.size, point.size = 1,
ncol4faceting, limits4colgradient,
lowercol = "darkblue", uppercol = "darkred",
add.label = FALSE, label.size = 2,
BMD_log_transfo = FALSE)
{
BMDtype <- match.arg(BMDtype, c("zSD", "xfold"))
if (missing(extendedres) | !is.data.frame(extendedres))
stop("The first argument of bmdplotwithgradient must be a dataframe
(see ?bmdplotwithgradient for details).")
if(scaling & !y0shift)
{
y0shift <- TRUE
warning(strwrap(prefix = "\n", initial = "\n",
"y0shift is forced to TRUE when scaling is TRUE."))
}
cnames <- colnames(extendedres)
if (BMDtype == "zSD")
{
if (scaling)
{
if (any(!is.element(c("id", "model", "b", "c", "d", "e", "f", "y0", "maxychange", "BMD.zSD"), cnames)))
stop("The first argument of bmdplotwithgradient must be a dataframe
containing at least columns named id, model, b, c, d, e, f, y0, maxychange and BMD.zSD.")
} else
{
if (any(!is.element(c("id", "model", "b", "c", "d", "e", "f", "y0", "BMD.zSD"), cnames)))
stop("The first argument of bmdplotwithgradient must be a dataframe
containing at least columns named id, model, b, c, d, e, f, y0 and BMD.zSD.")
}
BMD2plot <- data.frame(x = extendedres$BMD.zSD, id = extendedres$id)
}
else
{
if (scaling)
{
if (any(!is.element(c("id", "model", "b", "c", "d", "e", "f", "y0", "maxychange", "BMD.xfold"), cnames)))
stop("The first argument of bmdplotwithgradient must be a dataframe
containing at least columns named id, model, b, c, d, e, f, y0, maxychange and BMD.xfold.")
} else
{
if (any(!is.element(c("id", "model", "b", "c", "d", "e", "f", "y0", "BMD.xfold"), cnames)))
stop("The first argument of bmdplotwithgradient must be a dataframe
containing at least columns named id, model, b, c, d, e, f, y0 and BMD.xfold.")
}
BMD2plot <- data.frame(x = extendedres$BMD.xfold, id = extendedres$id)
}
if (BMD_log_transfo)
{
if (missing(xmin))
{
xmin <- min(BMD2plot$x[is.finite(BMD2plot$x) & BMD2plot$x != 0])
} else
{
if (xmin == 0)
{
warning(strwrap(prefix = "\n", initial = "\n",
"When using a log scale for the BMD plot, it is not possible to fix xmin at 0.
If the default value does not suit you, you can define a strictly positive value for xmin."))
xmin <- min(BMD2plot$x[is.finite(BMD2plot$x) & BMD2plot$x != 0])
}
}
} else
{
if (missing(xmin)) xmin <- 0
}
if (missing(xmax))
{
warning(strwrap(prefix = "\n", initial = "\n",
"We recommend you to define xmax at the maximal tested dose."))
xmax <- max(BMD2plot$x[is.finite(BMD2plot$x)])
}
if (!missing(shapeby))
{
if (!is.character(shapeby))
stop("shapeby should be a character string for the name of the column coding for the point shape.")
if (!is.element(shapeby, cnames))
stop("shapeby should be a character string corresponding to the name of a column of
extendedres, the dataframe given in input.")
BMD2plot$shapeby <- extendedres[, shapeby]
}
# calculation of ECDF by facetby
ntot <- nrow(BMD2plot)
if (!missing(facetby))
{
if (!is.character(facetby))
stop("facetby should be a character string for the name of the column used for facetting.")
if (!is.element(facetby, cnames))
stop("facetby should be a character string corresponding to the name of a column of
extendedres, the dataframe given in input.")
BMD2plot$facetby <- extendedres[, facetby]
if (!missing(facetby2))
{
if (!is.character(facetby2))
stop("facetby2 should be a character string for the name of the column used for facetting.")
if (!is.element(facetby2, cnames))
stop("facetby2 should be a character string corresponding to the name of a column of
extendedres, the dataframe given in input.")
BMD2plot$facetby2 <- extendedres[, facetby2]
BMD2plot$group <- paste(extendedres[, facetby], extendedres[, facetby2], sep = "_")
} else
{
BMD2plot$group <- BMD2plot$facetby
}
if (missing(line.size)) line.size <- 24 / max(table(BMD2plot$group))
uniqueby <- unique(BMD2plot$group)
n.uniqueby <- length(uniqueby)
BMD2plot$ECDF <- rep(0, ntot) # initialization
for (i in 1:n.uniqueby)
{
indi <- which(BMD2plot$group == uniqueby[i])
ntoti <- length(indi)
BMD2plot$ECDF[indi] <- (rank(BMD2plot$x[indi], ties.method = "first") - 0.5) / ntoti
}
g <- ggplot(data = BMD2plot, mapping = aes_(x = quote(x), y = quote(ECDF),
label = quote(id)))
if (missing(facetby2)) g <- g + facet_wrap(~ facetby) else g <- g + facet_grid(facetby2 ~ facetby)
} else
{
if (missing(line.size)) line.size <- 24 / nrow(BMD2plot)
BMD2plot$ECDF <- (rank(BMD2plot$x, ties.method = "first") - 0.5) / ntot
g <- ggplot(data = BMD2plot, mapping = aes_(x = quote(x), y = quote(ECDF),
label = quote(id)))
}
# Calculation of theoretical signal to color the lines
if (BMD_log_transfo)
{
x2plot <- 10^seq(log10(xmin), log10(xmax), length.out = npoints)
} else
{
x2plot <- seq(xmin, xmax, length.out = npoints)
}
ns <- nrow(extendedres)
N <- ns * npoints
curves2plot <- data.frame(x = rep(x2plot, ns),
id = rep(BMD2plot$id, each = npoints),
ECDF = rep(BMD2plot$ECDF, each = npoints),
signal = numeric(length = N))
for (i in 1:ns)
{
modeli <- extendedres$model[i]
if (modeli == "linear")
{
b <- extendedres$b[i]
d <- extendedres$d[i]
curves2plot$signal[(i-1)*npoints + 1:npoints] <- flin(x2plot, b = extendedres$b[i],
d = extendedres$d[i])
} else
if (modeli == "exponential")
{
curves2plot$signal[(i-1)*npoints + 1:npoints] <- fExpo(x2plot, b = extendedres$b[i],
d = extendedres$d[i], e = extendedres$e[i])
} else
if (modeli == "Hill")
{
curves2plot$signal[(i-1)*npoints + 1:npoints] <- fHill(x2plot, b = extendedres$b[i], c = extendedres$c[i],
d = extendedres$d[i], e = extendedres$e[i])
} else
if (modeli == "Gauss-probit")
{
curves2plot$signal[(i-1)*npoints + 1:npoints] <- fGauss5p(x2plot, b = extendedres$b[i], c = extendedres$c[i],
d = extendedres$d[i], e = extendedres$e[i],
f = extendedres$f[i])
} else
if (modeli == "log-Gauss-probit")
{
curves2plot$signal[(i-1)*npoints + 1:npoints] <- fLGauss5p(x2plot, b = extendedres$b[i], c = extendedres$c[i],
d = extendedres$d[i], e = extendedres$e[i],
f = extendedres$f[i])
}
if (y0shift)
{
if (scaling)
{
curves2plot$signal[(i-1)*npoints + 1:npoints] <-
(curves2plot$signal[(i-1)*npoints + 1:npoints] - extendedres$y0[i]) / extendedres$maxychange[i]
} else
{
curves2plot$signal[(i-1)*npoints + 1:npoints] <-
curves2plot$signal[(i-1)*npoints + 1:npoints] - extendedres$y0[i]
}
}
}
# no shape no facet
if (!missing(facetby))
{
curves2plot$facetby <- rep(extendedres[, facetby], each = npoints)
if (!missing(facetby2))
{
curves2plot$facetby2 <- rep(extendedres[, facetby2], each = npoints)
}
}
gg <- g + geom_line(data = curves2plot,
mapping = aes_(x = quote(x), y = quote(ECDF),
group = quote(id), color = quote(signal)),
size = line.size)
if (!missing(facetby))
{
if (!missing(facetby2))
{
gg <- gg + facet_grid(facetby2 ~ facetby)
}
else
{
if (missing(ncol4faceting))
{
gg <- gg + facet_wrap(~ facetby)
} else
{
gg <- gg + facet_wrap(~ facetby, ncol = ncol4faceting)
}
}
}
# Add of the color gradient
if (missing(limits4colgradient))
{
gg <- gg +
# scale_colour_gradient2(low = "darkblue", mid = "white",
# high = "darkred", midpoint = median(curves2plot$signal), space = "Lab",
# na.value = "grey50", guide = "colourbar", aesthetics = "colour")
scale_colour_gradient2(low = lowercol, mid = "white",
high = uppercol, midpoint = 0, space = "Lab",
na.value = "grey50", guide = "colourbar", aesthetics = "colour")
} else
{
gg <- gg +
scale_colour_gradient2(low = lowercol, mid = "white",
high = uppercol, midpoint = 0, space = "Lab",
na.value = "grey50", guide = "colourbar", aesthetics = "colour",
limits = limits4colgradient)
}
# Add of points (BMD values)
if (!missing(shapeby))
{
gg <- gg + geom_point(data = BMD2plot, mapping = aes_(shape = quote(shapeby)),
size = point.size)
} else
{
gg <- gg + geom_point(data = BMD2plot, size = point.size)
}
gg <- gg + theme_classic()
if(add.label)
{
if (!missing(shapeby))
warning(strwrap(prefix = "\n", initial = "\n",
"The type of points will not be seen when points are replaced by labels.
You should omit it in this case."))
gg <- gg + geom_label(size = label.size)
}
if (BMD_log_transfo)
gg <- gg + scale_x_log10()
gg <- gg + xlab("BMD")
return(gg)
}
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.