使用nls.lm无法获取95%置信区间：奇异Hessian矩阵问题

问题概述 #

使用minpack.lm::nls.lm拟合分段非线性模型后，显著参数a和c的95%置信区间计算失败，报错提示Hessian矩阵计算奇异（system is computationally singular）。从拟合结果可见，参数b、d、f、g的标准误差极大，p值为1，说明这些参数无法被现有数据有效估计，导致Hessian矩阵不可逆。

可行解决方案 #

1. 简化模型，移除无信息参数 #

由于b、d、f、g对模型解释力无贡献，直接移除这些参数，保留可被数据有效估计的a、c、e：

# 简化后的模型函数
custom_model_simplified <- function(day_of_year, params){
  params$c * exp(-1/2 * (((day_of_year - params$a)/params$e)^2))
}

# 残差函数
model_residuals_simplified <- function(params){
  observed_data$y - custom_model_simplified(observed_data$x, params)
}

# 重新拟合简化模型
model_fit_simplified  <- minpack.lm::nls.lm(
  par = list(a = 148, c = 0.02, e = 1), 
  fn = model_residuals_simplified, 
  lower = c(0, 0, 0)
)

# 查看拟合结果
summary(model_fit_simplified)
# 计算置信区间
confint(model_fit_simplified)

简化后模型的参数均能被数据约束，Hessian矩阵不再奇异，confint()可正常输出a和c的置信区间。

2. 采用参数剖面法计算置信区间 #

若需保留原模型结构，可使用参数剖面似然法绕开Hessian矩阵求逆步骤，这是非线性模型置信区间计算的稳健方法：

# 生成参数剖面
model_profile <- minpack.lm::profile(model_fit)
# 提取a和c的95%置信区间
confint(model_profile, parm = c("a", "c"))

参数剖面通过固定目标参数，调整其他参数使似然函数下降到临界值，无需依赖可逆的Hessian矩阵，适用于存在参数冗余的场景。

3. 正则化Hessian矩阵（备选方案） #

对奇异的Hessian矩阵添加微小对角扰动（岭回归思想），使其可逆后再计算置信区间：

# 提取模型的Hessian矩阵和参数估计值
hess <- model_fit$hessian
par_est <- model_fit$par

# 添加小的对角扰动（可根据情况调整大小）
ridge <- 1e-6 * diag(nrow(hess))
hess_reg <- hess + ridge

# 计算正则化后的协方差矩阵
cov_mat <- solve(hess_reg)

# 基于正态近似计算95%置信区间
se <- sqrt(diag(cov_mat))
ci <- cbind(
  lower = par_est - 1.96 * se,
  upper = par_est + 1.96 * se
)

# 提取a和c的置信区间
ci[c("a", "c"), ]

注意：该方法的结果依赖于扰动值的选择，需谨慎验证结果合理性，仅作为其他方法不可行时的备选。

内容的提问来源于stack exchange，提问作者Sophie Père