I am working with a data set that is comprised of three columns: patient ID (ID), TIME, and cervical dilation (CD). I apologize in advance for being unable to share my data, as it is confidential, but I have included a sample table below. Each patient CD was recorded in time as they progressed through labor. Time is measured in hours and CD can be 1-10cm. The number of time points/CD scores vary from patient to patient. In this model t is set in reverse, where 10 cm (fully dilated) is set as t=0 for all patients. This is done so that all patients can be aligned at time of full dilation. My dataset has no NA's and all patients have 2 or more time points.
| ID | TIME | CD |
|---|---|---|
| 1 | 0 | 10 |
| 1 | 3 | 8 |
| 1 | 6 | 5 |
| 2 | 0 | 10 |
| 2 | 1 | 9 |
| 2 | 4 | 7 |
| 2 | 9 | 4 |
I know for this problem I need to use nonlinear mixed effects model. I know from literature that the function that defines this biological process is modeled best as a biexponential function of the form CD= Cexp(-At)+(10-C)exp(-Lt), where A is the active labor rate [cm/hour], L is the latent labor rate [cm/hour], C is the diameter of the cervix [cm] at the point where the patient transitions from latent to active labor, and t is time in hours.
I have tried using both nlmer() and nlme() to fit this data, and I have used both the self-start biexponential function SSbiexp() as well as created my own function and its deriv(). Each parameter C, A, and L should have a random effect based on ID. Previous work has shown that C~4.98cm, A~0.41cm/hr, and L~0.07cm/hr. When using the SSbiexp(), there is a term for the second exponential component that is labeled here as C2, but should be the same as the (10-C) component of my self-made biexponential function.
When using nlme() with SSbiexp() I receive the error: Singularity in backsolve at level 0, block 1
nlme_ssbiexp<- nlme(CD~SSbiexp(TIME, C, A, C2, L),
data= df,
fixed = C+A+C2+L ~1,
random= C+A+C2+L ~1|ID,
start= c(C=4.98, A=0.41, C2= 5.02, L=0.07))
When using nlme() with my self-made biexp function I receive the error: Maximum number of iterations (maxIter = 50) reached without convergence
start<- c(C=4.98, A=0.41, L=0.07)
biexp<- ~ C*exp(-A*t) + (10-C)*exp(-L*t)
nfun <- deriv(biexp, namevec=c("C", "A", "L"),
function.arg=c("t","C", "A", "L"))
nlme_my_biexp<- nlme(CD~nfun(TIME, C, A, L),
data= df,
fixed = C+A+L ~1,
random= C+A+L ~1|ID,
start= c(C=4.98, A=0.41, L= 0.07))
When using nlmer() with SSbiexp() I receive an error: Error in devfun(rho$pp$theta) : Downdated VtV is not positive definite
nlmer_ssbiexp<-nlmer(CD~SSbiexp(TIME, C, A, C2, L)~(C|PTID)+(A|PTID)+(C2|PTID)+(L|PTID),
data=df,
start= c(T1= 4.98, R1=0.41, T2=5.02, R2=0.07))
When using nlmer() with my self-made biexp function I receive the error: Error in devfun(rho$pp$theta) : prss{Update} failed to converge in 'maxit' iterations
start<- c(C=4.98, A=0.41, L=0.07)
biexp<- ~ C*exp(-A*t) + (10-C)*exp(-L*t)
nfun <- deriv(biexp, namevec=c("C", "A", "L"),
function.arg=c("t","C", "A", "L"))
nlmer.fit.fun <- nlmer(CD ~ nfun(TIME, C, A, L) ~ (C|ID)+(A|ID)+(L|ID),
data = df,
start = start)
I had some success using the last combination of nlmer() and my self-made biexp fucntion, but only when I reduced my random effects to only including (C|ID).
start<- c(C=4.98, A=0.41, L=0.07)
biexp<- ~ C*exp(-A*t) + (10-C)*exp(-L*t)
nfun <- deriv(biexp, namevec=c("C", "A", "L"),
function.arg=c("t","C", "A", "L"))
nlmer.fit.fun <- nlmer(CD ~ nfun(TIME, C, A, L) ~ (C|ID),
data = df,
start = start)
I have tried increasing maxit and MaxIter, but both continued to fail to converge. I have been unable to find any solutions on Stack Overflow that help fixed these issues, though I have seen them discussed in several threads. I have also tried flipping the time scale so that t=0 is not always associated with CD=10, but it did not change my issues. I am new to R, so I am hoping someone that is an expert in nlmer() or nlme() might know fixes for these error messages.
Ben Bolker- if you are reading this- I would really love to talk with you!
