如何为非线性混合效应(nlme)模型添加随机效应？报错求解

问题背景 #

需要构建非线性混合效应模型描述x与y的关系，模型公式为：y = theta*(1-exp(-beta*x))，其中y随分组变量r随机变化。已通过nls()实现无随机效应的模型，但使用nlme()加入随机效应时，出现错误：

Error in eval(predvars, data, env) : object 'theta' not found

数据集构建代码 #

x = c(33,35,16,8,31,31,31,23,7,7,7,7,11,11,3,3,6,6,32,32,1,17,17,17,25,40,40,6,6,29,29,13,23,23,44,44,43,43,13,4,6,15,17,22,28,8,11,22,32,6,12,20,27,15,29,29,29,29,29,12,12,16,16,12,12,2,49,49,14,14,14,37,2.87,4.86,7.90,11.95,16.90,16.90,18.90,18.89,22.00,24.08,27.14,30.25,31.22,32.26,7,14,19,31,36,7,14,19,31,36,7,16,16,16,16,16,16,32,32,32,32,32,32,11,11,11,13,13,13,13,13,13,13,13,13,13,13,13,9,9)
y = c(39.61,32.66,27.06,21.74,22.18,38.19,35.02,23.13,9.70,14.20,13.40,15.30,18.80,19.00,3.80,4.90,15.00,14.20,24.90,16.56,1.76,29.29,28.49,18.64,27.10,9.47,14.14,10.27,8.44,26.15,25.43,22.00,19.00,13.00,73.19,67.76,32.34,36.86,8.00,1.57,8.33,16.20,14.69,18.95,20.52,4.92,8.28,15.27,18.37,6.60,10.98,12.56,19.04,5.49,21.00,12.90,17.30,11.40,12.20,15.63,15.22,33.80,17.78,19.33,3.86,8.57,30.40,13.39,11.93,4.55,6.18,12.70,2.71,7.23,5.61,22.74,15.71,16.95,18.31,20.78,17.64,20.00,19.52,24.86,30.06,24.92,4.17,11.02,10.08,14.94,25.98,0.00,3.67,3.67,6.69,11.90,5.06,13.21,10.33,0.00,0.00,6.47,8.38,28.57,25.26,28.67,27.92,33.69,29.61,6.11,7.13,6.93,4.81,15.34,4.90,14.94,8.88,10.24,8.80,10.46,10.48,9.19,9.67,9.40,24.98,50.79)
r = c("A","A","A","A","A","A","A","A","B","B","B","B","B","B","B","B","B","B","C","C","D","E","E","E","F","G","G","H","H","H","H","I","I","I","J","J","J","J","K","L","L","L","L","L","L","L","L","L","L","L","L","L","L","M","N","N","N","N","N","O","P","P","P","P","P","Q","R","R","S","S","S","T","U","U","U","U","U","U","U","U","U","U","U","U","U","U","V","V","V","V","V","V","V","V","V","V","W","X","X","X","X","Y","Y","Z","Z","Z","Z","Z","Z","AA","AA","AA","AB","AB","AB","AB","AB","AB","AB","AB","AB","AB","AB","AB","AC","AC")
df = data.frame(x,y,r)

无随机效应的nls模型代码 #

nls_test = nls(y~theta*(1-exp(-beta*x)),
              data = df,
              start = list(beta = 0.2, theta = 38), 
              trace = TRUE)

报错的nlme模型代码 #

nlme_test = nlme(y ~ theta*(1-exp(-beta*x)),
                 fixed = x~1,
                 random = r~1,
                 data = df,
                 start = c(theta = 38,
                           beta = 0.2))

错误原因 #

1. fixed参数定义错误：nlme()的fixed参数需要指定模型中固定效应的参数（theta和beta），而非将x作为响应变量。原代码中fixed = x~1完全不符合要求，导致模型无法识别theta和beta。
2. random参数格式错误：random参数需要明确指定哪个参数随分组r变化，格式应为参数名 ~ 1 | 分组变量，原代码random = r~1没有指定随机效应对应的模型参数。

修正方案 #

假设需求是让theta随分组r产生随机变化，beta为全局固定效应，修正后的代码如下：

nlme_test = nlme(y ~ theta*(1-exp(-beta*x)),
                 fixed = list(theta ~ 1, beta ~ 1),  # 声明theta和beta为固定效应
                 random = theta ~ 1 | r,  # theta随分组r随机变化
                 data = df,
                 start = c(theta = 38, beta = 0.2))  # 固定效应的初始值

扩展说明 #

• 如果beta也需要随分组随机变化，可将random参数调整为：

  random = list(theta ~ 1 | r, beta ~ 1 | r)
  

• nls()能正常运行是因为它不需要区分固定/随机效应，直接拟合全局参数；而nlme()对参数结构有严格要求，必须明确界定固定和随机效应的对应关系。

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