6
$\begingroup$

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)

enter image description here

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$.

$\endgroup$
3
+100
$\begingroup$

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)

enter image description here

$\endgroup$
3
  • $\begingroup$ Thanks Gi_F. Why does $xmax = xr + 0.01 * (xr - xl)$ and not $xmax = max(xs)$? Where does the $0.01$ come from? $\endgroup$ – LBogaardt Jan 11 '20 at 19:19
  • $\begingroup$ The output of Gamlss gives 23 knots, just as there are 23 coefficients. However, the splineDesign functions requires 27 knots. For evenly spaces knots, these can easily be obtained from the original 23 using: $$c(head(knots, n=2) - 2 * knots[2] + 2 * knots[1], \\knots, \\tail(knots, n=2) + 2 * knots[2] - 2 * knots[1])$$ However, could you explain why splineDesign requires 27? $\endgroup$ – LBogaardt Jan 11 '20 at 19:33
  • 1
    $\begingroup$ 1commemt)The xmax/xmin depends on the way the function pb() Is written in the gamlss package (it is a choice of the package authors)..it is not so important after all. 2comment)The number of knots, depend on the properties of the B-spline functions. In general, the total number of knots needed is ndx + 2*deg+ 1 (see also'value' session of splineDesign() keeping in mind ord = deg + 1). B-splines are presented in a really accessible way in par 2 of Eiles and Marx 1996. Hope this clarifies a bit my answer. $\endgroup$ – Gi_F. Jan 12 '20 at 7:41
0
$\begingroup$

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.

$\endgroup$
7
  • $\begingroup$ Not exactly what I am looking for (or I am misunderstanding it). This recreates the list of $y$ values for the originally inputted list of $x$ values (the xs in the code). This is similar to the fitted(m1). It does not seem to be able to give a $y$ for any new $x$. $\endgroup$ – LBogaardt Dec 20 '19 at 8:49
  • $\begingroup$ Sorry! I'll clarify. $\endgroup$ – eric_kernfeld Dec 20 '19 at 15:20
  • $\begingroup$ Actually, my answer was straight-up wrong, so I completely changed it. $\endgroup$ – eric_kernfeld Dec 20 '19 at 15:30
  • $\begingroup$ Thanks for changing your answer. It's pretty close to what I expected. Perhaps this is simply the best option available in the GAMLSS package. I had prefered not to need the m1 model, but only use the coefficients and the knots. This should be possible in principle, but it might not be implemented anywhere. $\endgroup$ – LBogaardt Dec 20 '19 at 17:49
  • 1
    $\begingroup$ It might be worth looking for the source code for the predict method if it's open source. $\endgroup$ – eric_kernfeld Dec 20 '19 at 19:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.