The immersed boundary method is a mathematical framework for modeling fluid-structure interaction. This formulation describes the momentum, viscosity, and incompressibility of the fluid-structure system in Eulerian form, and it uses Lagrangian coordinates to describe the structural deformations, stresses, and resultant forces. Integral transforms with Dirac delta function kernels connect the Eulerian and Lagrangian frames. The fluid and the structure are both typically treated as incompressible materials. Upon discretization, however, the incompressibility of the structure is only maintained approximately. To obtain an immersed method for incompressible hyperelastic structures that is robust under large structural deformations, we introduce a volumetric energy in the solid region that stabilizes the formulation and improves the accuracy of the numerical scheme. This formulation augments the discrete Lagrange multiplier for the incompressibility constraint, thereby improving the original method's accuracy. This volumetric energy is incorporated by decomposing the strain energy into isochoric and dilatational components, as in standard solid mechanics formulations of nearly incompressible elasticity. We study the performance of the stabilized method using several quasi-static solid mechanics benchmarks, a dynamic fluid-structure interaction benchmark, and a detailed three-dimensional model of esophageal transport. The accuracy achieved by the stabilized immersed formulation is comparable to that of a stabilized finite element method for incompressible elasticity using similar numbers of structural degrees of freedom.