nlme vs nlme4

Mon, Jan 23, 2017 1:46 AM

Terry,

I've also recently reported that same error, at
https://stat.ethz.ch/pipermail/r-sig-mixed-models/2016q4/025213.html.

BTW where is your function dfun defined? - should it be dfun2?

You say you had to add z to the start vector init2, but you've defined it
as a fixed effect when it's random. Try defining init2 as a list with
components nlpars for the fixed effects p1, p2, p3 and p4, and theta for
var(z), and it might work better (note the error message refers to theta).

Best wishes,
Tim Cole

--- 
Tim.cole at ucl.ac.uk Phone +44(0)20 7905 2666 Fax +44(0)20 7905 2381
Population, Policy and Practice Programme
UCL Great Ormond Street Institute of Child Health, London WC1N 1EH, UK

Date: Thu, 19 Jan 2017 12:08:51 -0600
From: "Therneau, Terry M., Ph.D." <therneau at mayo.edu>
To: r-sig-mixed-models at r-project.org
Subject: [R-sig-ME] nlme vs nlme4

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"),
              function.arg=tfun2)
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,
              na.action=na.omit)
nfit1
  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