Determining derivatives from GAM smooth object

Question

I have a quite simple time series data set consisting of annual averages of a singe variable ("AVERAGE"). I wish to investigate the rate of change (1st derivative) and acceleration (2nd derivative) and associated standard errors of the "trend" component of the time series. I have obtained the "trend" using the GAM and PREDICT functions of MGCV simply as follows:

A <- gam(AVERAGE ~ s(YEAR), data=DF, na.action=na.omit)
B <- predict(A, type="response", se.fit=TRUE)

I have determined derivatives through 2 separate methods, applying a high DoF cubic smooth spline and via first and second differences (lightly smoothed) and bootstrapping to approximate errors with both producing comparable results.

I note that the "gam.fit3" function facilitates determining upto 2nd order derivatives but is not called directly. I also note that using "predict.gam" with the type "lpmatrix" facilitates derivatives of smooths. I would like to use the "GAM" function directly to calculate the 1st and 2nd derivatives but am not sufficiently skilled to calculate or extract these derivatives. I tried to reconfigure Wood's example at the end of the "Predict.gam" help page for one variable but with no success. Any help to get me headed in the right direction would be terrific. Thanks Phil.

Other option is to use the emtrends function of the emmeans package - that one will give you the first derivative plus confidence intervals (and using argument max.degree also 2nd derivatives I believe if need be). — Tom Wenseleers, Jul 09 '21 at 15:06

score 6 · Answer 1 · answered Jan 08 '13 at 03:03

The example from predict.gam uses finite differences to approximate the derivatives of the smoothed terms

Here is an example to do this for a single predictor model. This is more straightforward that the example from the help.

A <- gam(AVERAGE ~ s(YEAR), data=DF, na.action=na.omit)
# new data for prediction
newDF <- with(DF, data.frame(YEAR = unique(YEAR)))
# prediction of smoothed estimates at each unique year value
# with standard error    
B <- predict(A,  newDF, type="response", se.fit=TRUE)


# finite difference approach to derivatives following
# example from ?predict.gam

eps <- 1e-7
X0 <- predict(A, newDF, type = 'lpmatrix')


newDFeps_p <- newDF + eps

X1 <- predict(A, newDFeps_p, type = 'lpmatrix')

# finite difference approximation of first derivative
# the design matrix
Xp <- (X0 - X1) / eps

# first derivative
fd_d1 <- Xp %*% coef(A)

# second derivative
newDFeps_m <- newDF - eps

X_1 <- predict(A, newDFeps_m, type = 'lpmatrix')
# design matrix for second derivative
Xpp <- (X1 + X_1 - 2*X0)  / eps^2
# second derivative
fd_d2 <- Xpp %*% coef(A)

If you are using boot strapping to get the confidence intervals, you should be able to get confidence intervals on these approximations.

It's so exciting to hear one talking about first differences. So very Newtonian. — IRTFM, Jan 08 '13 at 04:55
Many thanks to @mnel for the time and effort taken to consider the problem and provide sample script. The code works perfectly however, with the data I am using annual averages with relatively small accelerations, the second differences become unstable with the use of such a small time slice (eps = 1e-7). I have undertaken some trial and error and the script works perfectly with eps in the range 0.1 to 0.01. Results compare perfectly to smooth spline derivatives. Many thanks Phil. — Phil W, Jan 10 '13 at 06:31
This is a bit late, but is the following correct? Xp <- (X0 - X1) / eps Should this not be (X1 - X0)? — David M, Jan 21 '20 at 20:43

Determining derivatives from GAM smooth object

1 Answers1

Linked