# Use of splines in parameter estimation

The question here is about a passage in the very frequently (10,000+) cited paper by Storey and Tibshirani "Statistical significance for genomewide studies" (https://doi.org/10.1073/pnas.1530509100).

For their method they need to estimate a parameter named $$\pi_0$$, which is more or less the overall false discovery rate, for which they compute some $$\hat \pi_0(\lambda)$$ ($$\lambda$$ being some bias-variance tradeoff parameter) for various values of $$\lambda$$ from 0.01, 0.02, 0.03 and so on up to 0.95. Then they supposedly fit some natural cubic spline to that and extrapolate that to $$\lambda=1$$.

In their own words:

Consider Fig. 3, where we have plotted $$\hat π_0(λ)$$ versus λ for λ = 0, 0.01, 0.02,..., 0.95. By fitting a natural cubic spline to these data (solid line), we have estimated the overall trend of $$\hat π_0(λ)$$ as λ increases. We purposely set the degrees of freedom of the natural cubic spline to 3; this means we limit its curvature to be like a quadratic function, which is suitable for our purposes. It can be seen from Fig. 3 that the natural cubic spline fits the points quite well.

And this is Fig. 3:

Can anybody make sense of this description of the "natural cubic spline fit"? AFAIK the spline should interpolate the data, which it clearly doesn't. It looks more like just one section of a spline, but given natural boundary on the spline, that could only be a straight line IMHO, which it also isn't. I also don't understand the sentence about setting the degrees of freedom and how that "limits the curvature to a quadratic function". The curvature of a cubic spline is related to it's second derivative, which is obviously constant, and not quadratic.

I would shrug that off, if it weren't for those 10000+ citations and the authors being well respected in the field. Maybe I overlooked something? Anyway, the way the method is presented, I can't even implement it as I don't have a clue what's meant here.

REFERENCE

Storey, John D., and Robert Tibshirani. "Statistical significance for genomewide studies." Proceedings of the National Academy of Sciences 100.16 (2003): 9440-9445.

I found their code on the Wayback machine and they used the smooth.spline-function in R. The paper points to http://genomine.org/qvalue/results.html for code and data,which is defunct, but can still be found in a snapshot from 2004 which redirects to http://faculty.washington.edu/~jstorey, also defunct , but also with a snapshot, so here are the code and data:

https://web.archive.org/web/20040810001627/http://faculty.washington.edu/%7Ejstorey/qvalue/results.html

And here's the code for fig. 3 at the bottom of the code file:

#Figure 3
library(modreg)
lam <- seq(0,0.95,0.01)
pi0 <- rep(0,length(lam))
for(i in 1:length(lam)) {
pi0[i] <- mean(p>lam[i])/(1-lam[i])
}
spi0 <- smooth.spline(lam,pi0,df=3,w=(1-lam))
plot(lam,pi0,xlab=expression(lambda),ylab=expression(hat(pi)[0](lambda)))
lines(spi0)


library(modreg) is now part of the base R: https://stat.ethz.ch/pipermail/bioconductor/2010-June/034197.html

So the trick for making the fit so stable at the end seems to be, that they used $$1-\lambda$$ as weights. Also natural splines might have been a misnomer.

I'm not sure what we have learned here, except that the Wayback Machine is very useful an maybe worth a donation https://archive.org/

• I think you found the real problem (+1): that probably wasn't actually a natural cubic spline. With the smoothing spline, the df has to do with the amount of smoothing rather than the number of knots. Storey's current implementation in the Bioconductor qvalue package still uses smooth.spline() with 3 degrees of freedom as the default.
– EdM
Commented Jun 19, 2023 at 20:36
• Nice work digging that up! Commented Jun 19, 2023 at 20:54
• Great! As so often, the truth is only in the code (and also doesn't completely fit the description in the paper). Thank you so much! Commented Jun 19, 2023 at 21:25
• The answer basically gives a broader lesson to optimally use the Web Archive and for that only, +1. Commented Jun 20, 2023 at 4:14

Edit: In light of Lukas Lohse's answer (which I think should be the accepted one!), my original answer below is misleading.

Personally I learned about splines from Tibshirani's books, where he introduces "natural splines" in the context of regression splines first, and then treats smoothing splines separately. So I was assuming he used "natural cubic spline with df=3" here to mean "natural cubic regression spline, with 2 internal knots." Since the knots were not specified, I assumed they were probably chosen as quantiles by default, ie knots around the .33 and .67 quantiles of $$X$$.

But based on the code, it seems more likely that this paper meant "natural cubic smoothing spline, with 3 effective df," i.e. with $$\lambda$$ chosen so that the trace of the smoothing matrix is 3.

It may help to distinguish some terminology -- at least as it's used in Rob Tibshirani's books and papers. You seem to be familiar with smoothing splines (with a knot at every x-value), as in @F.Tusell's answer. But there are also regression splines (with typically just a few knots, often chosen at regular quantiles of x or sometimes chosen based on subject-matter knowledge).

I do wish these methods had been given names that made them easier to distinguish!

A natural cubic spline is specifically a regression spline with the added constraint that it must be linear beyond the outermost knots. Due to these constraints, a natural cubic spline with df=3 should have 2 internal knots (though that might possibly be off by one, depending on the software they used and whether the intercept is counted or not for df purposes).

See Ch 7.4-7.5 of An Intro to Stat Learning https://www.statlearning.com/ for more on regression splines vs smoothing splines.

• Well, if someone says "spline" to me, I hear "interpolating spline", which is IMHO the "default". Of course you can do regression with splines, no problem with that, but then you need to define the knots, which they don't do in the paper. Commented Jun 19, 2023 at 21:31
• @Elmar Zander, The fact that in Fig. 3 of your original question the spline did not interpolate the points, was what prompted me to assume that a smoothing spline was being used. As far as I have seen, when no mention of the knots is made in a paper, it is usually the case that all points are used as knots, with a penalty taking care of the smoothing. Commented Jun 20, 2023 at 5:15
• @ElmarZander I agree that Storey & Tibshirani should have been more careful in how they wrote up this paper! There seem to be different traditions: Interpolating splines or smoothing splines might be the "default" if you learned about splines from other books. But regression splines are the "default" in Elements of Stat Learning and Intro to Stat Learning (both co-authored by Tibshirani), so I originally assumed that this "natural cubic spline with df=3" must be a regression spline with 2 internal knots. Yet that's clearly not true in the code that Lukas Lohse unearthed... Commented Jun 20, 2023 at 12:49
• @GavinSimpson Standard regression splines do require specifying the knots. Wood seems to get around this "knot placement" problem, by some adaptive basis construction scheme (thanks for the link, BTW, quite interesting). But this is far from trivial, and would definitely need to be mentioned in the paper. Commented Jun 20, 2023 at 22:09
• @ElmarZander Sorry, I wasn't suggesting that something like Wood (2003) was being used in the paper, just a response to your comment "Of course you can do regression with splines, no problem with that, but then you need to define the knots". The terminology around splines, it's (mis)use online in fora like this, & other general comments has made trying to learn splines incredibly difficult.Comments here are an eg of this. All that is missing in the paper is the word "smoothing" (they used a natural cubic smoothing spline) plus a pointer that they implemented this via smooth.spline() in R Commented Jun 21, 2023 at 6:25

I have not read the paper you ask about, but may be I can help with some of the questions that puzzle you.

A cubic spline, natural or not, can be made to interpolate the data. This is usually of no statistical interest. What is of interest is a smoothing spline, the solution $$g(x)$$ to the problem $$\min\left\{\sum_{i=1}^n(y_i - g(x_i))^2 + \lambda\int [g''(x)]^2dx\right\}$$ where $$(x_i,y_i)$$, $$i=1,\ldots,n$$, are the points to be fitted, $$g(x)$$ is a spline with knots at $$x_1,\ldots,x_n$$, and $$\lambda$$ a smoothing parameter. If $$\lambda=0$$ you have an interpolating spline, but as $$\lambda$$ increases, it penalizes more and more the curvature and you may have something as shown in your Fig. 3.

The book Hastie and Tibshirani(1990) Generalized Additive Models, Chapman & Hall, might be, among many other sources, a good introduction to smoothing splines.

• Well, I know smoothing splines from Green and Silverman, but Storey and Tibshirani don't talk about "smoothing" anywhere, and then they'd have another parameter $\lambda$, which they would have to determine, which is also not trivial... Commented Jun 19, 2023 at 15:35
• I must respectfully disagree about such splines as being of "no statistical interest," at least with respect to the application is question. The problem involves an extrapolation of a presumably smooth underlying function through a region of high variability at high values of $\lambda$, to estimate the limit as $\lambda=1$. As much as I appreciate the benefits of a smoothing spline for some applications as in many GAM models, for this application the natural spline fit with its restriction to a linear fit beyond the outermost knots makes a lot of sense.
– EdM
Commented Jun 19, 2023 at 16:28
• I think that we will both be amused that the implementation in Storey's related qvalue() function in Bioconductor uses a default smooth.spline() with 3 degrees of freedom (trace of the smother matrix) in its pi0est() function, rather than a natural spline as reported in the PNAS paper.
– EdM
Commented Jun 19, 2023 at 20:24
• @LukasLohse great minds think along like channels :-)
– EdM
Commented Jun 19, 2023 at 20:37
• @EdM smooth.spline() is fitting a natural cubic smoothing spline in the sense that the second derivative is 0 at the boundaries (as one can confirm using the cars example from ?predict.smooth.spline and running predict(cars.spl, c(min(dist), max(dist)), deriv = 2). Commented Jun 20, 2023 at 17:22

The cited paper helped popularize the false discovery rate (FDR) as a way to evaluate the significance of results on thousands of features, as are examined in genome-wide studies. In that paper, Storey and Tibshirani modified the original Benjamini-Hochberg method "to improve upon a limitation of the FDR, namely that the FDR is not defined when there are no positive results."

Figure 3 shows a way to estimate the actual fraction of true negatives among the results, $$\pi_0$$.

The p-values for true negatives have a uniform distribution over [0,1]. If you look only at cases with high p values and thus those most likely to be true negatives, (1) the distribution of p values is relatively flat (see Figure 1 of the paper) and (2) the limiting height of the density distribution of those p values, at high p values, provides an estimate of $$\pi_0$$.

Beyond any p-value cutoff $$\lambda$$ among $$m$$ total p values, Storey and Tibshirani write the corresponding estimate $$\hat \pi_0(\lambda)$$ based on the number of values greater than the cutoff as:

$$\hat \pi_0(\lambda)=\frac{\# (p_i >\lambda;i=1, \dots , m)}{m (1-\lambda)} .$$

What you want is the estimate in the limit of a p-value cutoff ($$\lambda$$) of 1. Figure 3 shows the plot of values versus the cutoff for a particular data set. At very low p-value cutoffs the $$\pi_0(\lambda)$$ values will necessarily be high, as many true-positive cases will be included. With few cases having p values close to 1, the actual data points at the right of the plot are noisy.

The natural spline fit is just a simple attempt at extrapolation out to the limiting p value of 1. On a quick re-reading I didn't see where they specified the 4 knot positions for the 3 degrees of freedom used in the fit. It's possible that, as with the default of the ns() function in the R splines package, the outermost knots were set to the upper and lower limits of the $$\lambda$$ values (0, 0.95) so that you would not see much of a linear extrapolation at all.

In terms of ""limit[ing] the curvature to a quadratic function," that might just be an infelicity of expression. I need to think on that a bit more. In the meantime, just recognize that what's shown is a simple type of smoothing of the data to facilitate the extrapolation to the limiting value as $$\lambda \to 1$$.

• Thanks for your extensive answer. I know what they do in the rest of the paper and I really appreciate it. However, that part with the spline fit is really weak and this is my only question: given $\hat\pi_0(\lambda)$ for $\lambda$ in $[0, 0.01, \dots, 0.95]$, how exactly do they compute the spline fit and thus $\hat\pi_0(1)$. This is still completely unclear. Commented Jun 19, 2023 at 18:25
• BTW Even in the appendix, where they give the "complete" algorithm, they just write for this step "3. Let $\hat f$ be the natural cubic spline with 3 df of $\hat\pi_0(\lambda)$ on $\lambda$." and "4. Set the estimate of $\pi_0$ to be $\hat\pi_0=\hat f(1)$". I had quite some exposure to splines, but I have no clue what "natural cubic spline with 3 degrees of freedom" is supposed to be... Commented Jun 19, 2023 at 18:29
• @ElmarZander A "natural cubic spline" is a restricted cubic spline; see Wikipedia. Section 2.4.5 of Harrell's Regression Modeling Strategies notes that such a spline with $k$ knots involves estimating $k-1$ coefficients, $k-1$ is thus what I think about as the number of degrees of freedom (df) associated with the spline. I didn't see in the paper where they set the knot positions, however, or if they have a different interpretaiton of df as in the answer from @civilstat.
– EdM
Commented Jun 19, 2023 at 20:18
• @ElmarZander it's also possible that the paper didn't use a natural spline. Storey's related qvalue package uses a smoothing spline (as in the answer from @F.Tussell) with 3 degrees of freedom as the default in its underlying pi0est() function, not a restricted cubic regression spline. In that case, the df is the effective number of degrees of freedom, equal to the trace of the smoother matrix.
– EdM
Commented Jun 19, 2023 at 20:30
• I guess you're right here. Still, I would hold that in an algorithm in a paper, each step should uniquely specify a mathematical object, that is computed in that given step, and not rely on some specific implementation, which isn't even mentioned. Commented Jun 19, 2023 at 21:42