Help on plotting a regression plane?

John Williams <jwilliams at business.otago.ac.nz> writes:

I'm a bit embarrassed to ask this, but can anyone tell me (or point me
to a reference) on how to plot a regression plane in R?

I have a regression model with two independent variables and have
created a nice 3d scatterplot with scatterplot3d (thanks Uwe!) and now
would like to overlay the regression plane (gridded, preferably.)  Ay
pointers would be appreciated.

Doesn't look like something to be embarrassed about... (Uwe, you might
want to add this to one of the examples).

Thank you, Peter. If I have got some more time, I will do (see below).

As far as I can see the key would be to use the points3d() function
returned by scatterplot3d(). Something like:

model <- lm(....)
s <- scatterplot3d(....)

for (x in seq(...)) {
    line.x <- rep(x,2)
    line.y <- c(ymin, ymax)
    line.z <- predict(model, new = data.frame(x=line.x, y=line.y)
    s$points3d(line.x, line.y, line.z, type="l", lty="dotted")
}
for (y in seq(....)) {
    ...
}

Alternatively, you might generate the endpoints for a segments() call
by first producing the coordinates in 3d using predict and then
convert them with s$xyz.convert(). That'd be

d1<-rbind(data.frame(x=xmin,y=seq(...)), data.frame(x=seq(...),y=ymin)
d2<-rbind(data.frame(x=xmax,y=seq(...)), data.frame(x=seq(...),y=ymax)
p1 <- predict(model, new=d1)
p2 <- predict(model, new=d2)

line.begin <- s$xyz.convert(d1$x,d1$y,p1)
line.end <- s$xyz.convert(d2$x,d2$y,p2)

segments(line.begin$x,line.begin$y,line.end$x,line.end$y, lty="dotted")

[As you may have gathered, none of this is actually tested...]

Both ideas are very nice. I think predict() has to be used slightly
different.

Nevertheless, in the morning I begun writing a (now tested) example,
also using xyz.convert(). Here is the modified example 5 from the help:

 library(scatterplot3d)
 data(trees)
 attach(trees)
 s3d <- scatterplot3d(trees, type = "h", highlight.3d = TRUE,
     pch=16, main="scatterplot3d -- example 5", angle=60, scale.y=0.7)
 # I want to draw this regression plane as a grid:
 my.lm <- lm(Volume ~ Girth + Height)   
 my.pred <- function(x, y, model){
    my.coef <- coef(model)
    return(my.coef[1] + x * my.coef[2] + y * my.coef[3])
 }
 
 # Now drawing segments from (x1, y1, z1) to (x2, y2, z2). 
 # Looking at the existing plot helps to find the points at first:
 
 x1 <- seq(8, 22, 2)  # Girth
 y1 <- rep(60, 8)     # Height
 z1 <- my.pred(x1, y1, my.lm)
 x2 <- seq(8, 22, 2)
 y2 <- rep(90, 8)
 z2 <- my.pred(x2, y2, my.lm)

 # Now calculating the 3d -> 2d projection for the existing plot
 seg1 <- s3d$xyz.convert(x1, y1, z1)
 seg2 <- s3d$xyz.convert(x2, y2, z2)
 segments(seg1$x, seg1$y, seg2$x, seg2$y, lty="dashed")

 # As above, but the orthogonal lines of the grid:
 x1 <- rep(8, 7)        # Girth
 y1 <- seq(60, 90, 5)   # Height
 z1 <- my.pred(x1, y1, my.lm)
 x2 <- rep(22, 7)
 y2 <- seq(60, 90, 5)
 z2 <- my.pred(x2, y2, my.lm)
 seg1 <- s3d$xyz.convert(x1, y1, z1)
 seg2 <- s3d$xyz.convert(x2, y2, z2)
 segments(seg1$x, seg1$y, seg2$x, seg2$y, lty="dashed")



That's a workaround if you only want to plot it once, but surely
automating it would be nicer. One can also think about drawing
non-linear planes. Maybe this will become a additional function for the
package.

If there a any questions remaining, feel free to send a private mail.

Uwe Ligges
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._