Shape coexistence has been a subject of great interest in nuclear physics for many decades. In the context of the nuclear shell model, intruder excitations may give rise to remarkably low-lying excited 0+ states associated with different intrinsic shapes. In heavy open-shell nuclei, the dimension of the shellmodel configuration space that includes such intruder excitations becomes exceedingly large, thus requiring a drastic truncation scheme. Such a framework has been provided by the interacting boson model (IBM). In this article we address the phenomenon of shape coexistence and its relevant spectroscopy from the point of view of the IBM. A special focus is placed on the method developed recently which makes use of the link between the IBM and the self-consistent mean-field approach based on the nuclear energy density functional. The method is extended to deal with various intruder configurations associated with different equilibrium shapes. We assess the predictive power of the method and suggest possible improvements and extensions, by considering illustrative examples in the neutron-deficient Pb region, where shape coexistence has been experimentally studied.