Use of splines in parameter estimation

Question

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.

Answer 1

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/

Answer 2

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.

---

Original answer (deprecated):

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.

Answer 3

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.

Answer 4

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$.