如何在R中为SSbiexp/SSasympOFF自启动函数的NLS设置OFFSET=100？

问题描述 #

我在R里用自启动函数估计模型参数，但没法把x=0时的y值（偏移量）约束为100。用到的模型如下：

• 双组分指数模型：

model1=nls(y ~ SSbiexp(x, A1, lrc1, A2, lrc2) ,
    data = dat)

• 简单指数模型：

model2=nls(y ~ SSasympOff(x, Asym, lrc, c0), 
     data = dat)

希望约束偏移量为100，其他参数让自启动函数自动估算。尝试给SSbiexp加偏移量时，用Indometh数据集做了示例：

model_biexp=nls(conc ~ SSbiexp(time,A1, lrc1, A2, lrc2), data = Indometh)

# 尝试设置偏移量为10
model_biexpOFFSET=nls(conc ~ SSbiexp(time,A1, lrc1, A2, lrc2) + offset(10), data = Indometh)

运行后报错：

Error in nlsModel(formula, mf, start, wts) :
singular gradient matrix at initial parameter estimates
In addition: Warning message:
In nls(conc ~ SSbiexp(Indometh$time, A1, lrc1, A2, lrc2) + offset(10), :
No starting values specified for some parameters.
Initializing ‘A1’, ‘lrc1’, ‘A2’, ‘lrc2’ to '1.'.
Consider specifying 'start' or using a selfStart model

解决方法 #

直接用offset()会破坏自启动函数的参数初始化逻辑，导致nls无法生成合理初始值，进而报错。要约束x=0时的y值，得从模型数学形式入手，调整参数或构造新模型：

1. 简单指数模型（SSasympOff）的约束 #

SSasympOff的公式是：y = Asym + (c0 - Asym)*exp(-exp(lrc)*x)，x=0时y=c0，所以直接把c0固定为100即可：

# 直接固定c0为100，利用自启动函数估算其他参数
model2_constrained <- nls(y ~ SSasympOff(x, Asym, lrc, c0 = 100), 
                          data = dat)

如果遇到初始化报错，手动提取初始值再拟合：

# 获取无约束模型的初始参数
start_vals <- getInitial(y ~ SSasympOff(x, Asym, lrc, c0), data = dat)
# 固定c0为100
start_vals$c0 <- 100
# 用显式公式拟合约束模型
model2_constrained <- nls(y ~ Asym + (100 - Asym)*exp(-exp(lrc)*x), 
                          data = dat, start = start_vals)

2. 双组分指数模型（SSbiexp）的约束 #

SSbiexp的公式是：y = A1*exp(-exp(lrc1)*x) + A2*exp(-exp(lrc2)*x)，x=0时y=A1+A2。要让这个值等于100，只需令A2=100-A1，再调整拟合逻辑：

方法一：基于原自启动函数的初始值调整

# 获取无约束模型的初始参数
start_biexp <- getInitial(conc ~ SSbiexp(time, A1, lrc1, A2, lrc2), data = Indometh)
# 根据约束A1+A2=100，调整A2的初始值
start_biexp$A2 <- 100 - start_biexp$A1
# 拟合约束模型，公式中替换A2为100-A1
model_biexp_constrained <- nls(conc ~ A1*exp(-exp(lrc1)*time) + (100 - A1)*exp(-exp(lrc2)*time), 
                               data = Indometh, 
                               start = start_biexp[c("A1", "lrc1", "lrc2")]) # 只保留需要估算的参数

方法二：自定义带固定偏移量的自启动函数

如果需要重复使用，基于原SSbiexp构造新的自启动函数：

# 自定义带固定偏移量的双组分指数自启动函数
SSbiexpFixedOffset <- selfStart(
  function(x, A1, lrc1, lrc2, offset_val = 100) {
    A1*exp(-exp(lrc1)*x) + (offset_val - A1)*exp(-exp(lrc2)*x)
  },
  function(mCall, data, LHS) {
    # 调用原SSbiexp的初始值生成逻辑
    init <- getInitial(SSbiexp(x, A1, lrc1, A2, lrc2) ~ y, data = data)
    # 返回需要估算的参数初始值
    list(A1 = init$A1, lrc1 = init$lrc1, lrc2 = init$lrc2)
  },
  c("A1", "lrc1", "lrc2")
)

# 用自定义函数拟合，指定偏移量为100
model_biexp_constrained <- nls(conc ~ SSbiexpFixedOffset(time, A1, lrc1, lrc2, offset_val = 100), 
                               data = Indometh)

错误原因分析 #

你之前尝试的SSbiexp(...) + offset(10)相当于把模型改成y = A1*exp(...) + A2*exp(...) +10，这时候x=0时y=A1+A2+10，和你想要的“x=0时y=10”不符。更关键的是，这种写法会让nls无法识别这是自启动模型，只能给参数默认设为1，导致初始值不合理，进而出现奇异梯度矩阵错误。

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