[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"),
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
More information about the R-sig-mixed-models
mailing list