Recall that a stabilizer is an Abelian subgroup of the \(n\)-qubit Pauli group $$ \mathcal{G}_n = \{\pm 1,\pm i\}\cdot\{I,X,Y,Z\}^{\otimes n}. $$ We write \(\mathcal{S}=\langle g_1,\dots, g_{n-k}\rangle\) to mean the stabilizer \(\mathcal{S}\) is generated by the elements in angled brackets. If these elements are independent, the resulting subspace is \(2^k\)-dimensional, and hence can encode \(k\) qubits. The code corresponding to \(\mathcal{S}\) is the subspace \(C\subseteq (\mathbb{C}^2)^{\otimes n}\) defined through $$ C = \{\ket{\psi}\in(\mathbb{C}^2)^{\otimes n}: S\ket{\psi}=\ket{\psi}\ \forall S\in\mathcal{S}\}. $$ Clearly, we cannot have \(-I\in \mathcal{S}\). This rules out phases of \(\pm i\) appearing in the stabilizers.
If \(P\) is a Pauli operator that commutes with every stabilizer generator, then either (i) \(P\) is proportional to an element of the stabilizer, or (ii) \(P\) is in the normalizer, but is not proportional to a stabilizer itself. If \(d\) is the distance of our code and \(P=P_aP_b\), where \(P_a,P_b\) are Pauli operators of weight1 at most \((d-1)/2\), Option (ii) cannot apply. Why? Recall $$ d=\min_{P\in N(S)\backslash S}\textnormal{wt}(P) $$ where we are taking Paulis up to a phase, and so we know \(P_aP_b\) is not in \(N(S)\backslash S\). But since it commutes with all the generators, it is in \(\langle S, \pm i\rangle\). Or in words, it is proportional to a stabilizer. This is Option (i).
Hence we know \(P_aP_b = \alpha S\) for some \(S\). Also, we can always write \(P_a=\eta_a R S_a\) for some Pauli operator \(R\) with \(+1\) coefficient, stabilizer \(S_a\), and phase \(\eta_a\). This is because \(P_a\), viewed as an element of the quotient group (modding out phases), lies in some coset of the stabilizer subgroup with representative element \(R\) (which can be chosen to be Hermitian when viewed as an element of the Pauli group). Writing down the phase explicitly, we need the \(\eta_a\) factor. Crucially, we also know \(P_b = \eta_b R S_b\) with the same operator \(R\) because they lie in the same coset of the quotient group.
Now let's see what happens to an element of the codespace when the noise operation \(\mathcal{N}(\rho)=\sum_a E_a \rho E_a^\dagger\) takes place. We can write \(E_a=\sum_{b}c_{ab}P_b\) where \(P_b\) are Hermitian Pauli operators and therefore $$ \begin{align} \mathcal{N}(\rho)&=\sum_{a}\sum_{bb^\prime}c_{ab}c_{ab^\prime}^*P_b\rho P_{b^\prime}\\ &=\sum_{ab}\chi_{ab}P_a\rho P_{b} \end{align} $$ where \(\chi_{ab}\equiv \sum_{d}c_{da}c_{db}^*\) is a Hermitian, positive semidefinite matrix. If \(\Pi_x\) denotes a syndrome measurement for the stabilizer code, we have the (unnormalized) post-measurement state $$ \begin{align} \Pi_x\mathcal{N}(\rho)\Pi_x &= \sum_{ab}\chi_{ab}P_a\Pi_{x+s(P_a)}\rho\Pi_{x+s(P_b)}P_b \end{align} $$ where \(s(P)\in\{0,1\}^{n-k}\) denotes the syndrome of \(P\) and addition of bitstrings is mod 2. Now, if \(\rho=|\psi\rangle\langle\psi|\) where \(\ket{\psi}\) is in the codespace, the only way a term in the sum above can be nonzero is if both projections are onto the codespace, i.e., \(x+s(P_a)=x+s(P_b)=0\); otherwise, the projection is orthogonal to the codespace. But this means \(s(P_a)+s(P_b)=0\), or equivalently \(s(P_aP_b)=0\). Thus \(P_aP_b\) is an operator that commutes with every stabilizer generator and, so long as we can take for granted that it is not a logical operator, it is proportional to a stabilizer.
We can now use the above reasoning to arrive at the post-measurement state being equal to $$ \begin{align} \sum_{\substack{a,b\\P_a, P_b\ \textnormal{in the same coset}}} \chi_{ab}P_a\rho P_b&= \sum_R\sum_{(a,b)\leftarrow R} \chi_{ab}P_a\rho P_b\\ &= \sum_R\sum_{(a,b)\leftarrow R} \chi_{ab}\eta_a\eta_b^* RS_a\rho S_bR\\ &=\sum_R \alpha_R R\rho R \end{align} $$ where the coset on the left-hand side is referring to cosets of the quotient group, \((a,b)\leftarrow R\) is notation I made up for referring to \(P_a\) and \(P_b\) being in the same coset, and \(\alpha_R\equiv\sum_{(a,b)\leftarrow R} \chi_{ab}\eta_a\eta_b^*\) is nonnegative because \(\chi_{ab}\) is psd. Also, we have used $$ \begin{align} \bra{\psi}P_b&= (P_b \ket{\psi})^{\dagger}\\ &= (\eta_b R S_b\ket{\psi})^\dagger\\ &= \eta_b^*\bra{\psi}S_b R \end{align} $$ which in turn relies on the fact that we are working with Hermitian Paulis. (This also means \(\eta_b\in \{\pm 1\}\), but we don't really need to use this.)
In the end, we can conclude that \(\Pi_x\mathcal{N}(\rho)\Pi_x\) is proportional to a convex combination of random Pauli errors on the original state, and we have our result.