Event-by-event correlations of azimuthal anisotropy Fourier coefficients ($v_n$) in 8.16 TeV pPb collision data, collected by the CMS experiment at the LHC, are extracted using a subevent four-particle cumulant technique. The pseudorapidity range of the CMS tracker from $-$2.4 to 2.4 is divided into either two, three, or four distinct subevent regions. Each combination of four particles constructed from tracks with transverse momentum between 0.3 and 3.0 GeV is then analyzed in terms of how the particles populate the subevents. Using the subevent technique, correlations between $v_n$ of different orders are measured as functions of particle multiplicity and compared to the standard cumulant method without subevents. At high multiplicity, the $v_2$ and $v_3$ coefficients exhibit an anticorrelation; this behavior is observed consistently using various methods. The $v_2$ and $v_4$ correlation strength is found to depend on the number of subevents used in the calculation. At low-multiplicities, the results from different methods diverge because of different contributions of few-particle correlations. These findings significantly lower the multiplicity range that was previously established for the onset of collective behavior in small systems.