[R-sig-ME] sitar and p=splines

[Cole, Tim]

Hi,

My sitar mixed model package fits a natural spline mean curve to growth data, with three subject random effects reflecting random intercept (size), random time intercept (timing) and random time scaling (intensity).

I’m keen to modify the code to fit a penalised or p-spline rather than a b-spline, by including a smoother random effect, but I’m struggling to get the notation right.

Here is code to fit the sitar model using nlme with a b-spline (ns) smoother:

# first create data frame with centred x = age, y = height, id = id
 # install.packages('sitar')
 library(sitar)
  data <- setNames(heights[, 1:3], c('id', 'x', 'y'))
  data$x <- drop(scale(data$x, scale=FALSE))

# fit sitar model with 4 df
  { df <- 4
    knots <- with(data, quantile(x, (1:(df - 1)) / df))
    bounds <- with(data, range(x) + 0.04 * c(-1, 1) * diff(range(x)))
    mat <- matrix(rep(1, df), ncol = 1)
    start <- lm(y ~ splines::ns(x, knots = knots, Bound = bounds), data=data)
    start <- c(coef(start)[c(2:(df + 1), 1)], 0, 0)
    ey <- function(x,s1,s2,s3,s4,a,b,c) {
      x <- (x - b) * exp(c)     # subject-specific location and scale for x
      a + drop((cbind(s1,s2,s3,s4) * splines::ns(x,k=knots,B=bounds))%*%mat)
    }
    sitar1 <- nlme(y ~ ey(x,s1,s2,s3,s4,a,b,c),
         fixed = s1+s2+s3+s4+a+b+c ~ 1,
         random = a+b+c ~ 1 | id,                       ######## random term 1
         data = data, start = start)
  }

  # can fit the same sitar model using sitar package
  sitar2 <- sitar(x, y, id, data, 4)
  all.equal(sitar1, sitar2) # similar though not identical

  # now fit a single p-spline curve to the data using a random effect smoother
  # three p-spline functions from Paul Eilers

  tpower <- function(x, t, p)
    # Truncated p-th power function
    (x - t) ^ p * (x > t)

  bbase <- function(x, xl = min(x), xr = max(x), nseg = 10, deg = 3) {
    # Construct B-spline basis
    dx <- (xr - xl) / nseg
    knots <- seq(xl - deg * dx, xr + deg * dx, by = dx)
    P <- outer(x, knots, tpower, deg)
    n <- dim(P)[2]
    D <- diff(diag(n), diff = deg + 1) / (gamma(deg + 1) * dx ^ deg)
    B <- (-1) ^ (deg + 1) * P %*% t(D)
    B }

  makeZ <- function(x, xl = min(x), xr = max(x), nseg = 10) {
  # Construct smoothing random effect
    B <- bbase(x, xl, xr, nseg)
    D <- diff(diag(ncol(B)), diff = 2)
    Q <- solve(D %*% t(D), D)
    Z <- B %*% t(Q)
    Z
  }

  # creat a one group factor
  data$one <- 1

  # fit p-spline to group one with smoother random effect
  lm1 <- lme(y ~ x, random = list(one = pdIdent(~makeZ(x) - 1)), data = data)      ###### random term 2

  # plot the smooth mean curve for lm1 as points
  with(data, plot(x, fitted(lm1)))
  # and compare with sitar mean curve
  lines(sitar2, 'd', col=2)
  # rather different in shape as sitar adjusts for subject timing and intensity


The key question is how to combine the two approaches, i.e. to include the smooth random effect in sitar by combining these random terms:
  # random = list(one = pdIdent(~makeZ(x) - 1))
  # random = a+b+c ~ 1 | id
They obviously need to be a list, but the correct notation is not obvious.

A second key question is whether the random effect smoother makeZ(x) needs to be updated during the fit, as it depends on x, b and c which are being updated within ey. If it does I don’t think it can be done in nlme.

I’d value people’s thoughts.

Best wishes,
Tim
--
Population Policy and Practice
UCL Great Ormond Street Institute of Child Health,
30 Guilford Street, London WC1N 1EH, UK


	[[alternative HTML version deleted]]