When we plot data and then use nonlinear transformations in a regression model are we data-snooping?

Question

I've been reading up on data snooping, and how it can mean the in-sample error does not provide a good approximation of the out-of-sample error.

Suppose we are given a data set $(x_1,y_1),(x_2,y_2),...,(x_n,y_n)$, which we plot, and observe what appears to be a quadratic relationship between the variables. So we make the assumption that $$ y_i = \beta_0 + \beta_1 x_i + \beta_2 x_i^2 + \varepsilon, $$ where $\varepsilon$ is a variable representing noise.

Isn't this data snooping? We have let the data affect our model. So what implications does this have for the coefficients $\beta_0,\beta_1,\beta_2$ that we find; can they be considered reliable for making future predictions with different input variables?

I ask because there are countless notes/articles/books/etc.. on regression where they recommend looking at the data and then choosing a model that looks like it will fit well with the data. For example, here the author has some data, tries a linear model, and upon finding it unsatisfactory, he moves to a quadratic model which better fits the data. Similarly, here, people are discussing log transformations and the original poster is given the following advice:

If there is no theory to guide you, graphical exploration of the relationship between the variables, or looking at fitted vs observed plots both ways will tell you which model is appropriate.

So when we base our model on an observation of the plotted data, is this data snooping or not? If it isn't, then could someone give an explanation why this isn't data snooping?

If it is data snooping, then:

1. What are the consequences of this on the out-of-sample performance?
2. What should we do to avoid/overcome the data snooping issue in a regression model so that we will have good out-of-sample performance?

Answer 1

There is a way to estimate the consequences for out-of-sample performance, provided that the decision-making process in the modeling can be adequately turned into an automated or semi-automated process. That's to repeat the entire modeling process on multiple bootstrap re-samples of the data set. That's about as close as you can get to estimating out-of-sample performance of the modeling process.

Recall the bootstrap principle.

The basic idea of bootstrapping is that inference about a population from sample data (sample → population) can be modelled by resampling the sample data and performing inference about a sample from resampled data (resampled → sample). As the population is unknown, the true error in a sample statistic against its population value is unknown. In bootstrap-resamples, the 'population' is in fact the sample, and this is known; hence the quality of inference of the 'true' sample from resampled data (resampled → sample) is measurable.

Following that principle, if you repeat the full model building process on multiple bootstrap re-samples of the data, then test each resulting model's performance on the full data set, you have a reasonable estimate of generalizability in terms of how well your modeling process on the full data set might apply to the original population. So, in your example, if there were some quantitative criterion for deciding that quadratic rather than linear modeling of the predictor is to be preferred, then you use that criterion along with all other steps of the modeling on each re-sample.

It's obviously best to avoid such data snooping. There's no harm in looking at things like distributions of predictors or outcomes on their own. You can look at associations among predictors, with a view toward combining related predictors into single summary measures. You can use knowledge of the subject matter as a guide. For example, if your outcome is strictly positive and has a measurement error that is known to be proportional to the measured value, a log transform makes good sense on theoretical grounds. Those approaches can lead to data transformations that aren't contaminated by looking at predictor-outcome relationships.

