我正在尝试在 R 中使用 lme4 (lme4_1.1-31) 实现 NLME 模型
vW_start = 0.045/0.042
vK_start = exp(0.9)*vW_start
vW_start = log(vW_start-1)
vK_start = log(vK_start-1)
b2_start = 1.718
b3_start = 0.968
b4_start = -0.131
m_start = (b4_start + b3_start)/(-sign(b2_start))
nlin.mdl.results<-nlmer(log_ratio ~ mdl(L, t, vK, vW, b2, b3, m) ~ (vK | id) + (vW | id) + (b2 | id) + (b3 | id) + (m | id), mdl_data, start = c(vK = vK_start, vW = vW_start, b2 = b2_start, b3 = b3_start, m = m_start))
我收到以下错误消息,“devfun(rho$pp$theta) 中的错误:已更新的 VtV 不是正定的”。
非线性模型函数定义为:
mdl<-function(L, cell_type, vK, vW, b2, b3, m){
dvK = exp(vK)
vK = exp(vK)+1
dvW = exp(vW)
vW = exp(vW)+1
S = log(vW-1)-log(vK-1)+L
dS_vK = -1/(vK-1)
dS_vW = 1/(vW-1)
M2 = b2+log(vK)-log(vW)+S
dM2_b2 = 1
dM2_vK = 1/vK
dM2_vW = -1/vW
dM2_S = 1
M3 = b3+log(vK)-log(vW)+S
dM3_b3 = 1
dM3_vK = 1/vK
dM3_vW = -1/vW
dM3_S = 1
sgn_b2 = b2/sqrt(b2^2)
b4 = (-sgn_b2*exp(m))-b3
db4_b3 = -1
db4_m = -sgn_b2*exp(m)
M4 = b4+log(vK)-log(vW)+S
dM4_b4 = 1
dM4_vK = 1/vK
dM4_vW = -1/vW
dM4_S = 1
S = (-1/6)*(t-2)*(t-3)*(t-4)*S
M2 = (1/2)*(t-1)*(t-3)*(t-4)*M2
M3 = (-1/2)*(t-1)*(t-2)*(t-4)*M3
M4 = (1/6)*(t-1)*(t-2)*(t-3)*M4
f = S + M2 + M3 + M4
#calculate derivatives
df_vK = ((-1/6)*(t-2)*(t-3)*(t-4)*(dS_vK*dvK))+
((1/2)*(t-1)*(t-3)*(t-4)*((dM2_vK*dvK)+(dM2_S*dS_vK*dvK)))+
((-1/2)*(t-1)*(t-2)*(t-4)*((dM3_vK*dvK)+(dM3_S*dS_vK*dvK)))+
((1/6)*(t-1)*(t-2)*(t-3)*((dM4_vK*dvK)+(dM4_S*dS_vK*dvK)))
df_vW = ((-1/6)*(t-2)*(t-3)*(t-4)*(dS_vW*dvW))+
((1/2)*(t-1)*(t-3)*(t-4)*((dM2_vW*dvW)+(dM2_S*dS_vW*dvW)))+
((-1/2)*(t-1)*(t-2)*(t-4)*((dM3_vW*dvW)+(dM3_S*dS_vW*dvW)))+
((1/6)*(t-1)*(t-2)*(t-3)*((dM4_vW*dvW)+(dM4_S*dS_vW*dvW)))
df_b2 = (1/2)*(t-1)*(t-3)*(t-4)*(dM2_b2)
df_b3 = ((-1/2)*(t-1)*(t-2)*(t-4)*(dM3_b3))+
((1/6)*(t-1)*(t-2)*(t-3)*(dM4_b4*db4_b3))
df_m = (1/6)*(t-1)*(t-2)*(t-3)*(dM4_b4*db4_m)
gradient = cbind(df_vK, df_vW, df_b2, df_b3, df_m)
colnames(gradient)<-c("vK", "vW", "b2", "b3", "m")
attr(f, "gradient")<-gradient
return(f)
}
模型数据如下所示: 总共有 480 个观察值(120 个 id,每个在 t = 1、2、3 和 4 时有一个观察值)
在阅读了有关此错误消息的一些信息后,我仍然无法完全理解它是源于我的模型规格还是 nlmer 中的其他选项。我还尝试使用 nlme 包实现模型:
gData = groupedData(log_ratio ~ 1|id, mdl_data)
nlin.mdl.results<-nlme(log_ratio ~ mdl(L, t, vK, vW, b2, b3, m), data = gData, fixed = list(vK ~ 1, vW ~ 1, b2 ~ 1, b3 ~ 1, m ~ 1), random = vK + b2 + b3 + m ~ 1|id, start = c(vK = vK_start, vW = vW_start, b2 = b2_start, b3 = b3_start, m = m_start))
这在返回以下错误之前运行了一些迭代:“nlme.formula 中的错误(log_ratio ~ mdl(L,vK,vW,:0 级,块 1 的反向求解中的奇异性”)。任何建议将不胜感激!