# Re: [R] Rcmdr and scatter3d

From: John Fox <jfox_at_mcmaster.ca>
Date: Wed 05 Oct 2005 - 11:18:56 EST

Dear Naiara,

Combine the data sets and differentiate among them with a factor. Then use the groups argument to scatter3d (see ?scatter3d). If you're using the R Commander to make the plot, the 3D scatterplot dialog box as a plot by groups button. You can also fit colour-coded regression surfaces by group.

I've appended a new version of the scatter3d function, not yet in the Rcmdr package, which will also plot data ellipsoids (for the whole data set or by groups).

I hope this helps,
John

• snip --------------

ellipsoid <- function(center=c(0, 0, 0), radius=1, shape=diag(3), n=30){ # adapted from the shapes3d demo in the rgl package   degvec <- seq(0, 2*pi, length=n)
ecoord2 <- function(p) c(cos(p[1])*sin(p[2]), sin(p[1])*sin(p[2]), cos(p[2]))

```  v <- t(apply(expand.grid(degvec,degvec), 1, ecoord2))
v <- center + radius * t(v %*% chol(shape))
v <- rbind(v, rep(1,ncol(v)))
e <- expand.grid(1:(n-1), 1:n)
i1 <- apply(e, 1, function(z) z[1] + n*(z[2] - 1))
i2 <- i1 + 1
```

i3 <- (i1 + n - 1) %% n^2 + 1
i4 <- (i2 + n - 1) %% n^2 + 1
i <- rbind(i1, i2, i4, i3)
qmesh3d(v, i)
}

