# Clarification on how to perform nested cross-validation

I'm currently trying to get lifetables for our population based on death and population counts. My first idea was to follow this paper methodology but after some discussions we are planning to use different models, such as GLMMs and GAMs and explore the inclusion of different knots and families.

The author suggested to use train and test sets or cross-validation. I've been reading (great resources and questions on this site!) and I would like to know if the methodology I'm using is appropriate as I'm not entirely sure.

Here's what I have in mind:

1. Shuffle the data
2. Divide the data into training (80%) and test set (20%).
3. Test the different parameters in the training data, and compare them using the test data.
4. Once I have the chosen family and knots, shuffle the data again.
5. Do n iterations (20?) of 5-fold cross validation (size dataset=3440 observations) of one model to get the performance of the methodology (The idea is that the model chosen in 4 is the one able to better predict the data).

Now, based on the figure I've uploaded, here are my questions and confusion:

Should I use a nested cross-validation and in the inner loop test the model parameters (as the figure)? Each k-fold correspond to one set of parameters or I need to repeat a nested CV for each set of parameters? The theoretical number of knots may be too low or I may need to add a smooth factor to one of the variables, etc. I'm not sure if I can play with these parameters in the inner loop/s.

This brings me to my other question. Do I need to do one k-fold cross validation per method? If I use cross-validation for a GAM with a poisson family and a GAM with a negative binomial I may be on a fold where that model is able to better predict the data, where in another fold the same model may not perform so well, so instance the set of knots could be inappropriate.

Another thing that I'm unsure of is that one variable is calendar year. Do I need special care when shuffling the data? Perhaps make sure each fold has information about all the years or something else? I also thought of using a GLMM or GAMM, perhaps that would be enough?

My other question relates to the output parameters and how to evaluate the cross validation. Since my aim is to predict I was planning to look at the residuals of each model, AIC, and plot the predicted and observed values to see which method performed better. Should I use different parameters?

I appreciate any help or guidance as I just started reading about cross validation and have no experience.

If I understood your answer, basically nested cross validation would be something along those lines:

for each fold in the outer loop
for each parameter in the grid
run each parameter in the inner loop
evaluate each parameter in each k-fold
chose best parameter based on results from inner loop
run best parameter in each k-fold of the outer loop


The parameters I would like to test are:

1. Different set of knots for age (some fixed at the boundaries but I would like to test a few others for middle ages).

The reason to choose knots is due to extreme values at the boundaries, as there's a large mortality for the first point (infant mortality) and a large mortality above 85 (since data above 85 is grouped).

I just got Simon Wood's book so I can understand more about GAMs. I thought about using cubic splines but thin splines could be interesting. It would be interesting to look at the difference between thin plate splines and the use of predefined knots. If thin plate splines fit the data well it would be a great option.

1. GLM with a poisson and NB family and a GAM with poisson and NB family.
2. I would like to check the inclusion of a random effect for region.
3. I would also like to test whether an interaction between gender and age, and region and age improves the model.

I wasn't thinking about using ti for interactions but instead something along those lines:

gam(deads ~ region*age + gender*age + calendar_year + ti(age,k=7,bs="cr") + offset(log(pop)), knots=list(age=list_knots), family=poisson, data=data, method="REML")


But I also thought of trying something like this:

gam(deads ~ s(age,bs="fs",by=region) + s(age,bs="fs",by=gender) + calendar_year + ti(age,k=7,bs="cr") + offset(log(pop)), knots=list(age=list_knots), family=poisson, data=data, method="REML")


Obviously, I need to read more to check what makes sense and is more appropriate.

1. Finally, I was wondering whether I should add some factor smoothing to the factor variable.

Now, perhaps I should have a clearer idea of the final model. Let me hear your thoughts. But having all of these parameters to test was/is some part of my confusion about cross-validation in the inner loop.

I gave a general answer to this question, and here is what applies to your question:

Train and Validation Split:

First split the input into train and validation; but I'd also take the domain knowledge into account. In your case, I would take that year parameter into account and take the last few years of the data (not sure how many years your have, let say 2 out of 10, if you have 10) and assume that portion of the data as your validation set.

Nested Cross Validation and Parameter Search:

Now you can do what you explain in your diagram. Assume you have a method, which takes the input data, and the parameters (e.g. a parameter defining to use a GAM with a poisson family or a GAM with a negative binomial), and fit the corresponding model on the data. Let's call the set of all these parameters you're considering, a parameter grid.

Now for each of those outer folds, you do a whole grid search using the inner folds, to get a score for each parameter set. Then train your model using that best parameter set on the whole data given to you in the inner loop, and get its performance on the test portion of the outer loop.

Assume your parameter grid has 3 values in total (e.g. a GAM with a poisson family or a GAM with a negative binomial, and a [not regularized] linear model), and there's no other parameter involved. Then you'd do these many trainings:

$5 [\text{outer loop}] \times \left(3[\text{parameter grid}] \times 4[\text{inner loop}] + 1[\text{best parameters}]\right)$

