[R-sig-ME] nlme vs nlme4

[Vito Michele Rosario Muggeo]

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