Reviewer 1.6

[evogytis]

P 15+, I identify the following weaknesses in the simulation study, and I would candidly discuss these in the Discussion section:

• Dispersion parameter is set at 0.1 based on a previous branching process study, however there's not great support for re-using the value for a different epidemic. Was an effort made to estimate this parameter?
• The sampling model accounts for 'sequencing capacity', which is good. But is it the case that sequencing probability varied over the course of the epidemic?
• The 'model fitting' methods are ad-hoc. This is a situation where approximate Bayes computation would be appropriate. But I would try GEE/method of moments. I suspect the first three standardized moments of cluster sizes would be more useful than using mean/median/sd. - The moments could also be computed analytically if you dispense with the 'sequencing capacity' aspect of the sampling model, so would be very fast.
• Can you show a comparison of the empirical distribution of cluster sizes and the fitted distribution? A QQ plot? Are there any aspects of the empirical distribution that are not well characterised by the branching process model, eg at the tails?

[evogytis]

The reviewer wants candid discussion points, but I think we can re-run a few things.

• Our options are: arguing that dispersion is a difficult parameter to estimate with our data or doing an additional analysis with a different dispersion parameter. The latter isn't complicated, can be run while other stuff runs in the background and will find room in supp.
• Very good point. Not sure what the best way of addressing this is.
• The good news is that all simulations were saved as cluster sizes, rather than summary statistics. Can easily remake the figure with mean, stdev and skewness, though need to play around with which of those three factors make for a good pair of coordinates on the plot.
• I wanted to do a more direct comparison of empirical and simulated distributions, but couldn't think of a good vis for it. QQ is worth a shot.

[evogytis]

• Should be fixed now.
• Still not sure what do about this. Raise as point in discussion?
• Switched to using number of clusters, mean cluster size, cluster size standard deviation, cluster size skewness. Not sure what to do about GEE/method of moments. Anyone know how to do this?
• QQ plots are in the revisionsEpi Jupyter notebook. Shall we point to the notebook or include another supp figure?

[maxbiostat]

Right:

• I think what we've done is good for now. If there are any more qualms we have one more card up our sleeve: use the MLE of the dispersion from the sequence data as (I think) an upper bound.

• This is a great point. But as with most great comments, this opens a whole can of worms that is just outside the scope of the paper (imho).

• I know how to do both generalised estimating equations (GEE) and moments. What the reviewer doesn't know is that the mean is the sufficient statistic for R0. For k the story is a bit different, as the second moment is not exactly sufficient/ancillary but can be used to obtain an estimator. So it turns out that, under this model, mean and sd are the statistics to compute. We can throw some skewness in, but I'm willing to bet @evogytis 's scientific reputation that it won't improve "inference". The reviewer is correct about the methods being ad-hoc. What said reviewer also doesn't know is that ABC for instance is indeed quite suitable to the problem, but also requires a huge amount of tuning and checking, which I think, again, is outside the scope.

• I have no qualms either way. Letting @rambaut and @trvrb do their advisory thing here.

[trvrb]

I'm agreement about approaches here. Threading individual points:

1. Examining dispersion the same way we examined bias seems like absolutely the way to address this.
2. I think adding the caveat sentence (as you've done in e3078487cf3c61e39150d7e33c29c2a1d283f9fe) is the best we can hope to do.
3. I disagree a bit here in that I believe the Monte Carlo approach is a principled one. If you have some pieces of data you can simulate under a prior and see which simulations recover your data. Monte Carlo is equivalent in this case to an MCMC approach, just less efficient in that many simulations will be far off base. The heuristic bit was the choice of bounds for which simulations to include (the exact approach would be to keep only simulations that recover the exact distribution of sequence cluster sizes). The pursued approach of 95% HPDs on summary statistics seems reasonable, but is not exact. I would modify the text to better justify the use of Monte Carlo in this scenario.
4. The QQ plots seem good. I would include in one of the supplemental (extended data) figures. Probably just show a single QQ for bias=2, k=0.1 (our best fitting model). This is a "posterior predictive check". I would be more interpretable to me to have the simple distribution plotted as well (cluster size histogram) with empirical distribution along side mean and interval for simulations.