We present a multi-site formulation of mean-field theory applied to the disordered Bose-Hubbard model. In this approach the lattice is partitioned into clusters, each isolated cluster being treated exactly, with inter-cluster hopping being treated approximately. The theory allows for the possibility of a different superfluid order parameter at every site in the lattice, such as what has been used in previously published site-decoupled mean-field theories, but a multi-site formulation also allows for the inclusion of spatial correlations allowing us, e.g., to calculate the correlation length (over the length scale of each cluster). We present our numerical results for a two-dimensional system. This theory is shown to produce a phase diagram in which the stability of the Mott insulator phase is larger than that predicted by site-decoupled single-site mean-field theory. Two different methods are given for the identification of the Bose glass-to-superfluid transition, one an approximation based on the behaviour of the condensate fraction, and one of which relies on obtaining the spatial variation of the order parameter correlation. The relation of our results to a recent proposal that both transitions are non self-averaging is discussed.