Are splines overfitting the data?

My problem: I recently met a statistician that informed me that splines are only useful for exploring data and are subjected to overfitting, thus not useful in prediction. He preferred exploring with simple polynomials ... As I’m a big fan of splines, and this goes against my intuition I’m interested in finding out how valid these arguments are, and if there is a large group of anti-spline-activists out there?

Background: I try to follow Frank Harrell, Regression Modelling Strategies (1), when I create my models. He argues that restricted cubic splines are a valid tool for exploring continuous variables. He also argues that the polynomials are poor at modelling certain relationships such as thresholds, logarithmic (2). For testing the linearity of the model he suggests an ANOVA test for the spline:

$$H_0: \beta_2 = \beta_3 = \dotsm = \beta_{k-1} = 0$$

I’ve googled for overfitting with splines but not found that much useful (apart from general warnings about not using too many knots). In this forum there seems to be a preference for spline modelling, Kolassa, Harrell, gung.

I found one blog post about polynomials, the devil of overfitting that talks about predicting polynomials. The post ends with these comments:

To some extent the examples presented here are cheating — polynomial regression is known to be highly non-robust. Much better in practice is to use splines rather than polynomials.

Now this prompted me to check how splines would perform in with the example:

library(rms)
p4 <- poly(1:100, degree=4)
true4 <- p4 %*% c(1,2,-6,9)
days <- 1:70

set.seed(7987)
noise4 <- true4 + rnorm(100, sd=.5)
reg.n4.4 <- lm(noise4[1:70] ~ poly(days, 4))
reg.n4.4ns <- lm(noise4[1:70] ~ ns(days,4))
reg.n4.4rcs_ols <- ols(noise4[1:70] ~ rcs(days,5))

plot(1:100, noise4)
nd <- data.frame(days=1:100)
lines(1:100, predict(reg.n4.4, newdata=nd), col="orange", lwd=3)
lines(1:100, predict(reg.n4.4ns, newdata=nd), col="red", lwd=3)
lines(1:100, predict(reg.n4.4rcs_ols, newdata=nd), col="darkblue", lwd=3)

legend("top", fill=c("orange", "red","darkblue"),
legend=c("Poly", "Natural splines", "RCS - ols"))

Gives the following image: In conclusion I have not found much that would convince me of reconsidering splines, what am I missing?

1. F. E. Harrell, Regression Modeling Strategies: With Applications to Linear Models, Logistic Regression, and Survival Analysis, Softcover reprint of hardcover 1st ed. 2001. Springer, 2010.
2. F. E. Harrell, K. L. Lee, and B. G. Pollock, “Regression Models in Clinical Studies: Determining Relationships Between Predictors and Response,” JNCI J Natl Cancer Inst, vol. 80, no. 15, pp. 1198–1202, Oct. 1988.

Update

The comments made me wonder what happens within the data span but with uncomfortable curves. In most of the situations I'm not going outside the data boundary, as the example above indicates. I'm not sure this qualifies as prediction ...

Anyway here's an example where I create a more complex line that cannot be translated into a polynomial. Since most observations are in the center of the data I tried to simulate that as well:

library(rms)
cmplx_line <-  1:200/10
cmplx_line <- cmplx_line + 0.05*(cmplx_line - quantile(cmplx_line, .7))^2
cmplx_line <- cmplx_line - 0.06*(cmplx_line - quantile(cmplx_line, .3))^2
center <- (length(cmplx_line)/4*2):(length(cmplx_line)/4*3)
cmplx_line[center] <- cmplx_line[center] +
dnorm(6*(1:length(center)-length(center)/2)/length(center))*10

ds <- data.frame(cmplx_line, x=1:200)

days <- 1:140/2

set.seed(1234)
sample <- round(rnorm(600, mean=100, 60))
sample <- sample[sample <= max(ds$$x) & sample >= min(ds$$x)]
sample_ds <- ds[sample, ]

