##' Description ...
##'
##' Details ...
##'
##' @title Plot three variable in a 3D interactive graphics window
##'
##' @param x a parameter
##' @param y a parameter
##' @param z a parameter
##' @param xlab a parameter
##' @param ylab a parameter
##' @param axis.scales a parameter
##' @param zlab a parameter
##' @param revolutions a parameter
##' @param bg.col a parameter
##' @param axis.col a parameter
##' @param surface.col a parameter
##' @param neg.res.col a parameter
##' @param pos.res.col a parameter
##' @param square.col a parameter
##' @param point.col a parameter
##' @param text.col a parameter
##' @param grid.col a parameter
##' @param fogtype a parameter
##' @param residuals a parameter
##' @param surface a parameter
##' @param fill a parameter
##' @param grid a parameter
##' @param grid.lines a parameter
##' @param df.smooth a parameter
##' @param df.additive a parameter
##' @param sphere.size a parameter
##' @param threshold a parameter
##' @param speed a parameter
##' @param fov a parameter
##' @param fit a parameter
##' @param groups a parameter
##' @param parallel a parameter
##' @param ellipsoid a parameter
##' @param level a parameter
##' @param model.summary a parameter
##'
##' @import gWidgets2 gWidgets2RGtk2
##'
##' @return NULL
scatter3d <- function(x, y, z,
xlab=deparse(substitute(x)), ylab=deparse(substitute(y)),
axis.scales=TRUE,
zlab=deparse(substitute(z)), revolutions=0, bg.col=c("white", "black"),
axis.col=if (bg.col == "white") c("darkmagenta", "black", "darkcyan")
else c("darkmagenta", "white", "darkcyan"),
surface.col=c("blue", "green", "orange", "magenta", "cyan", "red", "yellow", "gray"),
neg.res.col="red", pos.res.col="green",
square.col=if (bg.col == "white") "black" else "gray", point.col="yellow",
text.col=axis.col, grid.col=if (bg.col == "white") "black" else "gray",
fogtype=c("exp2", "linear", "exp", "none"),
residuals=(length(fit) == 1), surface=TRUE, fill=TRUE, grid=TRUE, grid.lines=26,
df.smooth=NULL, df.additive=NULL,
sphere.size=1, threshold=0.01, speed=1, fov=60,
fit="linear", groups=NULL, parallel=TRUE, ellipsoid=FALSE, level=0.5,
model.summary = FALSE) {
## if (!"package:rgl" %in% search()) {
## ## Load "rgl" only when the "scatter3d" function is called ...
## ## ... but to get around R CMD CHECK's forceful behavious, we'll be clever ...
## cmd <- "require(rgl)"
## if (!eval(parse(text = cmd)))
## stop("Unable to load the RGL package ... is it installed?")
## }
if (!requireNamespace("rgl", quietly = TRUE)) {
gmessage("Please install the `rgl` package to use this module.")
return(invisible(NULL))
}
use.gams <- requireNamespace("mgcv", quietly = TRUE)
if (residuals == "squares"){
residuals <- TRUE
squares <- TRUE
}
else squares <- FALSE
summaries <- list()
if ((!is.null(groups)) && (nlevels(groups) > length(surface.col)))
stop(sprintf(gettext("Number of groups (%d) exceeds number of colours (%d)."),
nlevels(groups), length(surface.col)))
if ((!is.null(groups)) && (!is.factor(groups))) stop(gettext("groups variable must be a factor."))
bg.col <- match.arg(bg.col)
fogtype <- match.arg(fogtype)
if ((length(fit) > 1) && residuals && surface)
stop(gettext("cannot plot both multiple surfaces and residuals"))
## xlab cause these arguments to be evaluated
## ylab
## zlab
rgl::rgl.clear()
rgl::rgl.viewpoint(fov=fov)
rgl::rgl.bg(color=bg.col, fogtype=fogtype)
valid <- if (is.null(groups)) complete.cases(x, y, z)
else complete.cases(x, y, z, groups)
x <- x[valid]
y <- y[valid]
z <- z[valid]
minx <- min(x)
maxx <- max(x)
miny <- min(y)
maxy <- max(y)
minz <- min(z)
maxz <- max(z)
if (axis.scales){
## lab.min.x <- nice(minx)
## lab.max.x <- nice(maxx)
## lab.min.y <- nice(miny)
## lab.max.y <- nice(maxy)
## lab.min.z <- nice(minz)
## lab.max.z <- nice(maxz)
lab.min.x <- minx
lab.max.x <- maxx
lab.min.y <- miny
lab.max.y <- maxy
lab.min.z <- minz
lab.max.z <- maxz
minx <- min(lab.min.x, minx)
maxx <- max(lab.max.x, maxx)
miny <- min(lab.min.y, miny)
maxy <- max(lab.max.y, maxy)
minz <- min(lab.min.z, minz)
maxz <- max(lab.max.z, maxz)
min.x <- (lab.min.x - minx)/(maxx - minx)
max.x <- (lab.max.x - minx)/(maxx - minx)
min.y <- (lab.min.y - miny)/(maxy - miny)
max.y <- (lab.max.y - miny)/(maxy - miny)
min.z <- (lab.min.z - minz)/(maxz - minz)
max.z <- (lab.max.z - minz)/(maxz - minz)
}
if (!is.null(groups)) groups <- groups[valid]
x <- (x - minx)/(maxx - minx)
y <- (y - miny)/(maxy - miny)
z <- (z - minz)/(maxz - minz)
size <- sphere.size*((100/length(x))^(1/3))*0.015
if (is.null(groups)){
if (size > threshold) rgl::rgl.spheres(x, y, z, color=point.col, radius=size)
else rgl::rgl.points(x, y, z, color=point.col)
}
else {
if (size > threshold) rgl::rgl.spheres(x, y, z, color=surface.col[as.numeric(groups)], radius=size)
else rgl::rgl.points(x, y, z, color=surface.col[as.numeric(groups)])
}
if (!axis.scales) axis.col[1] <- axis.col[3] <- axis.col[2]
rgl::rgl.lines(c(0,1), c(0,0), c(0,0), color=axis.col[1])
rgl::rgl.lines(c(0,0), c(0,1), c(0,0), color=axis.col[2])
rgl::rgl.lines(c(0,0), c(0,0), c(0,1), color=axis.col[3])
rgl::rgl.texts(1, 0, 0, xlab, adj=1, color=axis.col[1])
rgl::rgl.texts(0, 1.05, 0, ylab, adj=1, color=axis.col[2])
rgl::rgl.texts(0, 0, 1, zlab, adj=1, color=axis.col[3])
if (axis.scales){
rgl::rgl.texts(min.x, -0.05, 0, lab.min.x, col=axis.col[1])
rgl::rgl.texts(max.x, -0.05, 0, lab.max.x, col=axis.col[1])
rgl::rgl.texts(0, -0.1, min.z, lab.min.z, col=axis.col[3])
rgl::rgl.texts(0, -0.1, max.z, lab.max.z, col=axis.col[3])
rgl::rgl.texts(-0.05, min.y, -0.05, lab.min.y, col=axis.col[2])
rgl::rgl.texts(-0.05, max.y, -0.05, lab.max.y, col=axis.col[2])
}
if (ellipsoid) {
dfn <- 3
if (is.null(groups)){
dfd <- length(x) - 1
radius <- sqrt(dfn * qf(level, dfn, dfd))
ellips <- ellipsoid(center=c(mean(x), mean(y), mean(z)),
shape=cov(cbind(x,y,z)), radius=radius)
if (fill) rgl::shade3d(ellips, col=surface.col[1], alpha=0.1, lit=FALSE)
if (grid) rgl::wire3d(ellips, col=surface.col[1], lit=FALSE)
}
else{
levs <- levels(groups)
for (j in 1:length(levs)){
group <- levs[j]
select.obs <- groups == group
xx <- x[select.obs]
yy <- y[select.obs]
zz <- z[select.obs]
dfd <- length(xx) - 1
radius <- sqrt(dfn * qf(level, dfn, dfd))
ellips <- ellipsoid(center=c(mean(xx), mean(yy), mean(zz)),
shape=cov(cbind(xx,yy,zz)), radius=radius)
if (fill) rgl::shade3d(ellips, col=surface.col[j], alpha=0.1, lit=FALSE)
if (grid) rgl::wire3d(ellips, col=surface.col[j], lit=FALSE)
coords <- ellips$vb[, which.max(ellips$vb[1,])]
if (!surface) rgl::rgl.texts(coords[1] + 0.05, coords[2], coords[3], group,
col=surface.col[j])
}
}
}
if (surface){
vals <- seq(0, 1, length.out=grid.lines)
dat <- expand.grid(x=vals, z=vals)
for (i in 1:length(fit)){
f <- match.arg(fit[i], c("linear", "quadratic", "smooth", "additive"))
if (! use.gams & f %in% c("smooth", "additive")) {
f <- "linear"
gmessage("To use 'smooth' or 'additive' fits you must install the 'mgcv' package for R.\n\nUsing a linear fit instead.",
title = "Error - mgcv not installed", icon = "error")
}
if (is.null(groups)){
mod <- switch(f,
linear = lm(y ~ x + z),
quadratic = lm(y ~ (x + z)^2 + I(x^2) + I(z^2)),
smooth = if (is.null(df.smooth)) mgcv::gam(y ~ s(x, z))
else mgcv::gam(y ~ s(x, z, fx=TRUE, k=df.smooth)),
additive = if (is.null(df.additive)) mgcv::gam(y ~ s(x) + s(z))
else mgcv::gam(y ~ s(x, fx=TRUE, k=df.additive[1]+1) +
s(z, fx=TRUE, k=(rev(df.additive+1)[1]+1)))
)
if (model.summary) summaries[[f]] <- summary(mod)
yhat <- matrix(predict(mod, newdata=dat), grid.lines, grid.lines)
if (fill) rgl::rgl.surface(vals, vals, yhat, color=surface.col[i], alpha=0.5, lit=FALSE)
if(grid) rgl::rgl.surface(vals, vals, yhat, color=if (fill) grid.col
else surface.col[i], alpha=0.5, lit=FALSE, front="lines", back="lines")
if (residuals){
n <- length(y)
fitted <- fitted(mod)
colors <- ifelse(residuals(mod) > 0, pos.res.col, neg.res.col)
rgl::rgl.lines(as.vector(rbind(x,x)), as.vector(rbind(y,fitted)), as.vector(rbind(z,z)),
color=as.vector(rbind(colors,colors)))
if (squares){
res <- y - fitted
xx <- as.vector(rbind(x, x, x + res, x + res))
yy <- as.vector(rbind(y, fitted, fitted, y))
zz <- as.vector(rbind(z, z, z, z))
rgl::rgl.quads(xx, yy, zz, color=square.col, alpha=0.5, lit=FALSE)
rgl::rgl.lines(xx, yy, zz, color=square.col)
}
}
}
else{
if (parallel){
mod <- switch(f,
linear = lm(y ~ x + z + groups),
quadratic = lm(y ~ (x + z)^2 + I(x^2) + I(z^2) + groups),
smooth = if (is.null(df.smooth)) mgcv::gam(y ~ s(x, z) + groups)
else mgcv::gam(y ~ s(x, z, fx=TRUE, k=df.smooth) + groups),
additive = if (is.null(df.additive)) mgcv::gam(y ~ s(x) + s(z) + groups)
else mgcv::gam(y ~ s(x, fx=TRUE, k=df.additive[1]+1) +
s(z, fx=TRUE, k=(rev(df.additive+1)[1]+1)) + groups)
)
if (model.summary) summaries[[f]] <- summary(mod)
levs <- levels(groups)
for (j in 1:length(levs)){
group <- levs[j]
select.obs <- groups == group
yhat <- matrix(predict(mod, newdata=cbind(dat, groups=group)), grid.lines, grid.lines)
if (fill) rgl::rgl.surface(vals, vals, yhat, color=surface.col[j], alpha=0.5, lit=FALSE)
if (grid) rgl::rgl.surface(vals, vals, yhat, color=if (fill) grid.col
else surface.col[j], alpha=0.5, lit=FALSE, front="lines", back="lines")
rgl::rgl.texts(1, predict(mod, newdata=data.frame(x=1, z=1, groups=group, stringsAsFactors = TRUE)), 1,
paste(group, " "), adj=1, color=surface.col[j])
if (residuals){
yy <- y[select.obs]
xx <- x[select.obs]
zz <- z[select.obs]
fitted <- fitted(mod)[select.obs]
res <- yy - fitted
rgl::rgl.lines(as.vector(rbind(xx,xx)), as.vector(rbind(yy,fitted)), as.vector(rbind(zz,zz)),
col=surface.col[j])
if (squares) {
xxx <- as.vector(rbind(xx, xx, xx + res, xx + res))
yyy <- as.vector(rbind(yy, fitted, fitted, yy))
zzz <- as.vector(rbind(zz, zz, zz, zz))
rgl::rgl.quads(xxx, yyy, zzz, color=surface.col[j], alpha=0.5, lit=FALSE)
rgl::rgl.lines(xxx, yyy, zzz, color=surface.col[j])
}
}
}
}
else {
levs <- levels(groups)
for (j in 1:length(levs)){
group <- levs[j]
select.obs <- groups == group
mod <- switch(f,
linear = lm(y ~ x + z, subset=select.obs),
quadratic = lm(y ~ (x + z)^2 + I(x^2) + I(z^2), subset=select.obs),
smooth = if (is.null(df.smooth)) mgcv::gam(y ~ s(x, z), subset=select.obs)
else mgcv::gam(y ~ s(x, z, fx=TRUE, k=df.smooth), subset=select.obs),
additive = if (is.null(df.additive)) mgcv::gam(y ~ s(x) + s(z), subset=select.obs)
else mgcv::gam(y ~ s(x, fx=TRUE, k=df.additive[1]+1) +
s(z, fx=TRUE, k=(rev(df.additive+1)[1]+1)), subset=select.obs)
)
if (model.summary) summaries[[paste(f, ".", group, sep="")]] <- summary(mod)
yhat <- matrix(predict(mod, newdata=dat), grid.lines, grid.lines)
if (fill) rgl::rgl.surface(vals, vals, yhat, color=surface.col[j], alpha=0.5, lit=FALSE)
if (grid) rgl::rgl.surface(vals, vals, yhat, color=if (fill) grid.col
else surface.col[j], alpha=0.5, lit=FALSE, front="lines", back="lines")
rgl::rgl.texts(1, predict(mod, newdata=data.frame(x=1, z=1, groups=group, stringsAsFactors = TRUE)), 1,
paste(group, " "), adj=1, color=surface.col[j])
if (residuals){
yy <- y[select.obs]
xx <- x[select.obs]
zz <- z[select.obs]
fitted <- fitted(mod)
res <- yy - fitted
rgl::rgl.lines(as.vector(rbind(xx,xx)), as.vector(rbind(yy,fitted)), as.vector(rbind(zz,zz)),
col=surface.col[j])
if (squares) {
xxx <- as.vector(rbind(xx, xx, xx + res, xx + res))
yyy <- as.vector(rbind(yy, fitted, fitted, yy))
zzz <- as.vector(rbind(zz, zz, zz, zz))
rgl::rgl.quads(xxx, yyy, zzz, color=surface.col[j], alpha=0.5, lit=FALSE)
rgl::rgl.lines(xxx, yyy, zzz, color=surface.col[j])
}
}
}
}
}
}
}
if (revolutions > 0) {
for (i in 1:revolutions){
for (angle in seq(1, 360, length.out=360/speed)) rgl::rgl.viewpoint(-angle, fov=fov)
}
}
if (model.summary) return(summaries) else return(invisible(NULL))
}
Rcmdr.select3d <-
function (...)
{
rgl::.check3d()
rect <- rgl::rgl.select(...)
llx <- rect[1]
lly <- rect[2]
urx <- rect[3]
ury <- rect[4]
if (llx > urx) {
temp <- llx
llx <- urx
urx <- temp
}
if (lly > ury) {
temp <- lly
lly <- ury
ury <- temp
}
proj <- rgl::rgl.projection()
function(x, y, z) {
pixel <- rgl::rgl.user2window(x, y, z, projection = proj)
apply(pixel, 1, function(p) (llx <= p[1]) && (p[1] <=
urx) && (lly <= p[2]) && (p[2] <= ury) && (0 <= p[3]) &&
(p[3] <= 1))
}
}
identify3d <-
function (x, y, z, axis.scales=TRUE, groups = NULL, labels = 1:length(x),
col = c("blue", "green", "orange", "magenta", "cyan", "red", "yellow", "gray"),
offset = ((100/length(x))^(1/3)) * 0.02)
{
valid <- if (is.null(groups))
complete.cases(x, y, z)
else complete.cases(x, y, z, groups)
labels <- labels[valid]
x <- x[valid]
y <- y[valid]
z <- z[valid]
minx <- min(x)
maxx <- max(x)
miny <- min(y)
maxy <- max(y)
minz <- min(z)
maxz <- max(z)
if (axis.scales){
# lab.min.x <- nice(minx)
# lab.max.x <- nice(maxx)
# lab.min.y <- nice(miny)
# lab.max.y <- nice(maxy)
# lab.min.z <- nice(minz)
# lab.max.z <- nice(maxz)
lab.min.x <- minx
lab.max.x <- maxx
lab.min.y <- miny
lab.max.y <- maxy
lab.min.z <- minz
lab.max.z <- maxz
minx <- min(lab.min.x, minx)
maxx <- max(lab.max.x, maxx)
miny <- min(lab.min.y, miny)
maxy <- max(lab.max.y, maxy)
minz <- min(lab.min.z, minz)
maxz <- max(lab.max.z, maxz)
min.x <- (lab.min.x - minx)/(maxx - minx)
max.x <- (lab.max.x - minx)/(maxx - minx)
min.y <- (lab.min.y - miny)/(maxy - miny)
max.y <- (lab.max.y - miny)/(maxy - miny)
min.z <- (lab.min.z - minz)/(maxz - minz)
max.z <- (lab.max.z - minz)/(maxz - minz)
}
x <- (x - minx)/(maxx - minx)
y <- (y - miny)/(maxy - miny)
z <- (z - minz)/(maxz - minz)
rgl::rgl.bringtotop()
identified <- character(0)
groups <- if (!is.null(groups))
as.numeric(groups[valid])
else rep(1, length(x))
repeat {
f <- Rcmdr.select3d(button="right")
which <- f(x, y, z)
if (!any(which))
break
rgl::rgl.texts(x[which], y[which] + offset, z[which], labels[which],
color = col[groups][which])
identified <- c(identified, labels[which])
}
unique(identified)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.