Aggregation of multiple simulations to stabilize p-values

[florianhartig]

When re-simulating residuals for fixed data and integer responses, p-values may spread within a certain range, due to the randomization procedure to smoothen out the integer values.

This phenomenon was discussed in https://github.com/florianhartig/DHARMa/issues/37

Particularly in low-data situation, this can result in p-values being quite variable. In this case, it might be desirable to have an option to obtain an aggregate p-value from multiple simulations.

Question is

a) how to best do this

b) what the properties of the resulting p-values are in terms of their distribution / type I error / power

[florianhartig]

Note: gamlss seems to suggest a similar procedure in https://www.rdocumentation.org/packages/gamlss/versions/5.1-4/topics/rqres.plot

[florianhartig]

D&S comment that they recommend several plots. I think what's easier to do is to do several simulations, and then plot with densities, and adjust weights in all tests or correct for multiple testing.