Mathematics behind how the normal to the plane is derived from a linear model when plotting with planes3d

Question

I've been experimenting with the rgl library to draw 3d plots of some linear models. The package works really nicely, but I have a question on the math behind how one draws a plane (using planes3d) that shows the 'fit' (predicted values ) of the linear model, for a model with 2 parameters.

The R package documentation says that the arguments -- a, b, c, and d -- of a call to planes3d must be chosen so they are the coordinates of the normal to the plane (given by the model). And there is an example on page 33 of this document: http://cran.r-project.org/web/packages/rgl/rgl.pdf, which I have slightly modified to use the data from the mtcars package (see R code below).

In the listing I create the model from the mtcars data at line 3, then I plot the data points at line 6. Next I plot the plane in lines 9 - 15. The a,b,c,d coordinates are chosen per instructions in the documentation. In particular,'c' is set to be -1, which is currently just a 'magic' value to me. I don't understand why -1 is chosen for 'c'.

But it works out. In line 18 - 26 I define vectors for the normal to the plane, and for a point lying within the plane. I take the dot product of these two and add the intercept (offset) and I get zero. It all works out as expected. But I am hoping someone can explain the theory behind >why< it works out.

R code:

 1  library(rgl)
 2  attach (mtcars)
 3  model <- lm(mpg ~  cyl + disp)
 4  
 5  # Plot data points 
 6  plot3d(cyl, disp, mpg, type="s", col="red", size=1) 
 7  
 8  
 9  # Plot the plane of points that make up the predicted values of the model
10  coefs <- coef(model)
11  a <- coefs["cyl"]
12  b <- coefs["disp"]
13  c <- -1
14  d <- coefs["(Intercept)"]
15  planes3d(a, b, c, d, alpha=0.5)
16  
17  
18  # define a normal vector to the plane
19  #
20  normal = c(a,b,c)
21  
22  # define a vector within the plane
23  #
24  newdata <- data.frame(cyl = 5, disp = 200)
25  predicted = predict(model,  newdata)
26  pointOnPlane = c(5, 200, predicted)
27  
28  
29  
30  # Verify that the dot products of the normal and the vector that defines pointOnPlane  is 0 
31  # (after we subtract the offset 'd', which came from  the intercept  of the model). 
32  
33  # This will print as 'TRUE'
34  sum(pointOnPlane * normal)  + d   == 0

epilogue: Got a great answer from Andre Silva, below. Thanks, Andre ! makes sense now ;^)

Answer 1

Let me give it a try.

What is the main idea here? Well, as you said, we want to show the plane that intersects the predictions of a linear model. So, you have this linear model model <- lm(mpg ~ cyl + disp) which means this:

$$z = \alpha x + \beta y + \text{intercept}$$

where $z = \text{mpg}$, $x = \text{cyl}$, $y = \text{disp}$

When you obtain the coefficients for this regression, what you're getting is a machine that predicts values for $z$ given some $x$ and $y$. Not only that, you get for free the equation of a plane:

$$\alpha x + \beta y + \text{intercept} - z = 0$$

Therefore, when you plot this plane using planes3d(a, b, c, d, alpha=0.5), you are effectively saying "Plot a point for every $x$, $y$ and $z$ that satisfies this equation." In this example, a and b is the $\alpha$ and $\beta$ you got from the regression. $z$ is actually given in the plane equation. That's why you need to set c to $-1$. d is not necessary in order to define a plane, but since it was added to the function planes3dand the regression depends on the intercept, you need to use it as d. Otherwise, the plane would be d units above (or below depending on the sign) the points it's supposed to fit.

The normal vector for the plane $\alpha x + \beta y + \text{intercept} - z = 0$ is $(\alpha, \beta, -1)$, so that is what they are doing with normal = c(a,b,c).

After this, they predict a point $z$ given some coordinates $x$ and $y$ and construct a vector on the plane:

newdata <- data.frame(cyl = 5, disp = 200)
predicted = predict(model,  newdata)
pointOnPlane = c(5, 200, predicted)

Since the dot product of a vector that lies on the plane with a normal vector to the plane is equal to 0 ($\vec{n}\cdot \vec{r} = 0$), then you're verifying that pointOnPlane is actually on the plane, which is equivalent to this part:

sum(pointOnPlane * normal)  + d   == 0 

And I think that is pretty much all there is to it.