Linear model selection after bootstrapping without overfitting I am trying to develop a model between a 19 year record of climate data and a 19 year record of ice-off dates on rivers. The two variables are linearly correlated. The goal is to build a linear model so that we can use the climate data to predict the ice-off dates in future years when we don't have ice data but do have climate data.
What I have done thus far is bootstrapping: I randomly select 14 years as training data and the remaining 5 years as testing data. I build the linear model on the 14 year training dataset, then apply it to the remaining 5 years, and evaluate the model performance using the nash-sutcliffe coefficient (https://en.wikipedia.org/wiki/Nash–Sutcliffe_model_efficiency_coefficient#targetText=The%20Nash–Sutcliffe%20model%20efficiency,Qm%20is%20modeled%20discharge.). I then repeat that 1000 more times, randomly sampling the 14 years of training data each time.
Now that I have done this, I want to pick the best model of the bunch. Should I take the model with the median nash-sutcliff coefficient, or the one with the best nash-sutcliff coefficient? What is the best next step here that avoids overfitting?
I'm a statistics beginner, so your help is greatly appreciated!
 A: AFAIK, bootstrapping is for getting the standard error estimates. If you want to validate models, look at N-fold cross-validation, or leave-one-out cross-validation (jackknife).
EDIT: since you said that years are not serially correlated, this simplifies your situation a lot! You can then treat the years as if they were just data points, and select the testing sets totally randomly. You can do e.g. 3-fold cross-validation by splitting the data in 3 randomly selected set of years, run the model 3 times (see how N-fold cross-validation works), and then put together all the 3 independently predicted testing sets from those 3 models and evaluate them with the Nash-Sutcliffe coefficient, which seems to be a good measure of the model efficiency.
Chose the model with the best coefficient. Cross-validation assures  that you do not overfit in those single models tested. However, since during the model selection procedure itself you use the whole dataset, there is a risk that you can overfit in that model selection procedure itself, if you overdo it, as fittingly pictured in the chart below:

So, be careful not to select among too many model variants. Another way could be to put apart yet another validation set to monitor how are you doing in the model selection procedure itself :)
Also see the answer here and this reference:
Cawley and Talbot "On Over-fitting in Model Selection and Subsequent Selection Bias in Performance Evaluation", JMLR, vol. 11, pp 2079−2107, 2010
A: In the end, you pick a model which is estimated on the whole data. This is not a model which is a weighted combination of inferior models estimated on only a subset of the data (which you seem to be implying). 
In the final, "optimal" model you need to have only statistically significant terms. For example, if the slope coefficient is non-significant according to the bootstrap-based p-value, the predictor has to be dropped from the model. This means that only the intercept remains provided it is statistically significant.
Finally, you are not performing bootstrap correctly. Each bootstrap sample must be created by sampling 19 data points with replacement. You, on the other hand, are sampling 14 data points without replacement. What you are doing is a randomized cross-validation, not bootstrap. Cross-validation is suitable for comparing different models but it is not suitable for estimating standard errors (and p-values) of the coefficients. 
