Create Spline from Coefficients and Knots in GAMLSS

Question

In the R package GAMLSS, it is possible to model a random variable $Y$ as a smoothed non-parametric function of some predictor $x$.

One option for such a function is the penalised spline using y~pb(x). This outputs a list of coefficients and knots which, combined with a set of basis splines, results in a smooth function of $x$.

How can one recreate the spline function, given the coefficients and the knots? (Preferably without having to write my own b-spline generating function).

For example:

library(gamlss)
set.seed(9876)
xstart <- 0
xend <- 100
datan <- 20000

seq <- seq(xstart, xend)
mean <- sapply(seq, function(x){0.5+0.2*sin(x/10)})
xs <- ceiling(runif(datan, xstart, xend))
ys <- sapply(xs, function(x){rnorm(1, mean = mean[x], sd = 0.1)})

m1 <- gamlss(ys~pb(xs))
plot(xs, ys)
lines(seq, mean, col="red")
lines(xs[order(xs)], fitted(m1)[order(xs)], col="green")

intercept <- m1$mu.coefficients[1] # 0.5495853
weight <- m1$mu.coefficients[2] # -0.0002851018
coefficients <- c(m1$mu.coefSmo[[1]]$coef) # c(-0.170704842, -0.066451626,  0.026591530,  0.119289203,  0.159657021,  0.149185418,  0.086505094,  0.003904402, -0.100156999, -0.188811997, -0.238717366, -0.237884900, -0.197802945, -0.090559794,  0.012576273,  0.101003289,  0.169210741,  0.181836117,  0.143883546,  0.061281663, -0.038608572, -0.136215586, -0.232871483)
knots <- m1$mu.coefSmo[[1]]$knots # c(-5.039,  0.010,  5.059, 10.108, 15.157, 20.206, 25.255, 30.304, 35.353, 40.402, 45.451, 50.500, 55.549, 60.598, 65.647, 70.696, 75.745, 80.794, 85.843, 90.892, 95.941 100.990 106.039)

How can I obtain the green function, knowing only the intercept, weight, coefficients and knots? I currently plot this function using fitted(m1). However, this is simply a list of $y$ values for the originally inputted list of $x$ values, it is not a function which gives $y$ for any new $x$.

Answer 1

the pb() function fits P-splines as described by Eilers and Marx (1996): B-splines on equally spaced knots and finite difference penalties. In the same paper there are some code chunks that show how to fit the model (in the appendix if I remember well). I do not know about the details of the fitted.gamlss method but the code below should be helpful (look at the blue line in the plot).

To compute the B-spline bases, I use the function splineDesign from splines package (the same function is used in the reference I menitoned above if I remember well).

To get $\hat{y}$ for a new value of $x$, you can just compute the corresponding value of the splineDesign function and use the coefficeints estimated before (see last line of the code and the green dot)

# B-splines stuffs - observed xs
ndx   = 20
deg   = 3
xr    = max(xs)
xl    = min(xs)
xmax  = xr + 0.01 * (xr - xl)
xmin  = xl - 0.01 * (xr - xl)
dt    = (xmax - xmin) / ndx
knots = seq(xmin - deg * dt, xmax + deg * dt, by = dt)
B     = splineDesign(knots = knots, x = xs, ord = deg + 1, derivs = 0,outer.ok = TRUE)

# Compute smooth
ff    = intercept + weight * xs + B %*% coefficients 

# New x-value
xn    = 35
Bn    =  splineDesign(knots = knots, x = xn, ord = deg + 1, derivs = 0,outer.ok = TRUE)
fn    = intercept + weight * xn + Bn %*% coefficients 

# Plot Results
lines(xs[order(xs)], ff[order(xs)], col = 'blue', lty = 2, lwd = 2)
points(xn, fn, col = 'green', pch = 16, cex = 1.5)

Answer 2

Here's a function call to generate new $y$ values for any $x$.

predict(m1, newdata = data.frame(xs = new_x), data = data.frame(xs=xs))

Note that my xs is from the xs in your example. In general, you'll need to modify this to provide the original data to the predict function, with the original names.