# Are there special techniques to combine estimates and standard errors from a complex discrete-event simulation?

My discrete-event simulation works as follows:

1. A sample is drawn from a simulated population
2. The simulation is run a few times for that sample
3. Each run produces k estimates and their corresponding standard errors
4. This whole process is repeated a few times (to account for the sampling variation)

In the end, I have a data frame that looks basically like this:

output <- data.frame(sample=c(1,1,1,2,2,2,3,3,3),
run  = c(1,2,3,1,2,3,1,2,3),
est  = c(1,2,2,3,1,2,3,2,1),
se   = c(2,3,2,4,1,2,3,4,3),
obs  = c(100,103,165,98,54,191,201,101,89))


Here sample refers to the index of the drawn sample at step 1, run refers to the index of the simulation run at step 2, est ist the estimate produced in the run (here, a mean value), se is the standard error of that estimate and obs is the number of obersvations used to calculate theses statistics.

My question is: How do I combine all of these to get one estimate with a corresponding confidence interval?

I realize that this is a multilevel structure. So far I used Rubin's Rule (https://bookdown.org/mwheymans/bookmi/rubins-rules.html#pooling-effect-estimates) to combine the estimates of each run inside a sample. This works, but it feels a little hacky. Especially since the equation for adjusted degrees-of-freedom (see Barnard, J., and D. B. Rubin. 1999. “Small-Sample Degrees of Freedom with Multiple Imputation.” Biometrika 86 (4): 948–55.) assumes a fixed n value, while I have different numbers of observations at each simulation run. I just used the mean value of the number of observations so far, with no justification for that.

Assuming this weird technique is correct, I would still only have managed to pool one level. Now I have a data frame like this:

output_pooled <- data.frame(sample=c(1, 2, 3),
est  = c(3, 2, 1),
se   = c(1.2, 3.4, 2.2))


(Not that all of the numbers here are completely made up and only used to illustrate my point)

Now I could use the same technique again, but I am unsure of what the obs value is now. Each sample itself is drawn without replacement (so there are no duplicates in one sample), but every sample is replaced after the full sample is drawn (so there might be duplicates accross samples).

This doesn't feel correct, but I haven't really found any other appropriate method. Papers on calculating confidence intervals for discrete-event simulations usually focus on steady-state simulations, which isn't what I am doing as each simulation run has a defined end (see for example: https://dl.acm.org/doi/10.5555/800108.803514).