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

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@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
> PLEASE do read the posting guide! 
> http://www.R-project.org/posting-guide.html

______________________________________________
R-help@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 Received on Wed Oct 05 11:29:23 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:40:35 EST