#' Graphic output of the PAMM function
#'
#' Provide graphic interpretation of the simulation results
#'
#' @param x a PAMM object
#' @param graphtype "group", "repl" or "both"
#' "group" give graphs with varying number of ID and with a fixed number
#' of replicates specified in nbrepl
#' "repl" give graphs with varying number of replicates and with a fixed
#' number of ID specified in nbgroup
#' "both" 3-D plot. see also fun3D argument. Note: useful only with multiple
#' group size and multiple number of replicates.
#'
#' @param nbgroup number of group for which plots the output. Necessary for "repl" type of graph
#' @param nbrepl number of replicates for which plots the output. Necessary for "group" type of graph
#' @param fun3D plot function used to plot the 3D graph.
#' "wireframe" uses lattice,
#' "persp" uses graphics
#' and "open3d" uses rgl
#' @param \dots potentially further arguments to pass to methods
#'
#'
#' @details
#' Parameters phi, theta, ticktype ("simple" or "detailed"), nticks, nbcol
#' from [persp()] function could also be specified for 3D plots.
#' In addition, color schemes ("grey", "cm.colors" and "rainbow") and
#' coltype ("restricted" or "0-1") parameters could also be specified for
#' 3D plots.
#'
#' @examples
#' \dontrun{
#' ours <- PAMM(numsim=10,group=c(seq(10,50,10),100),repl=c(3,4,6),
#' randompart=c(0.4,0.1,0.5,0.1),fixed=c(0,1,0.7))
#' plot(ours, "both")
#' plot(ours, "group",nbrepl=4)
#' plot(ours,"repl",nbgroup=20)
#' }
#'
#' @export
plot.PAMM <- function(x, graphtype = "both", nbgroup, nbrepl, fun3D = "wireframe",
...) {
if (!inherits(x, "PAMM"))
stop("use only with \"PAMM\" objects")
extra <- list(...) #not doing anything for the moment
xs <- c("int.pval", "sl.pval", "int.power", "sl.power")
ys <- c("ipv", "slpv", "ipo", "slpo")
labs <- c("int.p-value", "slope.p-value", "int.power", "slope.power")
mains <- c("Intercept P-value", "Slope P-value", "Intercept Power", "Slope Power")
if (graphtype == "group") {
par(mfrow = c(2, 2), mar = c(2, 2, 2, 2), oma = c(4, 4, 2, 0))
for (i in 1:4) {
plot(x$nb.ID[x$nb.repl == nbrepl], x[x$nb.repl == nbrepl, xs[i]], type = "l",
xlab = "", ylab = "", ylim = c(0, 1), bty = "L")
if (i <= 2)
abline(h = 0.05)
lines(x$nb.ID[x$nb.repl == nbrepl], x[x$nb.repl == nbrepl, paste("CIup.",
ys[i], sep = "")], lty = 2)
lines(x$nb.ID[x$nb.repl == nbrepl], x[x$nb.repl == nbrepl, paste("CIlow.",
ys[i], sep = "")], lty = 2)
}
mtext("Number of group", 1, line = 2, outer = TRUE, cex = 1.2)
mtext(c("P-value", "Power"), 2, at = c(0.75, 0.25), line = 2, outer = TRUE,
cex = 1.2)
mtext(c("Intercept", "Slope"), 3, at = c(0.25, 0.75), line = 0.5, outer = TRUE,
cex = 1.2)
}
if (graphtype == "repl") {
par(mfrow = c(2, 2), mar = c(2, 2, 2, 2), oma = c(4, 4, 2, 0))
for (i in 1:4) {
plot(x$nb.repl[x$nb.ID == nbgroup], x[x$nb.ID == nbgroup, xs[i]], type = "l",
xlab = "", ylab = "", ylim = c(0, 1), bty = "L")
if (i <= 2)
abline(h = 0.05)
lines(x$nb.repl[x$nb.ID == nbgroup], x[x$nb.ID == nbgroup, paste("CIup.",
ys[i], sep = "")], lty = 2)
lines(x$nb.repl[x$nb.ID == nbgroup], x[x$nb.ID == nbgroup, paste("CIlow.",
ys[i], sep = "")], lty = 2)
}
mtext("Number of replicates", 1, line = 2, outer = TRUE, cex = 1.2)
mtext(c("P-value", "Power"), 2, at = c(0.75, 0.25), line = 2, outer = TRUE,
cex = 1.2)
mtext(c("Intercept", "Slope"), 3, at = c(0.25, 0.75), line = 0.5, outer = TRUE,
cex = 1.2)
}
if (graphtype == "both") {
if (fun3D == "wireframe") {
par.set <- list(axis.line = list(col = "transparent"), clip = list(panel = "off"))
zl <- c("P-value", "Power", "P-value", "Power")
for (i in 1:4) {
p <- wireframe(as.formula(paste(xs[i], " ~ nb.ID + nb.repl")), x,
colorkey = FALSE, drape = TRUE, scales = list(arrows = FALSE, distance = c(2,
2, 2)), xlab = "Group", ylab = "Repl", main = mains[i], zlab = list(zl[i],
rot = 90), screen = list(z = -50, x = -70, y = 0), par.settings = par.set)
if (i == 1)
print(p, split = c(1, 1, 2, 2), more = TRUE) else if (i == 2)
print(p, split = c(1, 2, 2, 2), more = TRUE) else if (i == 3)
print(p, split = c(2, 1, 2, 2), more = TRUE) else if (i == 4)
print(p, split = c(2, 2, 2, 2))
}
}
if (fun3D == "persp") {
phi <- ifelse(is.null(extra$phi), 25, extra$phi)
theta <- ifelse(is.null(extra$theta), 30, extra$theta)
ticktype <- ifelse(is.null(extra$ticktype), "detailed", extra$ticktype)
nticks <- ifelse(is.null(extra$nticks), 4, extra$nticks)
nbcol <- ifelse(is.null(extra$nbcol), 100, extra$nbcol)
color <- ifelse(is.null(extra$color), "grey", extra$color)
coltype <- ifelse(is.null(extra$coltype), "restricted", extra$coltype)
par(mfrow = c(2, 2), mar = c(2, 2, 2, 2), oma = c(0, 4, 2, 0))
X <- unique(x$nb.ID)
Y <- unique(x$nb.repl)
if (color == "grey")
color <- grey(0:nbcol/nbcol) else if (color == "cm")
color <- cm.colors(nbcol) else if (color == "rainbow")
color <- rainbow(nbcol)
for (i in 1:4) {
Z <- matrix(x[, xs[i]], nrow = length(X), ncol = length(Y))
nrz <- nrow(Z)
ncz <- ncol(Z)
zfacet <- (Z[-1, -1] + Z[-1, -ncz] + Z[-nrz, -1] + Z[-nrz, -ncz])/4
if (coltype == "0-1")
facetcol <- facetcol <- cut(zfacet, seq(0, 1, 1/nbcol)) else if (coltype == "restricted")
facetcol <- cut(zfacet, nbcol)
persp(X, Y, Z, col = color[facetcol], zlim = c(0, 1), phi = phi,
theta = theta, ticktype = ticktype, nticks = nticks, xlab = "group",
ylab = "repl", zlab = "")
mtext(c("P-value", "Power"), 2, at = c(0.75, 0.25), line = 2, outer = TRUE,
cex = 1.2)
mtext(c("Intercept", "Slope"), 3, at = c(0.25, 0.75), line = 0.5,
outer = TRUE, cex = 1.2)
}
}
if (fun3D == "open3d") {
if (!requireNamespace("rgl", quietly = TRUE)) {
warning("rgl package needed for this function to work. Please install it. ",
call. = FALSE)
} else {
for (i in 1:4) {
rgl::open3d()
rgl::bg3d("white")
rgl::material3d(col = "white")
rgl::persp3d(unique(x$nb.ID), unique(x$nb.repl), matrix(x[, xs[i]],
nrow = length(unique(x$nb.ID)), ncol = length(unique(x$nb.repl))),
col = rainbow(10), box = FALSE, zlim = c(0, 1), xlab = "Nb group",
ylab = "Nb repl", zlab = labs[i])
}
}
}
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.