I have a set of linear regression problems, whose dimensionalities vary from 2D (line fitted to a single independent variable), to 3D (plane fitted to two independent variables) to nD where hyperplane is fitted to three or more independent variables. I am interested in obtaining the numerical values for the projections of the constituent datapoints to the best fitting plane (in the case of 3D problem) along all three principal axes (not orthogonal or shortest distance projection). The following page https://rpubs.com/ericroh/334817 gives an idea of what I'd like to achieve, but I'd like to calculate the projections not only along the y axis but along the x1 and x2 axes as well. My data and attempts are well represented by the code behind the link above and I'll reproduce it here
library(scatterplot3d)
dataset = cbind.data.frame(x1 = c(1.9,0.8,1.1,0.1,-0.1,4.4,4.6,1.6,5.5,3.4)
,x2 = c(66, 62, 64, 61, 63, 70, 68, 62, 68, 66)
,y = c(0.7,-1.0,-0.2,-1.2,-0.1,3.4,0.0,0.8,3.7,2.0))
scatterplot3d(x1,x2,y)
ls <- function(dataset, par) { with(dataset, sum((y-par[1]-par[2]*x1-par[3]*x2)^2)) }
result <- optim(par=c(0,0,0), ls, data=dataset)
plot3d <- scatterplot3d(x1, x2, y, angle = 55, pch = 16, color ="red")
my.lm <- lm(y ~ x1 + x2, data = dataset)
plot3d$plane3d(my.lm, lty.box = "solid")
Since the dimensionality of my regression varies, I'd be interested in a general solution that works for any nD problem. All my regressions have the form y ~ x1 + x2 + x3 + ... (no cross-terms) and I am using R's lm function for the regression. Is there R (base if possible) solution available?
I know I can solve a single unknown from equation such as y = A + B * x1 + C * x2 by plugging in the datapoint coordinates (all but the one of interest) but it seems a bit clunky.
I guess you can use Moore-Penrose generalized inverse of matrix, i.e., MASS::ginv, like below
A <- as.matrix(cbind(1, dataset[-length(dataset)]))
b <- dataset[[length(dataset)]]
p <- MASS::ginv(A) %*% b
and you will obtain
> p
[,1]
[1,] -11.4527879
[2,] 0.4502753
[3,] 0.1725176
which coincides with your approaches
> (result <- optim(par = c(0, 0, 0), f, data = dataset))
$par
[1] -11.4591902 0.4500494 0.1726229
$value
[1] 8.8654
$counts
function gradient
270 NA
$convergence
[1] 0
$message
NULL
> (my.lm <- lm(y ~ x1 + x2, data = dataset))
Call:
lm(formula = y ~ x1 + x2, data = dataset)
Coefficients:
(Intercept) x1 x2
-11.4528 0.4503 0.1725