R:非线性混合效应模型的 Nelder-Mead 优化

R: Nelder-Mead optimization for nonlinear mixed effects models

数据找到here

library(nlme)
library(dfoptim)
dat0 <- read.table("aids.dat2",head=T)
dat1 <- dat0[dat0$day<=90, ]   # use only first 90-day data
dat2 <- dat1[!apply(is.na(dat1),1,any),]  # remove missing data 
aids.dat <- groupedData(lgcopy ~ day | patid, data=dat2)
aids.dat$log10copy = log10(aids.dat$lgcopy)

myfun2 <- function(s.p1, s.b1, s.p2, s.b2){
  model = nlme(log10copy ~ exp(p1-b1*day) + exp(p2-b2*day + 1),
               fixed = list(p1 ~ 1, b1 ~ 1, p2 ~ 1, b2 ~ 1),
               random = list(patid = pdDiag(list(p1 ~ 1, b1 ~ 1, p2 ~ 1, b2 ~ 1))),
               start = list(fixed = c(p1 = s.p1, b1 = s.b1, p2 = s.p2, b2 = s.b2)),
               data =aids.dat) 
  return(model$logLik)
}

nmkb(fn = myfun2, par = c(1.13, 0.25, 13.1, 1.3), lower = c(0.8, -0.3, 5, -2), upper = c(1.6, 0.6, 20, 13))

我正在尝试 运行 Nelder-Mead 优化算法(在 dfoptim 包中)为我的非线性混合效应模型找到好的起始值。我的 objective 函数是我模型中的 logLik。然后我调用 nmkb,其中 fn 是我的函数,par 是我的 4 个参数的值及其 lowerupper 边界。但是,我 运行 正在进入错误

Error in nlme.formula(log10copy ~ exp(p1 - b1 * day) + exp(p2 - b2 * day +  : 
  argument "s.b1" is missing, with no default  

似乎是语法问题,但这里到底出了什么问题?

myfun2 <- function(p1b1p2b2){
    s.p1 <- p1b1p2b2[1]
    s.b1 <- p1b1p2b2[2]
    s.p2 <- p1b1p2b2[3]
    s.b2 <- p1b1p2b2[4]
    model = nlme(log10copy ~ exp(p1-b1*day) + exp(p2-b2*day + 1),
                 fixed = list(p1 ~ 1, b1 ~ 1, p2 ~ 1, b2 ~ 1),
                 random = list(patid = pdDiag(list(p1 ~ 1, b1 ~ 1, p2 ~ 1, b2 ~ 1))),
                 start = list(fixed = c(p1 = s.p1, b1 = s.b1, p2 = s.p2, b2 = s.b2)),
                 data =aids.dat) 
    return(model$logLik)
}