How can I use cubic splines for extrapolation?-CodePudding

I am looking to use natural cubic splines to interpolate between some data points using stats::splinefun(). The documentation states:

"These interpolation splines can also be used for extrapolation, that is prediction at points outside the range of ‘x’. Extrapolation makes little sense for ‘method = "fmm"’; for natural splines it is linear using the slope of the interpolating curve at the nearest data point."

I have attempted to replicate the spline function in Excel as a review, which is working fine except that I can't replicate the extrapolation approach. Example data and code below:

library(stats)

# Example data
x <- c(1,2,3,4,5,6,7,8,9,10,12,15,20,25,30,40,50)
y <- c(7.1119,5.862,5.4432,5.1458,4.97,4.8484,4.7726,4.6673,4.5477,4.437,4.3163,4.1755,4.0421,3.9031,3.808,3.6594,3.663)
df <- data.frame(x,y)

# Create spline functions
splinetest <- splinefun(x = df$x, y = df$y, method = "natural")

# Create dataframe of coefficients
splinetest_coef <- environment(splinetest)$z
splinetest_coefdf <- data.frame(i = 0:16, x = splinecoef_inf$x, a = splinecoef_inf$y, b = splinecoef_inf$b, c = splinecoef_inf$c, d = splinecoef_inf$d)

# Calculate extrapolated value at 51
splinetest(51)

# Result:
# [1] 3.667414

Question: How is this result calculated?

Expected result using linear extrapolation from x = 40 and x = 50 is 3.663 (51 - 50) x (3.663 - 3.6594) / (50 - 40) = 3.66336

The spline coefficients are as follows at i = 50: a = 3.663 and b = 0.00441355... Therefore splinetest(51) is calculated as 3.663 0.0441355 How is 0.0441355 calculated in this function?

CodePudding user response：

Linear extrapolation is not done by computing the slope between a particular pair of points, but by using the estimated derivatives at the boundary ("closest point" in R's documentation). The derivatives at any point can be calculated directly from the spline function, e.g. to calculate the estimated first derivative at the upper boundary:

splinetest(max(df$x), deriv = 1)
[1] 0.004413552

This agrees with your manual back-calculation of the slope used to do the extrapolation.

As pointed out in the comments, plotting the end of the curve/data set with curve(splinetest, from = 30, to = 60); points(x,y) illustrates clearly the difference between the derivative at the boundary (x=50) and the line based on the last two data points (i.e. (y(x=50) - y(x=40))/10)