Statistical analysis has become increasingly important in optically stimulated luminescence (OSL) dating since it has become possible to measure signals at the single-grain scale. The accuracy of large chronological datasets can benefit from the inclusion, in chronological modelling, of stratigraphic constraints and shared systematic errors. Recently, a number of Bayesian models have been developed for OSL age calculation; the R package "BayLum"presented herein allows different models of this type to be implemented, particularly for samples in stratigraphic order which share systematic errors. We first show how to introduce stratigraphic constraints in BayLum; then, we focus on the construction, based on measurement uncertainties, of dose covariance matrices to account for systematic errors specific to OSL dating. The nature (systematic versus random) of errors affecting OSL ages is discussed, based - as an example - on the dose rate determination procedure at the IRAMAT-CRP2A laboratory (Bordeaux). The effects of the stratigraphic constraints and dose covariance matrices are illustrated on example datasets. In particular, the benefit of combining the modelling of systematic errors with independent ages, unaffected by these errors, is demonstrated. Finally, we discuss other common ways of estimating dose rates and how they may be taken into account in the covariance matrix by other potential users and laboratories. Test datasets are provided as a Supplement to the reader, together with an R markdown tutorial allowing the reproduction of all calculations and figures presented in this study.