R中的非线性混合效应回归

机器算法验证 r 非线性回归 混合模式 增长模式
2022-03-12 21:21:10

令人惊讶的是，我无法使用 Google 找到以下问题的答案：

我有一些来自几个人的生物学数据，这些数据及时显示了大致的 sigmoid 生长行为。因此，我希望使用标准的逻辑增长对其进行建模

P(t) = k*p0*exp(r*t) / (k+p0*(exp(r*t)-1))

其中 p0 是 t=0 时的起始值，k 是 t->infinity 时的渐近极限，r 是增长速度。据我所见，我可以使用 nls 轻松建模（我缺乏理解：为什么我不能通过缩放时间和数据使用标准 logit 回归来建模类似的东西？编辑：谢谢尼克​​，显然人们这样做是为了比例，但很少http://www.stata-journal.com/article.html?article=st0147。关于此切线的下一个问题是该模型是否可以处理异常值 >1）。

现在我希望允许对三个参数 k、p0 和 r 产生一些固定的（主要是分类的）和一些随机的（个人 ID 和可能还有研究 ID）影响。nlme 是最好的方法吗？SSlogis 模型对于我正在尝试做的事情似乎是明智的，对吗？以下任何一个都是一个明智的模型吗？我似乎无法正确获取起始值，并且 update() 似乎只适用于随机效果，而不是固定效果 - 有什么提示吗？

nlme(y ~ k*p0*exp(r*t) / (k+p0*(exp(r*t)-1)), ## not working at all (bad numerical properties?)
            data = data,
            fixed = k + p0 + r ~ var1 + var2,
            random = k + p0 + r ~ 1|UID,
            start = c(p0=1, k=100, r=1))

nlme(y ~ SSlogis(t, Asym, xmid, scal), ## not working, as start= is inappropriate
            data = data,
            fixed = Asym + xmid + scal ~ var1 + var2, ## works fine with ~ 1
            random = Asym + xmid + scal ~ 1|UID,
            start = getInitial(y ~ SSlogis(Dauer, Asym, xmid, scal), data = data))

由于我是特别是非线性混合模型和一般非线性模型的新手，因此我将不胜感激一些阅读建议或指向新手问题的教程/常见问题解答的链接。

1个回答

我想分享一些自从提出这个问题以来我学到的东西。nlme 似乎是在 R 中对非线性混合效果进行建模的合理方法。从一个简单的基础模型开始：

library(nlme)
data <- groupedData(y ~ t | UID, data=data) ## not strictly necessary
initVals <- getInitial(y ~ SSlogis(t, Asym, xmid, scal), data = data)
baseModel<- nlme(y ~ SSlogis(t, Asym, xmid, scal),
    data = data,
    fixed = list(Asym ~ 1, xmid ~ 1, scal ~ 1),
    random = Asym + xmid + scal ~ 1|UID,
    start = initVals
)

然后使用更新来增加模型复杂度。start 参数使用起来有点棘手，可能需要一些修补才能弄清楚顺序。请注意 var1 对 Asym 的新固定效果如何遵循 Asym 的常规固定效果。

 nestedModel <- update(baseModel, fixed=list(Asym ~ var1, xmid ~ 1, scal ~ 1), start = c(fixef(baseModel)[1], 0, fixef(baseModel)[2], fixef(baseModel)[3]))

lme4 似乎对我的数据集中的异常值更健壮，并且似乎为更复杂的模型提供了更可靠的收敛。但是，缺点似乎是需要手动指定相关的似然函数。以下是 var1（二进制）对 Asym 具有固定效应的逻辑增长模型。您可以以类似的方式在 xmid 和 scal 上添加固定效果。请注意使用双重公式指定模型的奇怪方式作为结果 ~ 固定效应 ~ 随机效应。

library(lme4) ## careful loading nlme and lme4 concurrently
customLogitModel <- function(t, Asym, AsymVar1, xmid, scal) {
    (Asym+AsymVar1*var1)/(1+exp((xmid-t)/scal))
}

customLogitModelGradient <- deriv(
    body(customLogitModel)[[2]], 
    namevec = c("Asym", "AsymVar1", "xmid", "scal"), 
    function.arg=customLogitModel
)

## find starting parameters
initVals <- getInitial(y ~ SSlogis(t, Asym, xmid, scal), data = data)

# Fit the model
model <- nlmer(
    y ~ customLogitModelGradient(t=t, Asym, AsymVar1, xmid, scal, var1=var) ~ 
    # Random effects with a second ~
    (Asym | UID) + (xmid | UID) + (scal | UID), 
    data = data, 
    start = c(Asym=initVals[1], AsymVar1=0, xmid=initVals[2], scal=initVals[3])
)
其它你可能感兴趣的问题
上一篇使用 Metropolis-Hastings 算法的 MCMC：选择提案 下一篇满足详细平衡的 MCMC 是否会产生平稳分布？