According to Dyson's three fold way, from the viewpoint of global time reversal symmetry there are three circular ensembles of unitary random matrices relevant to the study of chaotic spectra in quantum mechanics. These are the circular orthogonal, unitary and symplectic ensembles, denoted COE, CUE and CSE respectively. For each of these three ensembles and their thinned versions, whereby each eigenvalue is deleted independently with probability $1-ξ$, we take up the problem of calculating the first two terms in the scaled large $N$ expansion of the spacing distributions. It is well known that the leading term admits a characterisation in terms of both Fredholm determinants and Painlevé transcendents. We show that modifications of these characterisations also remain valid for the next to leading term, and that they provide schemes for high precision numerical computations. In the case of the CUE there is an application to the analysis of Odlyzko's data set for the Riemann zeros, and in that case some further statistics are similarly analysed.