1

I have a scatter plot that I've fitted with a natural cubic spline with 4 knots using ns(), pictured below:

enter image description here

I would like to find the value of x where the slope is zero at the abline pictured below: enter image description here

Is there a straightforward way to do this in RStudio? I'm having trouble with trying to figure out what the natural cubic spline formula would be and haven't had experience with derivatives since college. If it's helpful, here is the summary output:

enter image description here

enter image description here

Thanks very much!

lilleyb
  • 13
  • 3

2 Answers2

0

In principle, one could simply differentiate the spline and solve for where the derivative is 0 (even if that solution requires some sort of root finding procedure). Alternatively, you could just use a bisection method with predictions from the spline.

More practically, you can evaluate the spline over a fine grid to approximate the derivative using finite differences. Here is an example

library(splines)

set.seed(0)


Simulate data


n <- 1000
x <- runif(n, min=0, max = 1)


We already know the derivative is 0 at x=0.5


y <- sin(pi*x) + rnorm(n, 0, 0.3)


plot(x, y)


knts <- c(0.2, 0.4, 0.6, 0.8)
fit <- lm(y~ns(x, knots = knts, Boundary.knots = c(0.1, 0.9)))


fine grid over which to evaluate the spline


h <- 1e-4
xx <- seq(0, 1, h)
yy <- predict(fit, newdata = list(x=xx))


lines(xx, yy, col='red', lwd=3)



est_yprime <- diff(yy)/h

Find the x which produces the smallest absolute yprime


x_min <- xx[order(abs(est_yprime))[1]]


Close enough, and hovers around the right answer when you run it


With different seeds


x_min
#> [1] 0.4611

Created on 2024-03-01 with reprex v2.0.2

Demetri Pananos
  • 36,121
0

Just piggybacking off the previous answer, there actually is a package in R designed to yield derivatives of natural splines, and you can pretty quickly obtain the derivatives of the fitted spline and then use a root finding package.

require(splines2)
require(rootSolve)

set.seed(1234)


n = 1000
x = rnorm(n);
basis = model.matrix(~naturalSpline(x,df=6,
                                    Boundary.knots = c(-4,4),
                                    knots = quantile(x,probs = seq(0.1,0.9,length.out=6))))
beta = rnorm(8,sd=3)


y = basis%*%beta + rnorm(n)
fit <- lm(y~-1+basis)


xeval = seq(min(x),max(x),0.01)
derivBasis = model.matrix(~-1+naturalSpline(xeval,
                                            Boundary.knots = c(-4,4),
                                            knots = quantile(x,probs = seq(0.1,0.9,length.out=6)),
                                            derivs = 1))
derivs = derivBasis%*%coef(fit)[-1]


deriv_func <- function(xpoint){
  naturalSpline(xpoint,
                Boundary.knots = c(-4,4),
                knots = quantile(x,probs = seq(0.1,0.9,length.out=6)),
                derivs = 1)%*%coef(fit)[-1]
}


zeroes = uniroot.all(deriv_func,interval = range(x))


plot(x,y,xlab='X',ylab='Y');abline(v=zeroes,lty=2)
lines(sort(x),fit$fitted.values[order(x)],col='red',lwd=5)
plot(xeval,derivs,xlab='X',ylab="df(X)/dx");abline(h=0);abline(v=zeroes,lty=2)

ObsFit FittedDeriv

user1848065
  • 345