能否用nlme包拟合带参数约束的混合效应模型？

问题背景 #

我需要拟合带有参数上下限约束的混合效应模型，由于后续要做联合建模（Joint Modeling），必须依赖JM包对应的nlme工具。nlme的lme()函数可以通过lmeControl()指定优化方法，其中optim()的L-BFGS-B算法支持设置参数上下限，但lmeControl()只能指定"optim"或"nlminb"这类字符串，没法直接传入自定义优化函数。我尝试给lme()和lmeControl()加upper/lower参数完全无效，自己写自定义优化函数传入也没成功（可能函数写法有问题），想知道有没有办法用nlme实现带参数约束的混合效应模型。

尝试过的代码 #

require(JM)
require(nlme)
require(tidyverse)

data(aids)
data(aids.id)

# 拟合无约束混合效应模型
fitLME <- lme(sqrt(CD4) ~ obstime + obstime:drug, random = ~ obstime | patient, data = aids)

fitLME %>% summary()

# 创建尝试加入约束的控制对象
ctrl <- lmeControl(optimMethod = "L-BFGS-B", opt = "optim",  lower = c(0,0, 0), upper = c(1,20,1))

# 用带约束的控制对象拟合模型
fitLME_ctrl <- lme(sqrt(CD4) ~ obstime + obstime:drug, random = ~ obstime | patient, data = aids, control = ctrl)

# 结果和无约束模型完全一致，约束未生效
fitLME_ctrl %>% summary()

可行解决方案 #

核心思路：通过局部环境替换optim函数传递约束 #

lmeControl()的opt参数仅接受函数名称字符串，但我们可以在局部环境中定义同名的optim函数，把约束参数嵌入其中，让lme()调用这个自定义版本的优化函数。

具体步骤 #

1. 定义带约束的包装函数
   写一个完全匹配optim()参数签名的函数，固定L-BFGS-B方法并传入你的约束上下限：

   custom_optim <- function(par, fn, gr = NULL, ..., method = "L-BFGS-B", 
                            lower = c(0, 0, 0), upper = c(1, 20, 1),
                            control = list(), hessian = FALSE) {
     # 调用原生optim并传入约束参数
     optim(par = par, fn = fn, gr = gr, ..., method = method, 
           lower = lower, upper = upper, control = control, hessian = hessian)
   }
   

2. 在局部环境中运行lme
   创建局部环境，将自定义函数命名为optim，然后在该环境内执行模型拟合：

   # 创建局部环境，继承全局环境的其他对象
   local_env <- new.env(parent = .GlobalEnv)
   # 在局部环境中替换optim函数
   assign("optim", custom_optim, envir = local_env)
   
   # 在局部环境中拟合带约束的模型
   fitLME_constrained <- eval(quote(
     lme(sqrt(CD4) ~ obstime + obstime:drug, random = ~ obstime | patient, data = aids,
         control = lmeControl(opt = "optim", optimMethod = "L-BFGS-B"))
   ), envir = local_env)
   
   # 查看约束后的模型结果
   summary(fitLME_constrained)
   

关键细节说明 #

• 自定义函数必须严格匹配原生optim()的参数顺序和名称，否则lme()调用时会因参数不匹配报错。
• lmeControl()中的optimMethod参数在这里只是占位，实际优化方法由自定义函数中的method指定。
• 局部环境替换不会影响全局的optim()函数，避免污染全局环境。

注意事项 #

• 要确保lower和upper的长度与模型待估计的参数数量一致（固定效应参数+方差协方差参数），可以先拟合无约束模型，通过coef(fitLME)和VarCorr(fitLME)查看参数数量。
• 如果约束过严导致模型无法收敛，需要调整约束范围，或者给lme()指定更合理的初始参数（通过start参数）。

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