  For ${\cal N}=1$ SQED with $N_f$ flavors regularized by higher derivatives in
the general $\xi$-gauge we calculate the three-loop anomalous dimension of the
matter superfields defined in terms of the bare coupling constant and
demonstrate its gauge independence. After this the four-loop $\beta$-function
defined in terms of the bare coupling constant is obtained with the help of the
NSVZ equation, which is valid for these renormalization group functions in all
loops. Next, we calculate the three-loop anomalous dimension and the four-loop
$\beta$-function defined in terms of the renormalized coupling constant for an
arbitrary subtraction scheme supplementing the higher derivative
regularization. Then we construct a renormalization prescription for which the
results coincide with the ones in the $\overline{\mbox{DR}}$-scheme and
describe all NSVZ schemes in the considered approximation. Also we demonstrate
the existence of a subtraction scheme in which the anomalous dimension does not
depend on $N_f$, while the $\beta$-function contains only terms of the first
order in $N_f$. This scheme is obtained with the help of a finite
renormalization compatible with a structure of quantum corrections and is NSVZ.
The existence of such an NSVZ scheme is also proved in all loops.
