2

I want to plot the regression surface from a model with an interaction term using rgl's interactive plotting system. It is easy to plot a regression plane for a model without an interaction term using: 

plot3d(x=x1, y=x2, z=y1, type="s", col="yellow", size=1)
planes3d(a=coef(mod1)[2], b=coef(mod1)[3], c=-1, d=coef(mod1)[1], alpha=.5)

However, when the plane twists, this seems to be more difficult. Following on this question: 3D equivalent of the curve function in r, I am trying: 

f2 <- function(x, y) as.vector(coef(mod2)%*%c(1, x, y, x*y))

curve_3d <- function(f2, x_range=c(0, 40), y_range=c(0, 40)){ 
  if (!require(rgl) ) {stop("load rgl")}

  xvec <- seq(x_range[1], x_range[2], by=1)
  yvec <- seq(y_range[1], y_range[2], by=1)
  fz   <- outer(xvec, yvec, FUN=f2)
  persp3d(xvec, yvec, fz, alpha=.5)
}
open3d()
plot3d(x=x1, y=x2, z=y2, type="s", col="yellow", size=1)
curve_3d(f2)

But, it's not working. (I've tried some other things as well, but I'm keeping this short.) My main problem so far seems to be with f2; however, I will also want this to look like planes3d, and I'm not sure if this is going to give me a wireframe.

Here's an example: 

set.seed(897)
x1 = rep(c(0, 10, 20, 30, 40), times=25)
x2 = rep(c(0, 10, 20, 30, 40), each=25)
y2 = 37 + 0.7*x1 + 1.2*x2 - 0.05*x1*x2 + rnorm(125, mean=0, sd=5)
mod2 = lm(y2~x1*x2)
open3d()
plot3d(x=x1, y=x2, z=y2, type="s", col="yellow", size=1)
curve_3d(f2)
Community
  • 1
  • 1
gung - Reinstate Monica
  • 11,583
  • 7
  • 60
  • 79
  • The usual way is to use `predict` with a model, rather than trying to calculate the prediction from a coefficient vector. You are not saying what errors you are getting, but it looks as though your arguments are not the same as the ones you are passing. `curve3d` is supposed to get three arguments and you are only passing one. – IRTFM Feb 18 '13 at 04:55
  • I'd be happy with that. How do you get the predictions from the model? I tried using `z=mod2$fitted.values`, but it also didn't work. (I've tried several things.) – gung - Reinstate Monica Feb 18 '13 at 04:58
  • The specific error that I've been able to track down is "Error in coef(mod2) %*% c(1, x, y, x * y) : non-conformable arguments". – gung - Reinstate Monica Feb 18 '13 at 04:59
  • For worked illustrations, I expect a data object. – IRTFM Feb 18 '13 at 04:59
  • I don't understand your comment. – gung - Reinstate Monica Feb 18 '13 at 05:00
  • http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example – IRTFM Feb 18 '13 at 05:11
  • Better solution: http://stackoverflow.com/questions/18147595/plot-3d-plane-true-regression-surface – Tal Galili Jun 05 '16 at 15:42

1 Answers1

4
grd <- expand.grid(x1=unique(x1), x2=unique(x2) )
grd$pred <-predict(mod2, newdata=grd)
persp3d(x=unique(grd[[1]]), y=unique(grd[[2]]), 
              z=matrix(grd[[3]],5,5), add=TRUE)

enter image description here

IRTFM
  • 258,963
  • 21
  • 364
  • 487