Causal discovery procedures aim to deduce causal relationships among variables in a multivariate dataset. While various methods have been proposed for estimating a single causal model or a single equivalence class of models, less attention has been given to quantifying uncertainty in causal discovery in terms of confidence statements. A primary challenge in causal discovery of directed acyclic graphs is determining a causal ordering among the variables, and our work offers a framework for constructing confidence sets of causal orderings that the data do not rule out. Our methodology specifically applies to identifiable structural equation models with additive errors and is based on a residual bootstrap procedure to test the goodness-of-fit of causal orderings. We demonstrate the asymptotic validity of the confidence set constructed using this goodness-of-fit test and explain how the confidence set may be used to form sub/supersets of ancestral relationships as well as confidence intervals for causal effects that incorporate model uncertainty.