使用R的nls()拟合口服药代动力学函数遇奇异梯度错误

问题背景 #

要拟合的口服药代动力学模型公式为：
[ C_{PO}(t) = \frac{k_a \cdot D \cdot (e^{-k \cdot t} - e^{-k_a \cdot t})}{V \cdot (k_a - k)} ]
其中D=100，使用R的nls()函数拟合时触发“singular gradient”错误，代码及对应数据如下：

dat <- data.frame(
        Time = c(10, 15, 20, 30, 40, 60, 90, 120, 180, 210, 240, 300, 360),
        C_PO = c(0,0.28,0.55,1.2,2,1.95,1.85,1.6,
                 0.86,0.78,0.6,0.21,0.18))

plot(dat$C_PO ~ dat$Time, data = dat, log = "y")

fit <- nls(C_PO~ka*100* (exp(-k*Time) - exp(-ka*Time))/(v*(ka- k)),
    data = dat,
    start = list(ka = 0.09, k = 0.01, v = 50))
# 报错：singular gradient

解决方法 #

1. 用minpack.lm包的nlsLM()替代nls() #

nls()默认的高斯-牛顿算法对初始值要求极高，而nlsLM()采用列文伯格-马夸尔特算法，对初始值鲁棒性更强，能有效规避奇异梯度问题。

代码示例：

# 首次运行需安装包
install.packages("minpack.lm")
library(minpack.lm)

# 重新拟合模型
fit_lm <- nlsLM(C_PO~ka*100* (exp(-k*Time) - exp(-ka*Time))/(v*(ka- k)),
                data = dat,
                start = list(ka = 0.09, k = 0.01, v = 50))

# 查看拟合结果
summary(fit_lm)

2. 优化初始值或重参数化模型 #

奇异梯度多因参数高度相关或初始值偏差过大导致，可通过以下方式优化：

• 分步拟合初始值：先取消除相（t>60）的数据拟合单指数模型C_PO ~ A*exp(-k*Time)，得到更准确的k初始值后，再代入原模型拟合ka和v。
• 参数重参数化：改写模型降低参数共线性，比如用清除率CL = V*k替换V（即V=CL/k），模型变为：

fit_reparam <- nlsLM(C_PO~ ka*100*(exp(-k*Time)-exp(-ka*Time))/( (CL/k)*(ka -k) ),
                     data = dat,
                     start = list(ka=0.09, k=0.01, CL=0.5))  # CL初始值=原V*k=50*0.01=0.5
summary(fit_reparam)

3. 处理数据中的0值 #

首个时间点的C_PO=0可能干扰梯度计算，可替换为极小值或直接删除该行数据：

# 方案1：替换0为极小值
dat$C_PO[dat$C_PO == 0] <- 1e-6

# 方案2：删除该行数据
# dat <- dat[-1, ]

# 重新拟合
fit_clean <- nlsLM(C_PO~ka*100* (exp(-k*Time) - exp(-ka*Time))/(v*(ka- k)),
                   data = dat,
                   start = list(ka = 0.09, k = 0.01, v = 50))

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