[R] Help on plotting a regression plane?

[Uwe Ligges]

Peter Dalgaard BSA wrote:
> 
> 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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._