ALICE is a large experiment at the CERN Large Hadron Collider. Located 52 meters underground, its detectors are suitable to measure muons produced by cosmic-ray interactions in the atmosphere. In this paper, the studies of the cosmic muons registered by ALICE during Run 2 (2015--2018) are described. The analysis is limited to multimuon events defined as events with more than four detected muons ($N_\mu>4$) and in the zenith angle range $0^{\circ}<~\theta<~50^{\circ}$. The results are compared with Monte Carlo simulations using three of the main hadronic interaction models describing the air shower development in the atmosphere: QGSJET-II-04, EPOS-LHC, and SIBYLL 2.3. The interval of the primary cosmic-ray energy involved in the measured muon multiplicity distribution is about $ 4 \times 10^{15}<~E_\mathrm{prim}<~ 6 \times 10^{16}$~eV. In this interval none of the three models is able to describe precisely the trend of the composition of cosmic rays as the energy increases. However, QGSJET is found to be the only model capable of reproducing reasonably well the muon multiplicity distribution, assuming a heavy composition of the primary cosmic rays over the whole energy range, while SIBYLL and EPOS-LHC underpredict the number of muons in a large interval of multiplicity by more than $20\%$ and $30\%$, respectively. The rate of high muon multiplicity events ($N_\mu>100$) obtained with QGSJET and SIBYLL is compatible with the data, while EPOS-LHC produces a significantly lower rate ($55\%$ of the measured rate). For both QGSJET and SIBYLL, the rate is close to the data when the composition is assumed to be dominated by heavy elements, an outcome compatible with the average energy $E_\mathrm{prim} \sim 10^{17}$~eV of these events. This result places significant constraints on more exotic production mechanisms.