sample_ds$$noise4 <- sample_ds$$cmplx_line + rnorm(nrow(sample_ds), sd=2)
reg.n4.4 <- lm(noise4 ~ poly(x, 6), data=sample_ds)
reg.n4.4rcs_ols <- ols(noise4 ~ rcs(x, 7), data=sample_ds)
AIC(reg.n4.4)

plot(sample_ds$$x, sample_ds$$noise4, col="#AAAAAA")
lines(x=ds$$x, y=ds$$cmplx_line, lwd=3, col="black", lty=4)

nd <- data.frame(x=ds$$x) lines(ds$$x, predict(reg.n4.4, newdata=ds), col="orange", lwd=3)
lines(ds\$x, predict(reg.n4.4rcs_ols, newdata=ds), col="lightblue", lwd=3)

legend("bottomright", fill=c("black", "orange","lightblue"),
legend=c("True line", "Poly", "RCS - ols"), inset=.05)

This gives the following plot: Update 2

Since this post I've published an article that looks into non-linearity for age on a large dataset. The supplement compares different methods and I've written a blog post about it.

• I don't really see where your statistician friend is coming from, to be honest. You can overfit with polynomials and splines just the same. Overfitting comes from your class of models having excessive capacity; what distinguishes the performance of various models is how they restrict their capacity. For (natural) splines, it is knot placement and number, for polynomials it is the degree.
– guy
Feb 1 '13 at 15:22
• @guy: That is also my belief, you can always overfit your data no matter what method you use. During my regression class my professor told me that polynomials bend where the majority of the data occur, thus making the extremes more unreliable. Though I haven't found any article supporting this claim. Feb 1 '13 at 16:07
• All of the curves in your first graph fail to fit the data at the far right hand side. Feb 5 '13 at 21:22
• Is the 'x' dimension in the above graphs time related? If it is, none of these methods are appropriate, because both are 'forward-looking' in the sense that they use neighboring points (on both sides) to model. Feb 7 '13 at 21:56
• @arielf: No the x is not intended as a time variable. I thought of it as some variable where we sample the maximum number of observations at the center. In my research we don't look into the future that much, I guess it's more in the field of inference than prediction. The variable is intended to be a cholesterol, blood pressure, BMI, or some other common continuous variable. Feb 8 '13 at 9:50

Overfitting comes from allowing too large a class of models. This gets a bit tricky with models with continuous parameters (like splines and polynomials), but if you discretize the parameters into some number of distinct values, you'll see that increasing the number of knots/coefficients will increase the number of available models exponentially. For every dataset there is a spline and a polynomial that fits precisely, so long as you allow enough coefficients/knots. It may be that a spline with three knots overfits more than a polynomial with three coefficients, but that's hardly a fair comparison.

If you have a low number of parameters, and a large dataset, you can be reasonably sure you're not overfitting. If you want to try higher numbers of parameters you can try cross validating within your test set to find the best number, or you can use a criterion like Minimum Description Length.

EDIT: As requested in the comments, an example of how one would apply MDL. First you have to deal with the fact that your data is continuous, so it can't be represented in a finite code. For the sake of simplicity we'll segment the data space into boxes of side $$\epsilon$$ and instead of describing the data points, we'll describe the boxes that the data falls into. This means we lose some accuracy, but we can make $$\epsilon$$ arbitrarily small, so it doesn't matter much.

Now, the task is to describe the dataset as sucinctly as possible with the help of some polynomial. First we describe the polynomial. If it's an n-th order polynomial, we just need to store (n+1) coefficients. Again, we need to discretize these values. After that we need to store first the value $$n$$ in prefix-free coding (so we know when to stop reading) and then the $$n+1$$ parameter values. With this information a receiver of our code could restore the polynomial. Then we add the rest of the information required to store the dataset. For each datapoint we give the x-value, and then how many boxes up or down the data point lies off the polynomial. Both values we store in prefix-free coding so that short values require few bits, and we won't need delimiters between points. (You can shorten the code for the x-values by only storing the increments between values)

