[R-sig-ME] nlme vs nlme4

Therneau, Terry M., Ph.D. therneau at mayo.edu
Thu Jan 19 19:08:51 CET 2017

I have a fit that works in nlme and fails in nlme4; I'd like to convert to the latter 
since I now need to use case weights and varFixed() + nlme leads to a lot of "wrong 
length" warnings.  (Google search suggests that this has sometimes been asked about but I 
don't find any anwsers).

As always, my first hypothesis is that I have some silly setup mistake.

# nlme code
tfun2 <- function(age, p1, p2,p3,p4, z)
     p1 + p2*(age+z- 60) + .5*(p4-p2)*((age+z-p3) + sqrt((age+z-p3)^2 +20))

dfun2 <- deriv(body(tfun2), c("p1", "p2", "p3", "p4", "z"),
init <- c(p1=1.2, p2=.03, p3=1, p4=1)

nfit1 <- nlme(y ~ dfun(age, p1, p2, p3, p4, z),
               fixed= p1 + p2 + p3 + p4  ~ 1,
               random= z~ 1|clinic, data=mydata, start=init,
   Fixed: p1 + p2 + p3 + p4 ~ 1
          p1          p2          p3          p4
  1.21516330  0.00237005 73.58059382  0.08424922

Random effects:
  Formula: z ~ 1 | clinic
                z   Residual
StdDev: 9.418258 0.06865345


# nlmer code
init2 <- c(p1=1.2, p2=.03, p3=1, p4=1, z=0)
mfit1 <- nlmer(y ~ dfun(age, p1, p2, p3, p4, z) ~
               0+ p1 + p2 + p3 + p4 + (0+z | clinic),
                data=mydata, start=init2, verbose=TRUE)
Error in devfun(rho$pp$theta) : Downdated VtV is not positive definite


I had to add z=0 to the start vector or else it complains that z is not found, apparently 
the derivative component of dfun is not enough to inform.

For those who want to know more:  This is essentially a change point model fit with a 
hyperbola.  The biomarker in question rises slowly with age, but for a substantial subset 
there is a fairly abrupt upward change in the slope sometime between the ages of 60 and 
90.  The trend line is 1.2 + .0024*(age-60), the average age of turn is 73.5, and the new 
slope is .086, and the std of turning age is 9.4.  The constant "20" in the formula is an 
arbitrary value that controls how sharp a turn the hyperbola makes.

I can provide an anonymized data set to someone (if this isn't just a silly error).  There 
are about 3000 data points.

Terry Therneau

More information about the R-sig-mixed-models mailing list