I am trying to plot a plane in a 3D environment in the context of a multivariate regression. I have a linear regression polynomial in the form of f(x)=beta_1 + beta_2*x_1 + beta_3 * x_3.
Obviously this formula describes a plane in R^3. I would like to plot this, but I fail to see how to do this efficiently. Obviously, what I need to do is to create a sort of grid on which I compute the values of my regression polynomial.
So far, this is what I have:
beta <- c(1, 1, 1) # Placeholder. This is the output of my regression analysis
p <- function(k) { t(beta) %*% k }
n <- 20
m <- 50
x <- seq(from=5, to=13, length.out=n)
z <- seq(from=20, to=50, length.out=m)
M <- mesh(x, z)
Now as far as I understand it, whether I'm using plot3D or rgl, I need to compute the value of p on these grid elements and store it in a matrix of dimension m * n. Obviously, in other languages I would just iterate over both matrices and fill the values in y by hand, but this seems very ugly and unfitting for R to me. So after some research, I stumbled upon the function mapply, which applies a function over two matrices. Note that the input to p requires a 1 to be appended to the argument before the call.
mapply(function (x1, x2) { p(c(1, x1, x2)) }, M$x, M$y)
But the output for this is very, very, very ugly and far from what I expected. Could someone help me out? It surprises me that plotting a plane in 3D seems to be so difficult, because after all, all that we need to uniquely determine a plane are three points in space. And yet all frameworks require me to use fairly complicated function calls?
I am sorry if the solution to my issue is trivial. I've checked out quite a few similar threads, but I still can't quite wrap my head around what I need to do as the issue at hand should be so simple.
In rgl it's much simpler: just use planes3d. That function defines planes using the parameterization a x + b y + c z + d = 0, so if we assume x1, f(x) and x3 in your notation are x, y and z in rgl notation, you would plot the plane using
planes3d(beta2, -1, beta3, beta1)
One complication is that planes are infinite, so you have to specify which part you want to plot in some other way. The usual way is to use plot3d() to set up a coordinate system, e.g.
plot3d(x1, y, x3, xlab = "x_1", ylab = "f(x)", zlab = "x_3")