The fundamental point here is the tradeoff. If I choose a 0-order polynomial (like f(x) = 3.4), then the model is very simple to store, but for the y-values, I'm essentially storing the distance to the mean. More coefficients give me a better fitting polynomial (and thus shorter codes for the y values), but I have to spend more bits describing the model. The model that gives you the shortest code for your data is the best fit by the MDL criterion.

(Note that this is known as 'crude MDL', and there are some refinements you can make to solve various technical issues).

• Thank you Peter for your answer. I've tried to wrap my head around MDL, especially how to apply it. It would be nice to have it explained based on one of my examples. As a non-statistician I like having things exemplified before I can understand the underlying logistics. The coin example in the Wiki-article didn't reach me... Feb 7 '13 at 19:39
• I've added an example. Feb 8 '13 at 10:16
• Thank you Peter for the example, it is much clearer to me now. Feb 8 '13 at 21:16

Statisticians have been arguing about polynomial fitting for ages, and in my experience, it comes down to this:

Splines are basically a series of different equations pieced together, which tends to increase the accuracy of interpolated values at the cost of the ability to project outside the data range. This is fine if you know your data is pure and coming from a consistent source and if you are trying to describe the likelihood of different values' presence within your range of values. However, we usually don't learn as much about the theoretical underpinnings driving our data, since a new spline starts when the old spline stops accurately describing the data. This makes prediction of values outside our data almost worthless.

Now, splines are not unique in this respect. Polynomial functions actually suffer from the same problem if we are just fitting the data and not using a theoretical framework for choosing the variables. Those who have a well-formed theory driving which variables to allow to vary and by how much will be more trusting of a complex polynomial function's ability to extrapolate predictions outside the data.

Many statisticians, though, are working with data without help from a pre-established theoretical framework, and this pushes some people towards simple polynomials. They reason that a less flexible function that fits the data is more likely to accurately predict values outside the data, because the function is less likely to be swayed by anomalies within the data. While I've had conversations about this with people who prefer simple polynomials, I've never gotten the feeling of an anti-spline group. It feels more like simple polynomials make some statisticians feel more comfortable about avoiding overfitting.

Disclaimer

Personally, I don't tend to use splines or simple polynomials with most of my data, because I work in a field with many pre-established theoretical frameworks. Also, I have usually observed the collection of the data and can get a decent grasp on what was driving the outcomes. In that case, I'm building more of a logical algorithm and testing the algorithm's fitness, rather than testing the fitness of a polynomial function. You can add this grain of salt to my answer.

• Polynomials are far more sensitive to anomalies within the data than splines are. An outlier anywhere in the data set has a massive global effect, while in splines the effect is local.
– guy
Feb 1 '13 at 15:18
• I see your point, and that is true if you are using a perfect information approach or do not have enough information about the nature of the data. Many statisticians (myself included) assume imperfect information and make an attempt to apply exclusion criteria based on known information before trying to fit the data. Dangerous outliers should then theoretically be excluded from the fitting attempt. If you don't have the known information about the nature of the data (and this is quite common), then you're stuck trying to work around the outliers. Feb 1 '13 at 15:32
• I would have to be better convinced that regression splines extrapolate more dangerously than polynomials. Feb 8 '13 at 13:59
• This isn't anything new. Rather, it's the difference seen between statistics done in early stages of understanding versus later stages of understanding. The more you understand a system, the less you rely on fitted functions and the more you rely on theoretical models. Feb 8 '13 at 14:40
• How about using restricted cubic splines, which constrain the functions to be linear in the outside of data points(I am reading Harrell's book). Anyway extrapolation is always suspicious. Think of an experiment which discovered the superconductivity or plasma. Theory should be proved via experiments! I think what functions to fit is more relevant to interpolation problem. Without theory, I guess you would not be able to pick just one model with predictors in error(also unknown distribution) and unknown distribution of y|x, even when you are give enough data. Nov 13 '18 at 12:58