0

I am attempting to implement a NLME model using lme4 (lme4_1.1-31) in R.

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))

I receive the following error, "Error in devfun(rho$pp$theta) : Downdated VtV is not positive definite".

The nonlinear model function is defined as:

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)
}

The model data looks like the following: model data table There are 480 observations total (120 ids, each with an observation at t = 1, 2, 3, & 4)

After reading a bit about this error message, I still can't quite decipher if it's originating from the specifications of my model or other options within nlmer. I also tried implementing the model using the nlme package:

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))

This ran for some iterations before returning the following error: "Error in nlme.formula(log_ratio ~ mdl(L, vK, vW, : Singularity in backsolve at level 0, block 1". Any suggestions would be appreciated!

asingh
  • 1
  • 1
  • That is a lot of parameters. I would use `nlme` for non-linear mixed-effects models. You should start simple, e.g., by fitting a model with `nlsList` and inspecting the fit, then passing that model to `nlme` and adding the random effects. – Roland May 02 '23 at 05:49

0 Answers0