  One of the primary challenges in enabling the scientific potential of 21 cm
intensity mapping at the Epoch of Reionization (EoR) is the separation of
astrophysical foreground contamination. Recent works have claimed that Gaussian
process regression (GPR) can robustly perform this separation, particularly at
low Fourier $k$ wavenumbers where the signal reaches its peak signal-to-noise
ratio. We revisit this topic by casting GPR foreground subtraction (GPR-FS)
into the quadratic estimator formalism, thereby putting its statistical
properties on stronger theoretical footing. We find that GPR-FS can distort the
window functions at these low k modes, which, without proper decorrelation,
make it difficult to probe the EoR power spectrum. Incidentally, we also show
that GPR-FS is in fact closely related to the widely studied optimal quadratic
estimator. As a case study, we look at recent power spectrum upper limits from
the Low Frequency Array (LOFAR) that utilized GPR-FS. We pay close attention to
their normalization scheme, showing that it is particularly sensitive to signal
loss when the EoR covariance is misestimated. This implies possible
ramifications for recent astrophysical interpretations of the LOFAR limits,
because many of the EoR models ruled out do not fall within the bounds of the
covariance models explored by LOFAR. Being more robust to this bias (although
not entirely free of it), we conclude that the quadratic estimator is a more
natural framework for implementing GPR-FS and computing the 21 cm power
spectrum.
