使用 lme4 中的 nlmer 时出现错误“Downdated VtV is not positive definite”

问题描述 投票:0回答:0

我正在尝试在 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)
}

模型数据如下所示: model data table 总共有 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 的反向求解中的奇异性”)。任何建议将不胜感激!

r lme4 mixed-models nlme
© www.soinside.com 2019 - 2024. All rights reserved.