R中非线性混合效应模型构建报错求助：对象未找到

需要在R中构建非线性混合效应模型，公式如下：
h = 1.30 + (dbh²/(b0 + b1*dbh)²)
其中b0由线性组合计算：b0 = a0 + a1*BA + a2*VH + a3*QMD + a4*NH，a0-a4为线性回归系数，且b0包含随分组变化的随机效应。

尝试的代码出现两类报错：

• Error in h(dbh) : argument "b0" is missing, with no default
• Error in eval(predvars, data, env) : object 'b1' not found

附原始数据（dput输出）：

your_data <- structure(list(dbh = c(22L, 87L, 110L, 27L, 69L, 34L, 49L, 89L, 
49L, 87L, 46L, 66L, 74L, 48L, 81L, 26L, 89L, 66L, 100L, 50L, 
74L, 29L, 82L, 58L, 75L, 63L, 100L, 54L, 54L, 64L), h = c(11.4, 
32.1, 31.3, 35.3, 40, 23.6, 34, 35, 30, 43, 36, 41.3, 40, 33.3, 
40.6, 21.4, 39.2, 42, 44, 36.5, 39.6, 24, 29, 39.5, 37, 42, 41.3, 
34.5, 37, 40), BA = c(27.71, 27.71, 43.35, 43.35, 26.48, 26.48, 
34.33, 34.33, 35.29, 35.29, 28.61, 28.61, 39.86, 18.1, 18.1, 
29.8, 29.8, 30.2, 30.2, 47.1, 47.1, 30.69, 30.69, 31.78, 31.78, 
45.49, 45.49, 9.34, 9.34, 34.65), VH = c(393.76, 393.76, 672.33, 
672.33, 388.41, 388.41, 526.63, 526.63, 538.45, 538.45, 396.61, 
396.61, 592.49, 249.89, 249.89, 445.99, 445.99, 466.5, 466.5, 
710.59, 710.59, 447, 447, 477.61, 477.61, 726.05, 726.05, 122.4, 
122.4, 553.31), QMD = c(36.8, 36.8, 53.9, 53.9, 43.3, 43.3, 54, 
54, 51.4, 51.4, 36.1, 36.1, 62.5, 27.2, 27.2, 43.5, 43.5, 46.2, 
46.2, 48, 48, 37.3, 37.3, 47.4, 47.4, 60.2, 60.2, 25.7, 25.7, 
60.6), NH = c(260L, 260L, 190L, 190L, 180L, 180L, 150L, 150L, 
170L, 170L, 280L, 280L, 130L, 310L, 310L, 200L, 200L, 180L, 180L, 
260L, 260L, 280L, 280L, 180L, 180L, 160L, 160L, 180L, 180L, 120L
)), row.names = c(NA, 30L), class = "data.frame")

• 函数名与响应变量冲突：自定义函数命名为h，与数据框中的响应变量h重名，导致调用混乱
• 公式语法错误：nlmer的公式格式不符合要求，未正确表达非线性模型与随机效应的关系
• 变量名不一致：自定义b0函数中使用了QM，但数据中对应的变量是QMD
• 缺失分组变量：混合效应模型需要分组标识，但原始数据中无ID列，无法定义随机效应的分组单元

步骤1：添加分组ID #

观察数据可知，每两行的环境变量（BA、VH、QMD、NH）完全一致，属于同一分组，因此添加ID列：

library(lme4)
# 为数据添加分组ID
your_data$ID <- rep(1:15, each = 2)

步骤2：定义非线性模型函数 #

避免与响应变量重名，将函数命名为height_model：

height_model <- function(dbh, b0, b1) {
  1.30 + (dbh^2 / (b0 + b1 * dbh)^2)
}

步骤3：设置合理初始值 #

非线性模型对初始值敏感，需根据数据范围调整初始参数：

start_vals <- list(
  b1 = 0.1,
  a0 = 10,
  a1 = 0.1,
  a2 = 0.01,
  a3 = 0.1,
  a4 = 0.01
)

步骤4：拟合非线性混合效应模型 #

直接在nlmer公式中表达b0的线性组合，并指定随机效应：

model <- nlmer(
  h ~ height_model(dbh, b0 = a0 + a1*BA + a2*VH + a3*QMD + a4*NH, b1) ~ (a0 | ID),
  data = your_data,
  start = start_vals,
  control = nlmerControl(optimizer = "bobyqa", maxfun = 1e5)
)

说明：

• (a0 | ID)表示a0随分组ID随机变化，满足b0包含随机效应的需求
• 使用bobyqa优化器并增加迭代次数，提升模型收敛性

步骤5：查看模型结果 #

summary(model)
# 绘制残差图评估拟合效果
plot(model)

• 若模型不收敛，可尝试调整初始值（如增大a0、修改b1的范围）
• 若需要多个随机效应（如a1也随ID变化），可修改随机部分为(a0 + a1 | ID)，但需确保数据量支持多参数估计

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