Nuclei located in the neutron-deficient Pb region have a complex structure, rapidly evolving as a function of neutron and proton numbers. The most famous example is $^{186}$Pb where the three lowest levels are $0^+$ states, the two excited $0^+$ ones being located at low excitation energy around 600 keV. Coexisting structures with different properties are found in the neighboring nuclei. Many experiments have been performed over the last few years in which in-band and out-of-band $\gamma$ transition probabilities have been measured. A detailed interpretation of experimental data requires the use of a method going beyond a mean-field approach that permits to determine spectra and transition probabilities. Such methods have already been applied to selected isotopes in this mass region. Our aim is to provide a systematic investigation of this mass region in order to determine how well experimental data can be understood using a state-of-the-art method for nuclear structure. The starting point of our method is a set of mean-field wave functions generated with a constraint on the axial quadrupole moment and using a Skyrme energy density functional. Correlations beyond the mean field are introduced by projecting mean-field wave functions on angular-momentum and particle number and by mixing the symmetry restored wave functions as a function of the axial quadrupole moment. A detailed comparison with the available data is performed for energies, charge radii, spectroscopic quadrupole moments, $E0$ and $E2$ transition probabilities for the isotopic chains of neutron deficient Hg, Pb, Po and Rn. The connection between our results and the underlying mean field is also analyzed. Qualitative agreement with the data is obtained although our results indicate that the actual energy density functionals have to be improved further to achieve a quantitative agreement.