Using the standard deviation in Cross Validation I'm running a Grid Search to find the optimal parameters for xgboost via sklearn.
I can see that the grid search picks the set of parameters with lowest mean MSE.
The problem is that upon inspecting the standard deviations they are in the range of the mean, which suggests that there is not statistical significance between the choice of parameters.
Can a means difference test be ran on these means? How do we make sure there is a statistically significant difference between the means?
And I find it weird that these packages work picking the lowest mean without paying attention to the standard deviation. Do you know any package that looks into this issue?
 A: 
"How do we make sure there is a statistically significant difference
  between the means?...these packages work picking the lowest mean
  without paying attention to the standard deviation..."

Searches for parameter values are looking for an optimum solution, based on specified criteria, in a particular learning/modeling task. There is no assurance that the parameter values at the optimum solution can be distinguished statistically from other sets of parameter values in the search space. One could argue that with finite data sets and smooth optimization criteria it would be impossible to find "significant" differences between the identified optimum values and sets of parameter values that are sufficiently close to that optimum.
What matters is how well the model based on the chosen parameter values performs for the purpose that you have in mind. You will be better served by focusing more on the significance of the model performance overall (via cross-validation, bootstrapping, etc.) and paying less attention to whether your optimum set of parameter values has a "statistically significant" difference from some other set of values.
A: Let me add some points to @EdM's excellent answer (which are too long for a comment).

 And I find it weird that these packages work picking the lowest mean without paying attention to the standard deviation. 



*

*I totally agree with you (OP) knowing uncertainty could (IMHO should) be used to stop the optimization process: no point in trying to tie down the optimum within a zone where we cannot distinguish performance/target function values. 

*From a modeling point of view, of course, EdM's "sufficiently close to the optimum" should hopefully be a comfortably large zone so even if we have uncertainty about the true optimum we can assume that our opimized model isn't doing too much worse than the model fit using the unknown true optimum of hyperparameters. 

*Of course, we can also turn around the question a bit. We do have (at least) two sources of variance to our performance evaluation: model instability and the finite number of cases we test. Model instability is IMHO something that crucially needs to be addressed during optimization (can be done e.g. by the 1-sd-rule mentioned above). Finite number of test cases is something we can influence rather directly, and IMHO we should check at the very least that the sample size we have allows to distinguish statistically what we consider practically significant differences in performance wrt. the application at hand.  

Do you know any package that looks into this issue?

No, but I've been working on and off (as time permits/optimizations are needed) towards this. In R, though. I'd be happy to join forces if you like.


*

*One variable at a time is not going to work unless you know the hyperparameter in question does not interact with other of your hyperparameters. Which is rare, but occasionally the case (random forest number of trees comes to my mind). 

