R中含双因子的非线性最小二乘回归拟合报错求助

问题根源 #

你遇到的singular gradient matrix错误主要来自两个核心问题：

• 模型设定逻辑错误：直接在米氏方程（y ~ m*x/(x+k)）后叠加虚拟变量的线性项，不符合因子对非线性模型参数的影响逻辑——因子应该改变米氏方程的参数m和k，而非给整体预测值加线性偏移。
• 共线性问题：四水平因子设置4个虚拟变量（f2c1-f2c4）会导致完全共线性（四个变量的和恒为1），加上两水平因子的虚拟变量后，进一步加剧矩阵奇异性，导致梯度无法计算。

正确的模型构建方式 #

针对你分析因子组合对m和k影响的需求，推荐两种合理的模型设定：

方式1：每个因子组合对应独立的m和k #

假设不同因子组合下，米氏方程的参数完全独立，适合分析每个组合的参数差异：

y ~ m[grp] * x / (x + k[grp])

其中grp是两个因子交叉生成的组合变量（共2*4=8个组合）。

方式2：因子对m和k产生线性效应 #

假设因子对m和k的影响是线性的，适合量化单个因子的主效应：

y ~ (m0 + m1*f1 + m2*f2c1 + m3*f2c2 + m4*f2c3) * x / (x + (k0 + k1*f1 + k2*f2c1 + k3*f2c2 + k4*f2c3))

注意：四水平因子使用3个虚拟变量（treatment编码，以某一水平为基准），避免共线性。

可运行代码示例 #

以下用模拟测试数据演示minpack.lm::nlsLM的正确拟合流程：

步骤1：生成测试数据 #

library(minpack.lm)
library(dplyr)

# 生成因子组合
set.seed(123)
f1 <- factor(rep(c("A", "B"), each = 40))  # 两水平因子，A为基准
f2 <- factor(rep(c("C1", "C2", "C3", "C4"), each = 10, times = 2))  # 四水平因子
x <- rep(seq(1, 20, length.out = 10), times = 8)

# 设定真实参数：f1和f2影响m和k
true_params <- expand.grid(f1 = c("A", "B"), f2 = c("C1", "C2", "C3", "C4")) %>%
  mutate(m = case_when(
    f1 == "A" & f2 == "C1" ~ 10,
    f1 == "A" & f2 == "C2" ~ 12,
    f1 == "A" & f2 == "C3" ~ 15,
    f1 == "A" & f2 == "C4" ~ 8,
    f1 == "B" & f2 == "C1" ~ 14,
    f1 == "B" & f2 == "C2" ~ 16,
    f1 == "B" & f2 == "C3" ~ 18,
    f1 == "B" & f2 == "C4" ~ 11
  ),
  k = case_when(
    f1 == "A" & f2 == "C1" ~ 5,
    f1 == "A" & f2 == "C2" ~ 6,
    f1 == "A" & f2 == "C3" ~ 8,
    f1 == "A" & f2 == "C4" ~ 4,
    f1 == "B" & f2 == "C1" ~ 7,
    f1 == "B" & f2 == "C2" ~ 9,
    f1 == "B" & f2 == "C3" ~ 10,
    f1 == "B" & f2 == "C4" ~ 6
  ))

# 合并数据并生成y值（加噪声）
data <- expand.grid(f1 = f1, f2 = f2, x = x) %>%
  left_join(true_params, by = c("f1", "f2")) %>%
  mutate(y = m*x/(x + k) + rnorm(nrow(.), 0, 0.5))

步骤2：拟合方式1的模型（独立参数） #

# 生成组合变量
data$grp <- interaction(data$f1, data$f2, drop = TRUE)
grp_levels <- levels(data$grp)

# 设置初始值：每个组合的m初始为10，k初始为5
start_vals <- c(
  setNames(rep(10, length(grp_levels)), paste0("m_", grp_levels)),
  setNames(rep(5, length(grp_levels)), paste0("k_", grp_levels))
)

# 拟合模型
model1 <- nlsLM(y ~ get(paste0("m_", grp)) * x / (x + get(paste0("k_", grp))),
                data = data,
                start = start_vals)

summary(model1)

步骤3：拟合方式2的模型（线性效应） #

# 对因子进行treatment编码（默认）
data$f1 <- relevel(data$f1, ref = "A")
data$f2 <- relevel(data$f2, ref = "C1")

# 设置初始值：基准水平的m0=10，k0=5，其余效应初始为0
start_vals2 <- c(m0 = 10, m1 = 0, m2 = 0, m3 = 0, m4 = 0,
                 k0 = 5, k1 = 0, k2 = 0, k3 = 0, k4 = 0)

# 拟合模型
model2 <- nlsLM(y ~ (m0 + m1*I(f1=="B") + m2*I(f2=="C2") + m3*I(f2=="C3") + m4*I(f2=="C4")) * x / 
                  (x + (k0 + k1*I(f1=="B") + k2*I(f2=="C2") + k3*I(f2=="C3") + k4*I(f2=="C4"))),
                data = data,
                start = start_vals2)

summary(model2)

分析因子对m和k的影响 #

• 对于方式1：直接提取每个组合的m和k估计值，用箱线图或方差分析比较不同因子水平下的参数差异。
• 对于方式2：查看summary(model2)中的系数显著性，判断单个因子水平对m和k的线性影响是否显著。

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