R中线性方程组重复测量数据最小二乘参数估计与标准误计算

问题定义 #

待估计的两类方程结构如下，所有$K_j$、$K_{ij}$均带有至少10次重复观测值，需要最小化残差平方和估计未知参数$X_i$、$X_{ij}$、$x_{ij}$，同时输出参数估计值的标准误：

• 形式1（线性加性模型）：$X_i + K_j + X_{ij} = K_{ij}$
• 形式2（带交互项的非线性模型）：$X_i + K_j + x_{ij} \times X_i \times K_j = K_{ij}$

实现逻辑 #

两类模型都可以通过R的最小二乘拟合函数直接实现，不需要手动编写优化与标准误计算逻辑：

1. 形式1是完全线性结构，将已知观测值$K_j$作为固定偏移项，用普通最小二乘lm()拟合即可
2. 形式2存在参数乘积项，属于非线性最小二乘问题，用收敛性更好的minpack.lm::nlsLM()拟合即可
   模型输出会直接返回参数估计值与基于残差的标准误，结果满足最小化残差平方和的要求。

完整实现代码 #

1. 示例数据预处理 #

先对原始构造代码做变量名匹配、因子转换，方便模型识别分组：

set.seed(123) # 设定随机种子保证结果可复现
i <- rep(1:3, each = 30)
j <- rep(rep(1:3, each = 10), 3)
K_j  <- rep(c(6, 5, 10), each = 30) + rnorm(90) # 匹配公式中的K_j
K_ij <- K_j + rnorm(90)

# 转换为因子，生成分组交互项
i <- as.factor(i)
j <- as.factor(j)
ij <- interaction(i, j) # 对应i-j交叉分组，匹配X_ij、x_ij的分组
data <- data.frame(i, j, ij, K_j, K_ij)

2. 形式1（线性模型）拟合 #

因为$K_j$是已知观测值、系数固定为1，所以作为偏移项（offset）加入模型，不参与系数估计：

# 模型拟合：去掉截距项，估计i组主效应、ij组交互效应
fit1 <- lm(K_ij ~ 0 + i + ij + offset(K_j), data = data)
# 提取参数估计值和标准误
result1 <- summary(fit1)$coefficients[, c("Estimate", "Std. Error")]
print(result1)

结果解释：

• 以i开头的行对应不同组的$X_i$估计值和标准误
• 以ij开头的行对应不同i-j交叉组的$X_{ij}$估计值和标准误
  示例数据中参数真实值为0，估计结果会在0附近波动。

3. 形式2（非线性模型）拟合 #

使用非线性最小二乘拟合，以形式1的估计结果作为参数初值提升收敛速度：

# 安装依赖包（首次运行执行）
# install.packages("minpack.lm")
library(minpack.lm)

# 提取线性模型结果作为初值
start_Xi <- coef(fit1)[paste0("i", levels(i))]
start_xij <- rep(0, length(levels(ij))) # x_ij初值设为0

# 拟合非线性模型
fit2 <- nlsLM(
  formula = K_ij ~ X_i[i] + K_j + x_ij[ij] * X_i[i] * K_j,
  data = data,
  start = list(
    X_i = setNames(start_Xi, levels(i)),
    x_ij = setNames(start_xij, levels(ij))
  )
)

# 提取参数估计值和标准误
result2 <- summary(fit2)$coefficients[, c("Estimate", "Std. Error")]
print(result2)

结果解释：

• X_i开头的行对应不同组的$X_i$估计值和标准误
• x_ij开头的行对应不同i-j交叉组的$x_{ij}$估计值和标准误

注意事项 #

• 重复测量样本量足够支撑标准误估计：示例数据共90个观测值，形式1待估参数共12个（3个$X_i$+9个$X_{ij}$），形式2待估参数共21个（多9个$x_{ij}$），残差自由度充足。
• 如果数据存在异方差，可以在拟合完成后用sandwich包计算稳健标准误，替换模型默认输出的标准误即可。
• 大规模数据下lm()运算效率完全满足要求，nlsLM的Levenberg-Marquardt算法收敛速度远优于基础nls，适合大样本场景。

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