Talking in code, here's how it'd look like:

parameter_grid = {'param1: ['binomial', 'poisson'],
'smoothing': ['yes', 'no']}
scores = []
for train, test in outer_folds:
model = GridSearchCV(estimator=my_custom_model,
param_grid=parameter_grid,
refit=True,
cv=4)
model.fit(train.X, train.y)
scores.append(model.score(test.X, test.y))
score = mean(scores)


For simplicity, I'm diverging from the actual API, so the above code is more like a psuedocode, but it gives you the idea.

This gives you an idea about how your parameter grid would perform on your data. Then you may think you'd like to add a regularization parameter to your linear model, or exclude one of your GAMs, etc. You do all the manipulations on your parameter set at this stage.

Final Evaluation:

Once you're done finding a parameter grid you're comfortable with, you'd then apply that on your whole train data. You can do a grid search on your whole train data with a normal 5 fold cross validation WITHOUT manipulating your parameter grid to find the pest parameter set, train a model with those parameters on your whole train data, and then get its performance on your validation set. That result would be your final performance and if you want your results to be as valid as thy can be, you should not go back to optimize any parameters at this point.

To clarify the parameter search at this stage, I'm getting help from scikit-learn API in Python:

parameter_grid = {'param1: ['binomial', 'poisson'],
'smoothing': ['yes', 'no']}
model = GridSearchCV(estimator=my_custom_model,
param_grid=parameter_grid,
refit=True,
cv=5)
model.fit(X_train, y_train)
model.score(X_validation, y_validation)


The above code does (model.fit(...)) a 5 fold cross validation (cv=5) on your training data, fits the best model on the whole data (refit=True), and finally gives your the score on the validation set (model.score(...)).

Deciding on what to out in your parameter_grid in this stage is what you do in the previous stage. You can include/exclude all the parameters you mention in there and experiment and evaluate. Once you're certain about your choice of parameter grid, then you move on to the validation stage.

• Thank you for your response! I'm a bit confused by this sentence You can do a grid search on your whole train data with a normal 5 fold cross validation WITHOUT manipulating your parameter grid to find the pest parameter set, train a model with those parameters on your whole train data, and then get its performance on your validation set. I've updated my question with a sketch of an algorithm and all the parameters I would like to test just to make sure I understood how to proceed. I haven't started yet as I'm trying to get the process as accurate as I can. – psoares Apr 4 '18 at 10:13
• I've updated the answer. – adrin Apr 4 '18 at 10:31
• Thank you but I'm still a bit confused. The purpose of the outer loop is only to evaluate the final parameters (or parameters?). We only use the test set of the outer loop once to validate the final model, is that it? All of the parameter tuning should be done in the inner loop, right? Now, the outer loop has 5-fold so there's 5 inner loops of 4-fold. Is the tunning done in all of these folds? And once the parameters are chosen, we repeat the nested CV? I hope this was not too confusing =/ sorry about that! – psoares Apr 4 '18 at 10:53
• It kinda is confusing. I don't understand the question. – adrin Apr 4 '18 at 10:55
• Sorry about that. When I'm confused it's even more difficult to explain the random thoughts in my mind :) Let me try in steps. In the image I uploaded the outer loop has 5-folds and the inner loop has 4-folds, so basically the number of folds for tunning parameters is 4*5=20. Correct? – psoares Apr 4 '18 at 11:05

What are you wanting to achieve by selecting knots? As a means of determining the complexity of the smooth terms in the model?

Generally there's not much to be gained from choosing knots in modern GAM theory; thin plate splines for example are the default in mgcv and have a knot at each unique data location (the subsequent basis is then eigen decomposed to give a lower rank approximation), and smoothness selection via ML or REML strategies typically works well enough to avoid over fitting. If not using thin plate splines, just spread the $k-1$ cubic regression spline knots evenly over the data interval or at $k-1$ evenly spaced quantiles of the data. If you are using other bases, such a p splines, it may not even make sense to choose knot locations beyond evenly spacing them over the interval of the data.

The key user choice is to select a sufficiently large basis dimension for each smooth, but there are tools to test for that. The penalty on wiggliness typically does the rest assuming the initial basis expansion is sufficient large/rich enough to either include the true function or a close approximation to it.

The GLM is a special case of the GAM where all terms are perfectly smooth (linear). Hence you could fit the GAM with smoothness selection and if linear effects are appropriate the smoothness penalty should become large (greater penalty on wiggliness). Alternatively, you could formulate the GAM as follows

gam(y ~ x + s(x, bs = "tp", m = c(2, 0)))


where the smooth term is specifically testing for a smooth effect of x in addition to the linear effect already included.

As for choosing the family, model diagnostics could motivate this also, as can the estimate for the parameter informing the extra variation of the Negative binomial over the Poisson distribution.

