How to find the derivative of bspline fit in R by using splines2

[CourseraXRay]

Greetings! I have a question concerns finding the derivative (or gradient) of bspline fit.

I am trying to estimate the 1st derivative of my bspline fitting. My code and results are included below.

require(splines2)

x = seq(0,10,by=0.1)
y = x^3+x^2+2*x+1
y_der = 3*(x^2)+2*x+2   # Manually calculated true derivative

mod <- lm(y ~ bSpline(x, knots = quantile(x, probs = seq(0, 1, 0.25)), degree = 3L, intercept = TRUE))
pdat <- data.frame(x = seq(min(x), max(x), length = 100))

yhat = predict(mod, newdata = pdat)

## Plot
plot(y ~ x)
lines(y_der ~ x, lwd = 2, col = "blue")
lines(yhat ~ x, data = pdat, lwd = 5, col = "green")

mod_der <- lm(y ~ bSpline(x, knots = quantile(x, probs = seq(0, 1, 0.25)), degree = 3L, deriv = 1L, intercept = TRUE))
yhat_der = predict(mod_der, newdata = pdat)
lines(yhat_der ~ x, data = pdat, lwd = 2, col = "red")

The green curve of bspline estimate almost perfectly fits the y=f(x) function. However the red curve of 1st derivative seems identical to the green curve, which is not what I expected. I thought it was somewhere near the blue parabola.

I understand that my R variable mod keeps the fitted bspline model. If I want to estimate the y value corresponds to a new x, I can just do predict(mod, data...).

I thought I just need to do the same thing if I want to estimate the derivative. I added deriv = 1L option in lm() to get the 1st derivatives of basis functions mod_der and then I called lm() again. I thought the result was the estimated curve of the gradient of bspline fit. But that doesn't work.

Is there an appropriate way to get the model of 1st (or even higher order) derivative of bspline fitting? Thanks in advance.

[wenjie2wang]

@CourseraXRay See the following code chunk for the revised example:

require(splines2)

x = seq(0,10,by=0.1)
y = x^3+x^2+2*x+1
y_der = 3*(x^2)+2*x+2   # Manually calculated true derivative

mat <- bSpline(x, df = 6, degree = 3L, intercept = TRUE)
mod <- lm(y ~ mat - 1)

## given new x
new_x <- seq(min(x), max(x), length.out = 1e3)
pmat <- predict(mat, new_x)
pmat_der <- deriv(pmat)
yhat <- c(pmat %*% coef(mod))
yhat_der <- c(pmat_der %*% coef(mod))

## Plot
plot(y ~ x)
points(y_der ~ x, lwd = 2, col = "blue")
lines(yhat ~ new_x, lwd = 5, col = "green")
lines(yhat_der ~ new_x, lwd = 2, col = "red")

There are several things that you need to be aware of.

• There is no derivs argument for bSpline() for splines2 v0.4.3. I will consider adding it for consistency with mSpline() in the future version.
• An intercept term is considered by default in lm(). Thus, you may either use y ~ mat - 1 or set intercept = FALSE when creating mat.
• The internal knots specified by the argument knots should be placed inside the boundary. Therefore, you may specify df instead to place internal knots based on quantiles of x.

[CourseraXRay]

Thank you so much Wenjie. That works like a charm! Keep up the good work! Raymond

…

On Mon, Aug 2, 2021 at 9:22 PM Wenjie Wang ***@***.***> wrote: @CourseraXRay <https://github.com/CourseraXRay> See the following code chunk for the revised example: require(splines2) x = seq(0,10,by=0.1)y = x^3+x^2+2*x+1y_der = 3*(x^2)+2*x+2 # Manually calculated true derivative mat <- bSpline(x, df = 6, degree = 3L, intercept = TRUE)mod <- lm(y ~ mat - 1) ## given new xnew_x <- seq(min(x), max(x), length.out = 1e3)pmat <- predict(mat, new_x)pmat_der <- deriv(pmat)yhat <- c(pmat %*% coef(mod))yhat_der <- c(pmat_der %*% coef(mod)) ## Plot plot(y ~ x) points(y_der ~ x, lwd = 2, col = "blue") lines(yhat ~ new_x, lwd = 5, col = "green") lines(yhat_der ~ new_x, lwd = 2, col = "red") There are several things that you need to be aware of. - There is no derivs argument for bSpline() for splines2 v0.4.3. I will consider adding it for consistency with mSpline() in the future version. - An intercept term is considered by default in lm(). Thus, you may either use y ~ mat - 1 or set intercept = FALSE when creating mat. - The internal knots specified by the argument knots should be placed inside the boundary. Therefore, you may specify df instead to place internal knots based on quantiles of x. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#10 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADA6GGX7D463YB6MGBAS473T25HH5ANCNFSM5BNNPYMA> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&utm_campaign=notification-email> .