[R-sig-ME] sitar and p=splines
Cole, Tim
t|m@co|e @end|ng |rom uc|@@c@uk
Fri Jul 19 00:35:42 CEST 2019
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]]
More information about the R-sig-mixed-models
mailing list