scatter3d <- function(x, y, z, xlab=deparse(substitute(x)), ylab=deparse(substitute(y)),

zlab=deparse(substitute(z)), revolutions=0, bg.col=c("white", "black"),     axis.col=if (bg.col == "white") "black" else "white",     surface.col=c("blue", "green", "orange", "magenta", "cyan", "red", "yellow", "gray"),

neg.res.col="red", pos.res.col="green", 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,

sphere.size=1, threshold=0.01, speed=1, fov=60,     fit="linear", groups=NULL, parallel=TRUE, ellipsoid=FALSE, level=0.5,     model.summary=FALSE){
require(rgl)
require(mgcv)
summaries <- list()
if ((!is.null(groups)) && (nlevels(groups) > length(surface.col)))

```        stop(sprintf(
gettextRcmdr("Number of groups (%d) exceeds number of colors
(%d)."),
nlevels(groups), length(surface.col)))
if ((!is.null(groups)) && (!is.factor(groups)))
stop(gettextRcmdr("groups variable must be a factor."))
```
bg.col <- match.arg(bg.col)
fogtype <- match.arg(fogtype)
if ((length(fit) > 1) && residuals && surface)

stop(gettextRcmdr("cannot plot both multiple surfaces and residuals"))

xlab # cause these arguments to be evaluated     ylab
zlab

```    rgl.clear()
rgl.viewpoint(fov=fov)
rgl.bg(col=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]
if (!is.null(groups)) groups <- groups[valid]
```    x <- (x - min(x))/(max(x) - min(x))
y <- (y - min(y))/(max(y) - min(y))
z <- (z - min(z))/(max(z) - min(z))
```

size <- sphere.size*((100/length(x))^(1/3))*0.015     if (is.null(groups)){

if (size > threshold) rgl.spheres(x, y, z, color=point.col, radius=size)

```            else rgl.points(x, y, z, color=point.col)

}
else {
if (size > threshold) rgl.spheres(x, y, z,
else rgl.points(x, y, z, color=surface.col[as.numeric(groups)])
}
rgl.lines(c(0,1), c(0,0), c(0,0), color=axis.col)
rgl.lines(c(0,0), c(0,1), c(0,0), color=axis.col)
rgl.lines(c(0,0), c(0,0), c(0,1), color=axis.col)
rgl.texts(1, 0, 0, xlab, adj=1, color=text.col)
```
rgl.texts(0, 1, 0, ylab, adj=1, color=text.col)     rgl.texts(0, 0, 1, zlab, adj=1, color=text.col)     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)),
lit=FALSE)
if (grid) 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)),
lit=FALSE)
if (grid) wire3d(ellips, col=surface.col[j], lit=FALSE)
coords <- ellips\$vb[, which.max(ellips\$vb[1,])]
if (!surface) rgl.texts(coords[1] + 0.05, coords[2],
coords[3], group,
col=surface.col[j])
}

}
}
if (surface){
vals <- seq(0, 1, length=grid.lines)
dat <- expand.grid(x=vals, z=vals)
for (i in 1:length(fit)){
f <- match.arg(fit[i], c("linear", "quadratic", "smooth",
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)) gam(y ~ s(x, z))
else gam(y ~ s(x, z, fx=TRUE, k=df.smooth)),
s(z))
else gam(y ~ s(x, fx=TRUE, k=df.additive[1]+1) +
)
if (model.summary) summaries[[f]] <- summary(mod)
yhat <- matrix(predict(mod, newdata=dat), grid.lines,
grid.lines)
if (fill) rgl.surface(vals, vals, yhat,
color=surface.col[i],
alpha=0.5, lit=FALSE)
if(grid) 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.lines(as.vector(rbind(x,x)),
as.vector(rbind(y,fitted)),
as.vector(rbind(z,z)),
color=as.vector(rbind(colors,colors)))
}
}
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)) gam(y ~ s(x, z) +
groups)
else gam(y ~ s(x, z, fx=TRUE, k=df.smooth) +
groups),
s(z) + groups)
else gam(y ~ s(x, fx=TRUE, k=df.additive[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.surface(vals, vals, yhat,
color=surface.col[j],
alpha=0.5, lit=FALSE)
if (grid) rgl.surface(vals, vals, yhat, color=if
(fill) grid.col
else surface.col[j], alpha=0.5, lit=FALSE,
front="lines", back="lines")
rgl.texts(0, predict(mod, newdata=data.frame(x=0,
z=0,
groups=group)), 0,
if (residuals){
yy <- y[select.obs]
xx <- x[select.obs]
zz <- z[select.obs]
fitted <- fitted(mod)[select.obs]
rgl.lines(as.vector(rbind(xx,xx)),
as.vector(rbind(yy,fitted)),
as.vector(rbind(zz,zz)),
col=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)) gam(y ~ s(x,
z),
subset=select.obs)
else gam(y ~ s(x, z, fx=TRUE, k=df.smooth),
subset=select.obs),
s(x) + s(z),
subset=select.obs)
else gam(y ~ s(x, fx=TRUE,
s(z, fx=TRUE,
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.surface(vals, vals, yhat,
color=surface.col[j],
alpha=0.5, lit=FALSE)
if (grid) rgl.surface(vals, vals, yhat, color=if
(fill) grid.col
else surface.col[j], alpha=0.5, lit=FALSE,
front="lines", back="lines")
rgl.texts(0, predict(mod,
newdata=data.frame(x=0, z=0, groups=group)), 0,
if (residuals){
yy <- y[select.obs]
xx <- x[select.obs]
zz <- z[select.obs]
fitted <- fitted(mod)
rgl.lines(as.vector(rbind(xx,xx)),
as.vector(rbind(yy,fitted)),
as.vector(rbind(zz,zz)),
col=surface.col[j])
}
}
}
}

}
}
if (revolutions > 0) {
for (i in 1:revolutions){
for (angle in seq(1, 360, length=360/speed))
rgl.viewpoint(-angle, fov=fov)

}
}
```

if (model.summary) return(summaries) else return(invisible(NULL))     }

John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
905-525-9140x23604
http://socserv.mcmaster.ca/jfox
```> -----Original Message-----
> From: r-help-bounces@stat.math.ethz.ch
> [mailto:r-help-bounces@stat.math.ethz.ch] On Behalf Of Naiara S. Pinto
> Sent: Tuesday, October 04, 2005 6:13 PM
> To: r-help@stat.math.ethz.ch
> Subject: [R] Rcmdr and scatter3d
>
> Hi folks,
>
> I'd like to use scatter3d (which is in R commander) to plot
> more than one dataset in the same graph, each dataset with a
> different color. The kind of stuff you would do with "holdon"
> in Matlab.
>
> I read a recent message that was posted to this list with a
> similar problem, but I couldn't understand the reply. Could
> someone give me one example? How do you plot subgroups using
> scatter3d?
>
> Thanks a lot!
>
> Naiara.
>
>
> --------------------------------------------
> Naiara S. Pinto
> Ecology, Evolution and Behavior
> 1 University Station A6700
> Austin, TX, 78712
>
> ______________________________________________
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help