DC22
DC22

Reputation: 23

Nonlinear mixed models: what am I doing wrong?

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!

Upvotes: 2

Views: 519

Answers (1)

Ben Bolker
Ben Bolker

Reputation: 226731

Here's how far I've gotten:

  • the exponential rates are supposed to be specified as logs of the rates (to make sure that the rates themselves stay positive, i.e. that we have exponential decay curves rather than growth curves)
  • I simplified the model significantly, taking out the random effects in T1 and T2.
nlmer.ssbiexp<- nlmer(CD~SSbiexp(TIME, T1, R1, T2, R2)~
                        (R1|ID) + (R2|ID),
                      data=df,
                      start= c(T1= 4.98, R1=log(0.41), T2=5.25, R2=log(0.07)))

At this point I would try:

  • specifying T2 as 10-T1 to simplify the model (not sure if this will actually work)
  • seeing if you can build back to more complex versions of the model that actually work

Upvotes: 2

Related Questions