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!