Another useful approach is to start with a highly flexible model (provided the model isn't at risk of overfitting), and pulling back from that toward a more parsimonious model. For example, with a continuous predictor you could start with a spline fit having multiple knots, then do an analysis of variance of nested models having progressively fewer knots to determine how few knots (down to even a simple linear term) can provide statistically indistinguishable results.

Frank Harrell's course notes and book provide detailed guidance for ways to model reliably without data snooping. The above process for validating the modeling approach can also be valuable if you build a model without snooping.

Answer 2

Here's a basic answer from a machine-learning perspective.

The more complex and large the model class you consider, the better you will be able to fit any dataset, but the less confidence you can have in out-of-sample performance. In other words, the more likely you are to overfit to your sample.

In data-snooping, one is engaging in a search through a possibly-very-large-and-flexible model space. So the chance of finding a model that overfits becomes more likely.

We can prove this doesn't happen (with high probability, under conditions) if the model space is limited enough, compared to the dataset size.

...

So the distinction between data-snooping and principled investigation can be as fine as: the space of models that, a priori, one is willing to consider.

For example, suppose that author finds no quadratic fit, so they move on to cubics, quartics, ..., and eventually they find a degree-27 polynomial that is a good fit, and claim this truly models the data-generating process. We would be very skeptical. Similarly if they try log-transforming arbitrary subsets of the variables until a fit occurs.

On the other hand, suppose the plan is to give up after cubics and say that the process is not explainable in this way. The space of degree-at-most-3 polynomials is quite restricted and structured, so if a cubic fit is indeed discovered, we can be pretty confident that it is not a coincidence.

...

Therefore, one way to generally prevent "false discovery", as we often call it, is to limit oneself a priori to a certain restricted set of models. This is analogous to pre-registering hypotheses in an experimental work.

In regression, the model space is already quite restricted, so I think one would have to try a lot of different tricks before being at risk of discovering a spurious relationship, unless the dataset is small.

Answer 3

Here is an answer from a physics perspective. If you are doing excessive "fitting," then you might be data snooping. However, if you are "modeling" in the way we mean in physics, then you are actually doing what you are supposed to do.

If you're response variable is decibels and your explanatory variables are things like power input and material properties, then if you didn't model in log space, you would be doing it wrong. This could be an exponential model, or a log transform.

Many natural phenomena result in not-normal distributions. In these cases, you should either use an analysis method that allows you to incorporate that distribution structure (Poisson regression, negative binomial, log-linear, lognormal, etc.) or transform the data keeping in mind that will also be transforming the variance and covariance structure.

Even if you don't have an example from the literature backing up the use of some particular distribution that is not normal, if you can justify your claim with a minimal explanation of why that distribution might make physical sense, or through a preponderance of similarly distributed data reported in the literature, then I think you are justified in choosing that given distribution as a model.

If you do this, then you are modeling, not fitting, and therefore not data snooping.

Answer 4

Finding iteratively the best analytical model that fits data that has an error term is acceptable within the constraints nicely explained in the article you quote.

But perhaps what you are asking is what is the effectiveness of such model when you use it to predict out-of-sample data that was not used to generate the model. If it is reasonable to assume that the data generating mechanism used to calculate the model and the mechanism that generates the new data are the same, there is nothing wrong with using the model you obtained.

But you may have some justifiable scepticism about this assertion which goes to the essence of frequentist statistics. As you develop the model, you obtain the parameters that best fit the data. To get a better model you add more data. But that does not help if you add data points that you do not know whether they belong to the same data-generating mechanism used to develop the model.

Here the issue is one of belief about how likely it is for the new data point(s) to belong to the same mechanism. This takes you directly to Bayesian analysis by which you determine the probability distribution of the parameters of the model and see how this distribution changes as you add more data. For an introductory explanation of Bayesian analysis see here. For a nice explanation of Bayesian regression see here.

Answer 5

We have let the data affect our model.

Well, all models are based on data. The issue is whether the model is being constructed from training data or testing data. If you make decisions of what type of model you want to look into based on plots of the training data, that's not data snooping.

Ideally, any metrics describing the accuracy of a model should be derived from completely "clean" data: that is, data that the model generation process is not in any way dependent on. There's a tension here, as the more data you train your model on, the more accurate it can be, but that also means there is less data to validate it on.

The difference between training a model, and choosing between two models based on their validation scores is, in some sense, a matter of degree rather than kind. It can be a very large degree, however. If you're choosing between two different models, then looking at validation scores gives you at most one bit of data leakage. But as you add more and more hyperparameters, the distinction between them and regular parameters can start to blur.

As you build a model, you should gradually transition from exploration, in which you prioritize fitting your model to the training data as much as possible, to validation, where you prioritize estimating out of sample accuracy. If you want to be absolutely sure that you aren't engaging in data snooping, you should find someone to run your model on data that you have no access to.

Answer 6

Another physics perspective (see also @albalter's nice answer).
In the analysis of physics data understanding the "error" bars in measurement is of paramount importance. If you cannot account for the size of the measured errors that may reveal an unknown causative phenomenon or, more usually and sadly, some unforeseen problem with the experiment. Physicists are also concerned with recognizing systematic errors as opposed to random errors.

Fitting a model to a dataset whose values have been transformed carries the problem that the statistical distribution of the data measurement error is distorted. If the measurement errors were normally distributed they are no longer normally distributed after applying the transformation and so least squares optimization of the model may not be appropriate. Such a fit is convenient as a data summary, but the fitted function may have no underlying physical basis.

Whether this matters or not depends on your use of the fit: data summary or physics.

An outstanding example in astronomy comes from the measurement of the radial distribution of light elliptical galaxies (roundish objects with no gas) where, during the 1960's the popular fitting function was $e^{r^{1/4}}$, a "law" proposed by G. de Vaucouleurs in 1953 and still in use this century (doi: 10.1093/mnras/113.2.134).

There has been no satisfactory physical explanation for the "de Vaucouleurs law". There are several physically motivated light profiles based on assumptions about internal stellar relaxation processes and the influence of external tidal forces exerted by neighboring galaxies. The irony today is that better data has led to fitting generalizations of the de Vaucouleurs profile as in $e^{r^{1/n}}$, the Sersic profile where the value of $n$ reflects a distortion of the profile away from the de Vaucouleurs profile.

There is a fine summary of astronomical light profiles in Mark Whittle's UVa lecture notes at

http://people.virginia.edu/~dmw8f/astr5630/Topic07/Lecture_7.html

and in Amina Helmi's Kapteyn Institute lectures (with pictures!) at:

http://www.astro.rug.nl/~ahelmi/galaxies_course/class_VII-E/ellip-06.pdf