For the random effect, yuo can fit this as a random effect smoother in gam() or bam() via bs = "re" and then use the summary() output to decide if this smooth is required — the test there has good properties for random effect terms. I also believe AIC as implemented in recent mgcv versions is suitable for use with simple random effect terms such as this.

Interactions can be tested using inference via a ANOVA-like decomposition into "main effects" and "interaction" via ti() terms:

y ~ ti(x) + ti(z) + ti(x, z)


where one would focus on the ti(x, z) to determine whether the interaction was required or not.

As for the factor-smooth question; fit both models using all the data (one with and one without that factor-level smooth) and compare the model fits via AIC and the Wald-like tests reported in the summary() output.

My general point here is that it may be better when fitting GAMs, which learn from the data, to use all the data to inform the estimates of the smooths. Then use statistical inference to decide which terms are needed.

Or just fit the most complex model and include extra penalties (via select = TRUE with method = "REML" or method = "ML") to do the selection for you.

• Thank you for your exhaustive and helpful comment! I have edited my answer to address some of your comments. Could you please let me know your thoughts? The other points I agree with you and I’m planning to do those. I already have an idea of the model and family that could work based on visual inspection of the data and theoretical knowledge but I would still like to see whether interactions and smooth terms are needed. You make a valid point regarding using all the data. Do you have any literature for the use or not of GAM with CV? – psoares Apr 5 '18 at 8:36

• About the nested cross validation I've come for myself to a description that I find makes it easy to get the calculations right while (by) not describing the situation as nested at all.
The basic idea is that I need two things: a way of training my model as function of my data that takes care of all autotuning/selection/hyperparameter optimization train_tuned_model (training_data) and a validation/verification of the model obtained by train_tuned_model (training_data). The latter may be done by a simple straightforward cross validation. And internally, train_tuned_model (training_data) may use another out-of-bootstrap or cross validation or whatever heuristic I deem suitable.
This way, I still get the correct calculations the same way I get them with nested cross validation, but conceptually, I can describe the whole thing by easier and commonly available building blocks. Incidentally, this also answers the question, whether the final model trained from the whole data set should use the optimized hyperparameter set found inside the nested cross validation or whether hyperparameters should be optimized again (my answer: yes).
• Regarding whether you do:

a) outer CV loop - hyperparameter loop - inner CV loop or
b) outer CV loop - inner CV loop - hyperparameter loop:

Both are sensible. b) is a paired design and a) is its unpaired analogon. So as the paired design is possible, this may give a little bit more power to your optimization.

• One point that has not yet been mentioned: the success of the cross validations (both of them) depends crucially on your splits generating independent subsets of your data. Your data has a much more complicated structure than a simple "rows are [assumed to be] independent". And you know a lot about this structure. Do take this structure into account:

• "calendar year. Do I need special care when shuffling the data?" As always: this depends.

• There are data/models where year is de facto a random factor (This year's harvest data doesn't help at all predicting next year's harvest, and the model will be used for unknown future years). In that case, you need to make sure that test years are independent of training years.
• Other data/models may be longitudinal in nature (i.e. knowing the last years helps predicting the next year(s)) and the production use scenario always has the last $n$ years available. In this case you need to split so that training is always the last n years before test year. (Same principle may needed for spatial data)
• With calender year I'm running out of fantasy formulating a scenario where year is a normal fixed factor. Maybe an approximation would be that you find that it is best to recognise whether the year was cold-wet, wet-hot or dry-hot. You do not need to take any special precaution for fixed factors.
• In general, splitting needs to make sure random factor levels in test data do not occur in the training data.

• Fixed factors to not need special care in the splitting,
• but you need to be careful and make sure the factor is really fixed wrt. your prediction scenario. Fixed factors in data collection may be random factors in prediction scenarios (see above).
• Factors can be quite deeply nested: in that case, correct splitting at the uppermost random factor will automatically lead to correct splitting for the lower levels.
• If random factors are crossed, splitting into independent levels for training and testing implies that some data is excluded from both training and testing for this split.
• Thank you your answer! Optimisation/tuning and verification are good terms and (to me) helps to avoid some confusion. We only have 5 years of data (2011 to 2015), so I don’t expect a large variability among the number of deaths per year. There’s an increase in mortality for older ages, but there’s also an increase of the older population, which is used as an offset. If I understood correctly you’re suggesting keeping 2015 as test set and the remaining years as training. We also have information about the region, which also depends on the population (large regions more deaths). – psoares Apr 6 '18 at 7:58
• Can I clarify a few details? How do I perform the k-fold CV if the training is always the last n years before test data? If the test data is 2012 then the training should be 2011,2013 until 2015. Or should I perform CV in a ladder form? Train on 2011 and test on 2012, next k-fold, train on 2011 and 2012 and test on 2013, and so on. – psoares Apr 6 '18 at 16:19
• @pavid, Bob Hyndman has a blog post about this: robjhyndman.com/hyndsight/tscv (I think 5 years is probably just not enough to train a longitudinal model.) – cbeleites unhappy with SX Apr 6 '18 at 19:18