Shear-keyed inter-modular connections (IMCs) are integral components of high-rise modular steel structures (MSSs), providing robust interconnectivity to support grouped tubular columns across modules, thereby introducing column discontinuities and distinctive structural behavior. This study conducted a comprehensive numerical assessment and theoretical analysis of the axial compression behavior of grouped tubular columns based on a validated finite element model (FEM), which captured the member-to-structural level behavior of steel hollow section (SHS) columns and accommodated geometric imperfections. An FEM was initially developed and validated using 28 axial compression tests documented in the literature, comprising 15 tests on cold-formed and 13 on hot-rolled steel hollow section (SHS) columns. The primary parameters explored in tests included material properties (stainless/carbon), processing methods (cold-formed/hot-rolled), cross-section sizes (D/B), cross-sectional or member slenderness ratios (D/tc, B/tc, or Lc/r), and the number of columns (1, 7, and 11). A comprehensive parametric numerical study involving 103 grouped tubular column FEMs then investigated the influence of initial imperfection, shear-key height (Lt), thickness (tt), steel tube length (D), width (B), thickness (tc), and height (Lc) alongside the effects of space between tube and key, and the gap between tubes. The results indicated that the load-shortening behavior of the grouped columns consists of linear elastic, inelastic, and recession stages. The failure modes observed primarily displayed an S-shaped pair of inward and outward local buckling on the outer sides and double S-shaped local buckling on the interior sides. The buckling arose near the shear key or at 1/4 or 1/2 of the column height. None of the considered models experienced global buckling. Increasing tt, Lt, tc, D, or B enhances strength and stiffness, while Lc or Lc/r linearly affects stiffness and ductility. The columns’ nominal axial strength was reduced because of the shear keys, which decreased compression yielding and caused localized elastic buckling. Subsequently, the theoretical analysis revealed that the design codes do not capture this behavior, and thus, their capacity estimate yields inaccurate findings. This discrepancy renders existing code prediction equations, including those from Indian (IS800), New Zealand (NZS400), European (EC3:1-1), Canadian (CSA S16), American (AISC360-16), and Chinese (GB50017) standards, as well as the model proposed by Li et al., non-conservative. To assure conservative results, the paper recommended modification of existing standards and proposed prediction equations based on a fourth-order differential equation that describes the actual behavior of modular steel columns grouped with shear keys. The proposed design approach accurately predicted the axial compression capacity of modular steel-grouped columns, proving conservative yet effective. This provides valuable data that could transform design and construction techniques for MSSs, extending to various column and IMC forms through adaptable design parameters. This enhancement in structural performance and safety significantly contributes to the advancement of modular construction practices.