[R] Rcmdr and scatter3d
John Fox
jfox at mcmaster.ca
Wed Oct 5 03:18:56 CEST 2005
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,
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){
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,
color=surface.col[as.numeric(groups)], radius=size)
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)),
shape=cov(cbind(x,y,z)), radius=radius)
if (fill) shade3d(ellips, col=surface.col[1], alpha=0.1,
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)),
shape=cov(cbind(xx,yy,zz)), radius=radius)
if (fill) shade3d(ellips, col=surface.col[j], alpha=0.1,
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",
"additive"))
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)),
additive = if (is.null(df.additive)) gam(y ~ s(x) +
s(z))
else 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.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),
additive = if (is.null(df.additive)) gam(y ~ s(x) +
s(z) + groups)
else 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.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,
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]
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),
additive = if (is.null(df.additive)) gam(y ~
s(x) + s(z),
subset=select.obs)
else 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.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,
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)
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
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox
--------------------------------
> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Naiara S. Pinto
> Sent: Tuesday, October 04, 2005 6:13 PM
> To: r-help at 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 at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
More information about the R-help
mailing list