library(nlme)
mydat <-
structure(list(class = c(91L, 91L, 91L, 91L, 91L, 91L, 92L, 92L,
92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L, 92L,
92L, 93L, 93L, 93L, 94L, 94L, 94L, 94L, 94L, 94L, 94L, 94L, 94L,
94L, 94L, 94L, 94L, 95L, 95L, 95L, 95L, 95L, 95L, 95L, 95L, 95L,
95L, 95L, 95L), days = c(5.7, 11.1, 13.9, 15.3, 18.3, 18.9, 1.9,
2.1, 2.9, 3.4, 4.4, 5, 6.9, 10.4, 11.6, 13, 13.4, 15.7, 15.9,
17.3, 17.7, 19.4, 2.3, 12.6, 15.4, 2, 2.4, 4.9, 5.4, 5.7, 7,
7.1, 7.7, 9.1, 9.6, 12.6, 16.6, 17.1, 1, 2, 4.4, 5.6, 5.7, 10.4,
12.1, 12.9, 13, 15.6, 16.1, 18.6), outcome = c(3.31586988521325,
3.26226473964254, 3.26098046236747, 3.26086828734987, 3.26081582971007,
3.2608136610639, 1.54175217273846, 1.74336277564818, 2.48010804039677,
2.73455940271066, 2.86602619542132, 2.85753365781511, 2.81667209739959,
2.80430238913247, 2.80395479988755, 2.80383291979961, 2.8038189189449,
2.80379174103878, 2.80379113213262, 2.80378896928776, 2.80378871890839,
2.80378827755366, 1.96537220180574, 3.00124636046136, 3.00096700482166,
2.05608815148142, 2.44248026102198, 4.03918455971327, 4.08570704450138,
4.09781416453829, 4.09869791687544, 4.09759058744364, 4.09084045815843,
4.07921200433542, 4.07668896637335, 4.07081047825795, 4.06993724272757,
4.06991925387706, 1.07225026715462, 3.72090308875724, 4.93988448353623,
4.7209971449984, 4.70681751285687, 4.49435510591282, 4.48824648431355,
4.4870870783591, 4.48698191392076, 4.48574152736339, 4.48566703246688,
4.48551396890595), s = c(0.0693, 0.0693, 0.0693, 0.0693, 0.0693,
0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693,
0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693,
0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693,
0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693,
0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693, 0.0693,
0.0693, 0.0693, 0.0693, 0.0693, 0.0693)), .Names = c("class",
"days", "outcome", "s"), row.names = 1001:1050, class = "data.frame")
library(nlme)
obj_NM <- function(arg){
model = nlme(outcome ~ exp(beta1) * s * days / (s * days + exp(1 - beta3 * s * days)) +
exp(beta2) * exp(1 - beta4 * s * days) / (s * days + exp(1 - beta4 * s * days)),
data = mydat,
fixed = list(beta1 ~ 1, beta2 ~ 1, beta3 ~ 1, beta4 ~ 1),
random = list(class = pdDiag(list(beta1 ~ 1, beta2 ~ 1, beta4 ~ 1))),
start = list(fixed = c(beta1 = arg[1], beta2 = arg[2], beta3 = arg[3], beta4 = arg[4])), verbose = FALSE)
return(model$logLik)
}
control = list(fnscale = -1)
optim(par = c(1.310239, -4.668217, 17.01345, 2.402943), fn = obj_NM, hessian = TRUE, control = control)
运行上面的代码会给我错误和警告：
Error in chol.default((value + t(value))/2) :
the leading minor of order 2 is not positive definite In addition: Warning messages:
1: In nlme.formula(outcome ~ exp(beta1) * s * days/(s * days + exp(1 - :
Singular precision matrix in level -1, block 1
2: In nlme.formula(outcome ~ exp(beta1) * s * days/(s * days + exp(1 - :
Singular precision matrix in level -1, block 1
我的目标是获取 hessian 矩阵并调查为什么我的nlme模型可能不合适。我试图最大化我的目标函数，因此我设置fnscale = -1（文档说它应该是负数以便optim执行最大化）。但是，我不确定如何处理错误消息。有没有办法optim输出粗麻布矩阵？似乎是一个错误nlme阻止了它这样做。