6
$\begingroup$

I am trying to fit some time series data to a smoothing spline in R. However, it seems like the spline is fitting the data too perfectly, meaning overfitting. I was trying to figure out what settings to change to try and adjust the level of smoothing. I don't want to manually set the $\lambda$ parameter, because it seems like that should automatically be set according to some metric. I believe the default is generalized cross-validation, so that should work okay.

Here is some data and code. Can anyone tell me the correct way to apply the splines to the dataset.

I can manually change the number of Knots, but that seems a bit manual. I have a number of datasets to fit, so I don't want to manually fit the number of knots each time. Is there a better way to determine the penalty. I suppose having a twice differentiable function is optimal, so no sharp edges.

library(npreg)

y <- c(23.0, 27.0, 25.0, 25.0, 25.0, 22.0, 22.0, 21.0, 18.0, 16.0, 17.0, 17.0, 19.0, 19.0, 19.0, 20.0, 19.0, 18.0, 20.0, 19.0, 17.0, 21.0, 20.0, 16.0, 15.0, 16.0, 14.0, 14.0, 12.0, 14.0)

x <- c(1979, 1980, 1981, 1982, 1983, 1984, 1985, 1986, 1987, 1988, 1989, 1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008)

mod <- ss(x, y)
pred_y <-predict(mod, x)$y
plot(x,y)
lines(x, pred_y, lwd=2)

The corresponding picture is.

enter image description here

Improve this question
$\endgroup$
16
  • $\begingroup$ What is the default argument for the number of knots? $\endgroup$
    – Galen
    Sep 6, 2022 at 4:09
  • $\begingroup$ The number of knots is a hyperparameter. You can tune it. $\endgroup$
    – Galen
    Sep 6, 2022 at 4:10
  • 1
    $\begingroup$ Yeah, so I want to do a smoothed interpolation--so not fitting all of the points exactly. Then I will compute the derivative at each data point, and use those derivatives to fit a differential equation. So I need the derivatives at each point. $\endgroup$
    – krishnab
    Sep 6, 2022 at 4:20
  • 1
    $\begingroup$ You would probably want to specify one of df, spar or lambda in order to make the function smoother. $\endgroup$
    – Glen_b
    Sep 6, 2022 at 4:43
  • 1
    $\begingroup$ You may need to start with one or more plausible descriptions of the process you'd like to approximate and look at performance on those. Presumably you're interested in optimizing typical performance rather than trying to optimize based on the noise in one sample. $\endgroup$
    – Glen_b
    Sep 6, 2022 at 23:44

2 Answers 2

Reset to default
4
$\begingroup$

More specifically, you are right that something is being done to select the smoothing parameters for the spline and that by default this is GCV. It is known (from the spline / GAM literature) that GCV can undersmooth and I believe this is what you are seeing here.

Choosing another method, such as REML smoothness selection leads to a more reasonable fit:

# reusing objects from your post
m_reml <- ss(x, y, method = "REML")
p_reml <- predict(m_reml, x)$y
plot(x, y)
lines(x, pred_y, col = "red", lwd = 2)
lines(x, p_reml, col = "blue", lwd = 2)

produces

enter image description here

As you have a variable that is ordered in time, there is the additional complication that if there is some signal that is contaminated by autoregressive noise, this violates the assumptions used to select the smoothing parameters and can lead you to over fit. One option in that case would be to model the autocorrelated noise:

library("mgcv")
df <- data.frame(y = y, x = x)
m_gam <- gamm(y ~ s(x, k = 6), data = df, method = "REML",
              correlation = corAR1(form = ~ x))

p_gam <- predict(m_gam$gam, newdata = df)

plot(x, y)
lines(x, pred_y, col = "red", lwd = 2)
lines(x, p_reml, col = "blue", lwd = 2)
lines(x, p_gam, col = "green", lwd = 2)

enter image description here

