The estimation of the rock mass deformation modulus ( \({\mathrm{D}}_{\mathrm{f}}\) ) via an empirical model has an approximately half-century history. However, reliable estimation of \({\mathrm{D}}_{\mathrm{f}}\) has been a challenging task in terms of the theoretical concept of input parameters and data analysis methods. Analytical models present the principal input parameters; however, the concept of principal input parameters (PIP) will develop with an emphasis on in situ stress by participating in the confined Young's modulus and shear and normal joint stiffness at a specified normal stress. A review of seventy empirical models revealed that the majority of existing empirical relationships suffer from a lack of PIP. In this study, based on the geological strength index (29 < GSI < 83), confined Young's modulus, and shear and normal joint stiffness at specified normal stresses, the deformation modulus (1 < \({\mathrm{D}}_{\mathrm{f}}\) < 39.8 GPa) is forecasted by an empirical model. The database copes quite well with eighty-two data sets of different rock types. A new multigene genetic program (MGP) integrates five genes with a maximum depth of three as an optimal mathematical relationship in terms of fitness functions. A comparison of the estimated deformation modulus with several existing empirical models based on the same database confirms the superiority of the new MGP. The integration of the analytical base PIP improves the global acceptability of empirical models in analytical or numerical analysis.