[R-sig-ME] nlme, nlmer and gradient
Cole, Tim
tim.cole at ucl.ac.uk
Fri Nov 18 12:09:47 CET 2016
Resending with correct r-sig-mixed-models address
Hi Ben,
My CRAN sitar library implements Mary Lindstrom's 1995 growth curve model using nlme, a natural B-spline mean curve and 3 subject random effects. I'm now trying to port it to lme4 by adding the gradient attribute to the function. I'm fairly sure the derivatives are correct, but it fails with error: step factor reduced below 0.001 without reducing pwrss.
I've seen your Algal (non)-linear mixed model example<http://www.rpubs.com/bbolker/3423> where nlmer fails with the same error, which you attribute to the fragility of nlmer. Is that likely to be the explanation here too?
The curious thing is that the same model fits fine in nlme, unless I include the gradient attribute in which case it fails with the equivalent error: step halving factor reduced below minimum in PNLS step. I'm wondering why the "fragility" of nlmer should also show itself in nlme when the gradient is provided, but not otherwise.
Some example code follows, using data from the heights dataframe in sitar and a spline curve with 6 df.
library(sitar)
library(lme4)
# nlme function without gradient
fitnlme <- function(x,s1,s2,s3,s4,s5,s6,a,b,c) {
c(
(cbind(s1,s2,s3,s4,s5,s6) * ns((x - b) * exp(c), knots=knots, B=bounds))
%*% matrix(rep(1,df), ncol=1)
) + a
}
# nlmer function with gradient
fitnlmer <- function(x,s1,s2,s3,s4,s5,s6,a,b,c) {
nsx <- function(xx, bb, cc) ns((xx - bb) * exp(cc), knots=knots, B=bounds)
ey <- function(mat1, mat2) c((mat1 * mat2) %*% matrix(rep(1, df), ncol=1))
css <- cbind(s1,s2,s3,s4,s5,s6)
cns <- nsx(x, b, c)
dimnames(cns)[[2]] <- dimnames(css)[[2]]
y <- ey(css, cns)
delta <- 0.1
dydb <- (ey(css, nsx(x, b + delta, c)) - y) / delta
dydc <- (ey(css, nsx(x, b, c + delta)) - y) / delta
attr(y, 'gradient') <- cbind(cns, a=1, b=dydb, c=dydc)
y + a
}
# set up model
x <- with(heights, age - mean(age))
bounds <- range(x) + 0.04 * c(-1,1) * diff(range(x))
df <- 6
knots <- quantile(x, 1:(df-1)/df)
# starting values
start <- setNames(coef(lm(height ~ ns(x, knots=knots, B=bounds), data=heights)), c('a', paste0('s', 1:df)))
start <- c(start[c(1:df+1, 1)], b=0, c=0)
# nlme fits using fitnlme without gradient
nlme1 <- nlme(height ~ fitnlme(x,s1,s2,s3,s4,s5,s6,a,b,c),
fixed = s1+s2+s3+s4+s5+s6+a+b+c ~ 1, random = a+b+c ~ 1 | id, data = heights, start = start)
# the equivalent nlme model in the sitar library also fits
sitar1 <- sitar(age, height, id, heights, 6)
# nlme fails using fitnlmer with gradient
nlme2 <- nlme(height ~ fitnlmer(x,s1,s2,s3,s4,s5,s6,a,b,c),
fixed = s1+s2+s3+s4+s5+s6+a+b+c ~ 1, random = a+b+c ~ 1 | id, data = heights, start = start)
# nlmer fails
nlmer1 <- nlmer(height ~ fitnlmer(x,s1,s2,s3,s4,s5,s6,a,b,c) ~ s1+s2+s3+s4+s5+s6+a+b+c + (a+b+c)|id, data=heights, start=start)
I'd appreciate your thoughts. And would adding starting values for theta likely make any difference?
Best wishes,
Tim
---
tim.cole at ucl.ac.uk<mailto: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
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list