R:Nelder-Mead 优化到 select 起始值时出错
R: Error with Nelder-Mead optimization to select starting values
可以找到数据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(arg){
s.p1 <- arg[1]
s.b1 <- arg[2]
s.p2 <- arg[3]
s.b2 <- arg[4]
model = nlme(log10copy ~ exp(p1-b1*day) + exp(p2-b2*day),
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(10,0.5,6,0.005), lower = c(5, 0.1, 5, 0.001), upper = c(15, 1, 10, 0.1))
运行上面的代码,我运行进了几个错误:
Error in nlme.formula(log10copy ~ exp(p1 - b1 * day) + exp(p2 - b2 * day), :
step halving factor reduced below minimum in PNLS step
In addition: Warning message:
In nlme.formula(log10copy ~ exp(p1 - b1 * day) + exp(p2 - b2 * day), :
Singular precision matrix in level -1, block 1
该模型与 par = c(10,0.5,6,0.005)
中的起始值非常吻合。但是,我认为随着随机算法开始使用 lower = c(5, 0.1, 5, 0.001), upper = c(15, 1, 10, 0.1)
中的其他起始值,nlme
调用 运行s 解决上述问题,因为它对起始值非常敏感。因此,nmkb
调用从来没有任何意义。
有什么办法可以避免这种情况吗?
模型的log-liklihood需要最大化,但是R中的很多优化过程都给出了最小化结果。所以要优化的函数必须是负对数似然。所以它应该是这样的:
myfun2 <- function(arg){
s.p1 <- arg[1]
s.b1 <- arg[2]
s.p2 <- arg[3]
s.b2 <- arg[4]
model = nlme(log10copy ~ exp(p1-b1*day) + exp(p2-b2*day),
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)
}
而且虽然还是有很多warning,但是在我的机器上已经没有报错了,算法收敛成功。
$par
[1] 13.460199068 0.848526199 7.764024099 0.001513636
$value
[1] -358.6631
$feval
[1] 197
$restarts
[1] 0
$convergence
[1] 0
$message
[1] "Successful convergence"
Warning messages:
1: In nlme.formula(log10copy ~ exp(p1 - b1 * day) + exp(p2 - b2 * day), :
Singular precision matrix in level -1, block 1
可以找到数据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(arg){
s.p1 <- arg[1]
s.b1 <- arg[2]
s.p2 <- arg[3]
s.b2 <- arg[4]
model = nlme(log10copy ~ exp(p1-b1*day) + exp(p2-b2*day),
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(10,0.5,6,0.005), lower = c(5, 0.1, 5, 0.001), upper = c(15, 1, 10, 0.1))
运行上面的代码,我运行进了几个错误:
Error in nlme.formula(log10copy ~ exp(p1 - b1 * day) + exp(p2 - b2 * day), :
step halving factor reduced below minimum in PNLS step
In addition: Warning message:
In nlme.formula(log10copy ~ exp(p1 - b1 * day) + exp(p2 - b2 * day), :
Singular precision matrix in level -1, block 1
该模型与 par = c(10,0.5,6,0.005)
中的起始值非常吻合。但是,我认为随着随机算法开始使用 lower = c(5, 0.1, 5, 0.001), upper = c(15, 1, 10, 0.1)
中的其他起始值,nlme
调用 运行s 解决上述问题,因为它对起始值非常敏感。因此,nmkb
调用从来没有任何意义。
有什么办法可以避免这种情况吗?
模型的log-liklihood需要最大化,但是R中的很多优化过程都给出了最小化结果。所以要优化的函数必须是负对数似然。所以它应该是这样的:
myfun2 <- function(arg){
s.p1 <- arg[1]
s.b1 <- arg[2]
s.p2 <- arg[3]
s.b2 <- arg[4]
model = nlme(log10copy ~ exp(p1-b1*day) + exp(p2-b2*day),
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)
}
而且虽然还是有很多warning,但是在我的机器上已经没有报错了,算法收敛成功。
$par
[1] 13.460199068 0.848526199 7.764024099 0.001513636
$value
[1] -358.6631
$feval
[1] 197
$restarts
[1] 0
$convergence
[1] 0
$message
[1] "Successful convergence"
Warning messages:
1: In nlme.formula(log10copy ~ exp(p1 - b1 * day) + exp(p2 - b2 * day), :
Singular precision matrix in level -1, block 1