[R] Is there a regression surface demo?

[Joshua Wiley]

Hi All,

Does anyone know of a function to plot a regression surface for two
predictors?  RSiteSearch()s and findFn()s have not turned up what I
was looking for.  I was thinking something along the lines of:
http://mallit.fr.umn.edu/fr5218/reg_refresh/images/fig9.gif

I like the rgl package because showing it from different angles is
nice for demonstrations.  I started to write my own, but it has some
issues (non functioning code start below), and I figured before I
tried to work out the kinks, I would ask for the list's feedback.

Any comments or suggestions (about functions or preferred idioms for
what I tried below, or...) are greatly appreciated.

Josh


RegSurfaceDemo <- function(formula, data, xlim, ylim, zlim,
                           resolution = 100) {
  require(rgl)
  ## This cannot be the proper way to extract variable names from formula
  vars <- rownames(attr(terms(formula), "factors"))

  ## if no limits set, make them nearest integer to
  ## .75 the lowest value and 1.25 the highest
  ranger <- function(x) {
    as.integer(range(x) * c(.75, 1.25))
  }
  if(is.null(xlim)) {xlim <- ranger(data[, vars[2]])}
  if(is.null(ylim)) {ylim <- ranger(data[, vars[3]])}
  if(is.null(zlim)) {zlim <- ranger(data[, vars[1]])}

  ## This does not actually work because the data frame
  ## does not get named properly (actually it throws an error)
  ##f <- function (x, y) {
  ##  predict(my.model, newdata = data.frame(vars[2] = x, vars[3] = y))
  ##}

  ## Fit model
  my.model <- lm(formula = formula, data = data)

  ## Create X, Y, and Z grids
  X <- seq(from = xlim[1], to = xlim[2], length.out = resolution)
  Y <- seq(from = ylim[1], to = ylim[2], length.out = resolution)
  Z <- outer(X, Y, f)

  ## Create 3d scatter plot and add the regression surface
  open3d()
  with(data = data,
       plot3d(x = vars[2], y = vars[3], z = vars[1],
              xlim = xlim, ylim = ylim, zlim = zlim))
  par3d(ignoreExtent = TRUE)
  surface3d(X, Y, Z, col = "blue", alpha = .6)
  par3d(ignoreExtent = FALSE)
  return(summary(my.model))
}



-- 
Joshua Wiley
Ph.D. Student, Health Psychology
University of California, Los Angeles
http://www.joshuawiley.com/