[R-sig-ME] nlme vs nlme4
Vito Michele Rosario Muggeo
vito.muggeo at unipa.it
Sat Jan 21 20:25:53 CET 2017
dear Terry,
I am sorry, but I don't have an answer to your question. However, as
it could be of interest for you, I would like to point out a paper of
mine dealing with piecewise (segmented) mixed models (or random
changepoint) in a likelihood-based framework. The method uses no
smooth approximation and it also allows random effects in the
changepoint parameter. I have written some functions to fit such
models in R. They should be included in the segmented package, but
currently are available on RG. Here the links
http://journals.sagepub.com/doi/abs/10.1177/1471082X13504721 #paper
https://www.researchgate.net/publication/263222159
https://www.researchgate.net/publication/292629179 #some notes about
the R functions
https://www.researchgate.net/publication/292986444 #R code with example
The code on RG works for basic examples, but I have an updated
version. Let me know if you are interested in.
I hope this helps you,
best,
vito
"Therneau, Terry M., Ph.D." <therneau at mayo.edu> ha scritto:
> 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
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
More information about the R-sig-mixed-models
mailing list