how do i plot a linear model lm using scatterplot3d via the plane3d function?

scatterplot3d() will not be able to plot models with larger dimensionality (than 2 input dimensions and 1 output dimension) in 3D. In fact, such a plot would not be valid since the values in the additional dimensions will presumably be different for the different observations. They would therefore influence how closely the model fits and a plot that neglects these would be misleading.

That said, s3d$plane3d does not handle malformed input very well. For instance, if the dimensionality of the model is not as expected, it will return confusing error messages (as you have seen). There is also no help for this function and in fact the function is nested in another function in the package and has no comments. As a result this will all be fairly difficult to understand, but if you want to go deeper you have to read the code of the package, which you can find here.

You can absolutely have your plot show a partial regression surface, but you have to tell plot3d, which dimensions you want. Essentially you’d be plotting a plane in 3d space where you should have a hyperplane in higher dimensional space.

Your attempt 2 was on the right track. But you do not hand over the right argument. The function wants x.coef and y.coef etc. but not xyz.coords and therefore it apparently tries to interpret the vectors you hand over as coefficient and fails. You could do this instead:

s3d$plane3d(Intercept=lms$coefficients["(Intercept)"][[1]],
            x.coef=lms$coefficients["elem$pH"][[1]],
            y.coef=lms$coefficients["elem$water_Loss"][[1]],
            draw_polygon = TRUE, 
            draw_lines = TRUE, 
            polygon_args = list(col = rgb(0.8, 0.8, 0.8, 0.8)))

However, it is unlikely, however, that you will even see the regression surface in your plot, because the influence of the dimensions you do not plot will shift it out of the visible area of your figure. If you want to pull it back by force, you have to modify your intercept:

average_intercept <- lms$coefficients["(Intercept)"][[1]] + lms$coefficients["elem$elev"][[1]] * mean(elem$elev)
s3d$plane3d(Intercept=average_intercept,
            x.coef=lms$coefficients["elem$pH"][[1]],
            y.coef=lms$coefficients["elem$water_Loss"][[1]],
            draw_polygon = TRUE, 
            draw_lines = TRUE, 
            polygon_args = list(col = rgb(0.8, 0.8, 0.8, 0.8)))

But the plane you see is actually only a 2d-slice through the 3d-surface that is your regression and accurately only represents the observations you that happen to have exactly the average value in that third dimension (elev in your case).

In fact, this is exactly what you would get if you ran the regression without the additional dimension(s); so you might as well do and plot that.