In this instance it doesn't help (for some definition of "help") as the model has assigned all the variation to the autocorrelation process (the AR(1)) and the resulting trend is linear with some aspect of the model fit becoming no-positive definite - a sure sign that the model is over-fit or too complex. This often happens because a trend and an autocorrelation process like the AR(p) are not always identifiable from data

Improve this answer
$\endgroup$
4
  • $\begingroup$ thanks so much for this detailed response. This is very helpful. I will keep in mind what you said about REML. The point about autocorrelation is also important. I am not as familiar with the mgcv package, but does the gamm command first remove the linear trend and then try to fit an AR(1) process? I imagine this linear trend is estimated and then the series is detrended before fitting the AR(1), but just wanted to confirm that. $\endgroup$
    – krishnab
    Sep 17, 2022 at 3:36
  • $\begingroup$ It is interesting that the process fits the AR(1) process so well, but I get what you mean. So from the code it seems like it tried to fit up to AR(6), as per the $k$ parameter, but then columns 2-6 are all linearly dependent and hence the matrix is singular? Is that why the corresponding line is purely linear trend, since there is no trendiness left in the residuals to assign anywhere else. $\endgroup$
    – krishnab
    Sep 17, 2022 at 3:53
  • $\begingroup$ Only an AR(1) was fitted in the correlation matrix of residuals; the k relates to the size of the basis expansion of x, which resulted in 5 columns (basis functions) being added to the model matrix plus a constant term for the intercept. This is basically a Generalised Least Squares-like model, except we're treating the smooth random effects, so the linear part of the basis and the intercept are in the model matrix while the other 4 basis functions are actually in the model matrix for the random effects. All parameters are being estimated simultaneously, so the AR(1) isn't estimated later $\endgroup$
    – Gavin Simpson
    Sep 17, 2022 at 12:19
  • 1
    $\begingroup$ Ahh okay. I understand what you mean now. Thanks so much for your help on this question. I really appreciate it. $\endgroup$
    – krishnab
    Sep 17, 2022 at 23:04
8
$\begingroup$

The number of knots is a hyperparameter. You can tune it.

There are many approaches although for time series you will need to be cautious about how you partition the data. See this online textbook passage for guidance on splitting. Essentially, you'll have to sequentially break the data up into blocks.

Improve this answer
$\endgroup$
5
  • $\begingroup$ Oh yes, I don't need to worry about data leakage, as in the Hyndman book. I don't have a test set that I am trying to evaluate on, in this case. But definitely something to keep in mind when I get to that part. Thanks. $\endgroup$
    – krishnab
    Sep 6, 2022 at 4:19
  • $\begingroup$ @krishnab Okay. I am speculating that grid search with leave-one-out CV might be suitable since you don't appear to have lot of data. It might be worth looking up differentiation-specific smoothing methods (if you have not already). $\endgroup$
    – Galen
    Sep 6, 2022 at 4:28
  • $\begingroup$ yeah that makes sense. I can do that. In terms of penalties on the derivative, yeah I can include that. Usually I can penalize the difference between the values at $x_i$ and $x_{i+1}$, with a tuning parameter. This is similar to total variation denoising, if you are familiar with signal processing. I get what you mean though, I think I have enough to go on now. $\endgroup$
    – krishnab
    Sep 6, 2022 at 4:31
  • $\begingroup$ @kirshnab you mention that you have multiple datasets to fit. Are you expecting them to all have the same "number of knots"? If yes, consider holding out some of these datasets as test sets and using the rest as train sets to figure out a good parameter $\endgroup$
    – Shadi
    Sep 6, 2022 at 19:29
  • $\begingroup$ @Shadi yes, that makes sense. Yeah, I am just working on preprocessing the data at this point, but you are correct as well as Galen. I will have to think of a way to hold out some datasets or hold out some observations, etc. Something like that makes sense. $\endgroup$
    – krishnab
    Sep 6, 2022 at 20:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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