[R-sig-ME] nlme, nlmer and gradient

[Cole, Tim]

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]]