The detection of a stochastic background of low-frequency gravitational waves by pulsar-timing and astrometric surveys will enable tests of gravitational theories beyond general relativity. These theories generally permit gravitational waves with non-Einsteinian polarization modes, which may propagate slower than the speed of light. We use the total-angular-momentum wave formalism to derive the angular correlation patterns of observables relevant for pulsar timing arrays and astrometry that arise from a background of subluminal gravitational waves with scalar, vector, or tensor polarizations. We find that the pulsar timing observables for the scalar longitudinal mode, which diverge with source distance in the luminal limit, are finite in the subluminal case. Furthermore, we apply our results to $f(R)$ gravity, which contains a massive scalar degree of freedom in addition to the standard transverse-traceless modes. The scalar mode in this $f(R)$ theory is a linear combination of the scalar-longitudinal and scalar-transverse modes, exciting only the monopole and dipole for pulsar timing arrays and only the dipole for astrometric surveys.