Answer from the literature
Background
The question of the validity of the incoherent approximation is explored in detail in this publication:
Bravyi, S., Englbrecht, M., König, R. et al. Correcting coherent errors with surface codes. npj Quantum Inf 4, 55 (2018). https://doi.org/10.1038/s41534-018-0106-y
They used a mapping of the surface code to Majorana fermions in order to do an efficient simulation that includes coherent noise, with code distance up to 49, which is actually really cool.
The noise model they apply is $U=(e^{iZ\theta})^{\otimes n}$, namely, an equal $Z$ rotation with angle $\theta$ on all qubits. I think this model is a bit artificial since I don't believe a control error will be so uniformly biased. In fact, I'm not even sure why a uniform Z rotation on all the qubits causes any error, and why it doesn't just amount to changing the starting time of the experiment. But let's go with it...
Results
First of all, they show the threshold angle associated with the logical error threshold is about $0.1\pi$ which comes out consistent with the well known 1% logical error rate, using $p=\sin^2\theta$.
However, here's the ratio of the Pauli twirled (incoherent) noise with $p=\sin^2\theta$, and you can see it's pretty bad below the threshold.
In any case you can see the result is quite pessimistic in this specific case.
My sketch of the general principles
As you model it - no, it isn't, at least formally speaking. Your control error model is coherent whereas your statistical Pauli error implementation is not. However, at the limit of small errors, they are effectively close enough.
You can see the difference if you look at the Kraus matrices. For example, if your target gate is $V$ and your (near identity) control error is $U=1+\sqrt{\epsilon} X$, then the Kraus matrix, which operates on your state as $\rho\to E\rho E^\dagger$, is $E = UV $, and so the transformation for the imperfect gate is:
$$
\rho\to\rho_c = (1+\sqrt{\epsilon} X)V\rho V^\dagger(1-\sqrt{\epsilon} X)\\=V\rho V^\dagger+\sqrt\epsilon(XV\rho V^\dagger-V\rho V^\dagger X) - \epsilon(XV\rho V^\dagger X)
$$
However, in the statistical bitflip error, the transformation is:
$$
\rho\to \rho_{ic} = V\rho V^\dagger+\epsilon XV\rho V^\dagger X
$$
Basically, coherent and incoherent errors do not behave in the same way - coherent errors add up in amplitude, and incoherent errors add in up in squares of amplitudes.
However, as you wrote in your question, the ancilla measurement will project the state onto the even/odd parity subspaces. I.e., if we denote $\rho_c=UV\rho V^\dagger U^\dagger$, and $P_{\pm}$ to be our even/odd parity eigenstates, Then after the ancilla meaurement we will have, in the coherent error model, the state:
$$
pP_+\rho_c P_+ + (1-p)P_- \rho_c P_-.
$$
However, if this is the only error in the circuit, then we can assume $\rho_1 = V\rho V^\dagger$ is in a definite parity state, and therefore the cross terms such as $X\rho_1$ and $\rho_1 X$ cancel out, because if e.g $\rho_1$ is in an even parity state, then:
$$
P_+ X\rho_1 P_+ = X P_- \rho_1 P_+ = X P_- P_+ \rho P_+ ^2 = 0
$$
and the same for the other cross terms. Thus, after the ancilla measurement:
$$
\rho_{c} \to \rho_{ic}
$$.
That being said, this assumes that the only error which occured is the single coherent error. If however state will not have definite parity before this error, e.g if you have two coherent errors that affect the same qubit, this will no longer be correct, and the errors in general will be larger.
