For particle physics observables at colliders such as the LHC at CERN, it has been common practice for many decades to estimate the theoretical uncertainty by studying the variations of the predicted cross sections with a priori unpredictable scales. In astroparticle physics, this has so far not been possible, since most of the observables were calculated at Born level only, so that the renormalization scheme and scale dependence could not be studied in a meaningful way. In this paper, we present the first quantitative study of the theoretical uncertainty of the neutralino dark matter relic density from scheme and scale variations. We first explain in detail how the renormalization scale enters the tree-level calculations through coupling constants, masses and mixing angles. We then demonstrate a reduction of the renormalization scale dependence through one-loop SUSY-QCD corrections in many different dark matter annihilation channels and enhanced perturbative stability of a mixed on-shell/$barrm DR$ renormalization scheme over a pure $barrm DR$ scheme in the top-quark sector. In the stop-stop annihilation channel, the Sommerfeld enhancement and its scale dependence are shown to be of particular importance. Finally, the impact of our higher-order SUSY-QCD corrections and their scale uncertainties are studied in three typical scenarios of the phenomenological Minimal Supersymmetric Standard Model with eleven parameters (pMSSM-11). We find that the theoretical uncertainty is reduced in many cases and can become comparable to the size of the experimental one in some scenarios.