沿各坐标轴投影数据点至回归超平面的通用R实现问询

我有一系列维度不同的线性回归任务，涵盖：

• 2D场景：拟合单自变量的直线
• 3D场景：拟合双自变量的平面
• nD场景：拟合三个及以上自变量的超平面

我需要获取每个数据点沿所有主轴（x₁、x₂、…、y）到最优拟合超平面的投影数值（注意：不是正交/最短距离投影）。

目前能找到的实现思路仅支持沿y轴计算投影，但我需要同时处理所有坐标轴。我的测试数据和初步尝试代码如下：

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")

由于回归维度可变，我需要一个适用于任意nD场景的通用方案，所有回归都采用y ~ x1 + x2 + x3 + ...的形式（无交叉项），且使用R的lm函数执行回归。是否有基于R基础包的解决方案？我知道可以通过固定其他坐标求解单个未知量的方式计算，但这种方法太繁琐。

核心思路 #

线性回归的超平面方程可表示为：
$$\hat{y} = \beta_0 + \beta_1x_1 + \beta_2x_2 + ... + \beta_kx_k$$
沿某一主轴的投影，本质是固定该轴外的所有变量值，代入超平面方程求解该轴的投影值。我们可以通过提取lm模型的系数，编写通用函数自动处理所有坐标轴的投影计算。

实现代码 #

# 通用投影计算函数
compute_axis_projections <- function(lm_model, data) {
  # 提取模型系数、自变量名和因变量名
  coefs <- coef(lm_model)
  var_names <- names(coefs)[-1]  # 剔除截距项对应的名称
  response_name <- as.character(lm_model$terms[[2]])
  
  # 初始化结果数据框，保留原始数据
  proj_results <- data
  
  # 计算y轴投影（即模型拟合值）
  proj_results[[paste0("proj_", response_name)]] <- fitted(lm_model)
  
  # 计算每个x轴的投影
  for (var in var_names) {
    proj_col <- paste0("proj_", var)
    # 变形超平面方程：x_var = (y - β0 - Σ(β_j x_j))/β_var （j≠var）
    numerator <- data[[response_name]] - coefs[1]
    # 减去其他自变量的β*x项
    for (other_var in setdiff(var_names, var)) {
      numerator <- numerator - coefs[other_var] * data[[other_var]]
    }
    # 计算当前x轴的投影值
    proj_results[[proj_col]] <- numerator / coefs[var]
  }
  
  return(proj_results)
}

# 测试3D场景（用户提供的数据集）
my.lm <- lm(y ~ x1 + x2, data = dataset)
proj_data_3d <- compute_axis_projections(my.lm, dataset)
print(proj_data_3d)

# 测试2D场景
set.seed(123)
data_2d <- data.frame(x1 = rnorm(10), y = 2 + 3*x1 + rnorm(10))
lm_2d <- lm(y ~ x1, data = data_2d)
proj_data_2d <- compute_axis_projections(lm_2d, data_2d)
print(proj_data_2d)

代码说明 #

1. 自动适配维度：函数会从lm模型中自动识别自变量数量，无需手动修改即可处理2D、3D或更高维度的回归场景。
2. y轴投影：直接调用fitted()获取模型拟合值，这就是数据点沿y轴到拟合面的投影。
3. x轴投影：对每个自变量，将超平面方程变形为以该变量为求解目标的形式，代入数据点的实际y值和其他自变量值，计算得到该轴的投影值。
4. 结果输出：返回的结果数据框包含原始数据和所有轴的投影值，方便后续分析。

内容的提问来源于stack exchange，提问作者Suola