Long-range correlations for pairs of charged particles with two-particle angular correlations are studied in $pp$ at $\sqrt{\text{s}}=13$ TeV with various Monte Carlo generators. The correlation functions are constructed as functions of the relative azimuthal angle $\Delta \phi$ and the pseudorapidity separation $\Delta \eta$ for pairs of different particle species with the identified hadrons such as $\pi$, $K$, $p$, and $\Lambda$ in wide $\Delta \eta$ ranges. Fourier coefficients are extracted for the long-range correlations in several-multiplicity classes using a low-multiplicity template fit method. The method allows one to subtract the enhanced away-side jet fragments in high-multiplicity events with respect to low-multiplicity events. However, we found that due to a kinematic bias on jets and differing model implementation of flow and jet components, subtracting the nonflow contamination in small systems can bias the results. It is found that the pythia8 default model where the presence of the collective flow is not expected but the bias results in very large flow. Also extracting flow signal from the epos4 and pythia8 string shoving models is not possible because of the flow signal introduced in the low-multiplicity events. Only a multiphase transport string melting model among studied model calculations is free from this bias and shows a mass ordering at low $p_T$ and particle type grouping in the intermediate-$p_T$ range. This feature was first observed in large systems, but the mass ordering in small systems differs from that observed in large collision systems.