We extend the self-energy functional theory to the case of interacting lattice bosons in the presence of symmetry breaking and quenched disorder. The self-energy functional we derive depends only on the self-energies of the disorder-averaged propagators, allowing for the construction of general non-perturbative approximations. Using a simple single-site reference system with only three variational parameters, we are able to reproduce numerically exact quantum Monte Carlo (QMC) results on local observables of the Bose–Hubbard model with box disorder with high accuracy. At strong interactions, the phase boundaries are reproduced qualitatively but shifted with respect to the ones observed with QMC due to the extremely low condensate fraction in the superfluid phase. Deep in the strongly-disordered weakly-interacting regime, the simple reference system employed is insufficient and no stationary solutions can be found within its restricted variational subspace. By systematically analyzing thermodynamical observables and the spectral function, we find that the strongly interacting Bose glass is characterized by different regimes, depending on which local occupations are activated as a function of the disorder strength. We find that the particles delocalize into isolated superfluid lakes over a strongly localized background around maximally-occupied sites whenever these sites are particularly rare. Our results indicate that the transition from the Bose glass to the superfluid phase around unit filling at strong interactions is driven by the percolation of superfluid lakes which form around doubly occupied sites.