2017-06-27 31 views
0

的数据可以发现hereR:错误与内尔德 - 米德优化选择初始值

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调用会遇到上述问题,因为它对起始值非常敏感。结果,nmkb调用从不算什么。

有没有办法规避这种情况?

回答

0

模型日志似乎需要最大化,但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) 
} 

虽然还是有很多的警告,有我的机器上没有更多的错误,并且算法收敛成功。

$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