I'm trying to prove that the following equality holds, where $X$ is an arbitrary vector field on a Riemannian manifold (since I already know the left side is also equal to $\nabla_{X} R$, if I proved it the second Bianchi identity would then follow):
$$-2 \Delta \nabla_{i} X^{i}-2 R_{i j} \nabla_{i} X^{j}+\nabla_{i} \nabla_{j} \nabla_{i} X^{j}+\nabla_{i} \nabla_{j} \nabla_{j} X^{i} =2 X^{i} \nabla_{j} R_{i j}$$
Now, by my computations:
$$\begin{aligned} \nabla_{i} \nabla_{j} \nabla_{j} X^{i} &= \nabla_{i} \nabla_{j} \nabla_{j} X^{i} - \nabla_{j} \nabla_{i} \nabla_{j} X^{i} + \nabla_{j} \nabla_{i} \nabla_{j} X^{i} \\ &= R_{ijk}^{i} \nabla_{j} X^{k} + \nabla_{j} \nabla_{i} \nabla_{j} X^{i} \\ &= R_{ij} \nabla_{i} X^{j} + \nabla_{j} \nabla_{i} \nabla_{j} X^{i}\end{aligned}$$
$$\begin{aligned} \nabla_{j} \nabla_{i} \nabla_{j} X^{i} &= \nabla_{j} \nabla_{i} \nabla_{j} X^{i} - \nabla_{j} \nabla_{j} \nabla_{i} X^{i} + \nabla_{j} \nabla_{j} \nabla_{i} X^{i} \\ &= \nabla_{j}(R_{ijk}^{i} X^{k}) + \Delta \nabla_{i} X^{i} \\ &= \nabla_{j} (R_{jk} X^k) + \Delta \nabla_{i} X^{i} \\ &= \nabla_{i}(R_{ij} X^{j}) + \Delta \nabla_{i} X^{i} \end{aligned}$$
So we have:
$$\begin{aligned} &-2 \Delta \nabla_{i} X^{i}-2 R_{i j} \nabla_{i} X^{j}+\nabla_{i} \nabla_{j} \nabla_{i} X^{j}+\nabla_{i} \nabla_{j} \nabla_{j} X^{i} \\ &= -2 \Delta \nabla_{i} X^{i}-2 R_{i j} \nabla_{i} X^{j} + \nabla_{j} \nabla_{i} \nabla_{j} X^{i} +\nabla_{i} \nabla_{j} \nabla_{j} X^{i} \\ &= -2 \Delta \nabla_{i} X^{i}-2 R_{i j} \nabla_{i} X^{j} + \nabla_{i}(R_{ij} X^{j}) + \Delta \nabla_{i} X^{i} + R_{ij} \nabla_{i} X^{j} + \nabla_{j} \nabla_{i} \nabla_{j} X^{i} \\ &= -2 \Delta \nabla_{i} X^{i}-2 R_{i j} \nabla_{i} X^{j} + 2\nabla_{i}(R_{ij} X^{j}) + 2\Delta \nabla_{i} X^{i} + R_{ij} \nabla_{i} X^{j} \\ &= 2\nabla_{i}(R_{ij} X^{j}) - R_{i j} \nabla_{i} X^{j} \\ &= 2X^{j} \nabla_{i}R_{ij} + R_{ij} \nabla_{i} X^{j}\end{aligned}$$
What's going on here? Where did I make a mistake so that an extra $R_{ij} \nabla_{i} X^{j}$ could come in? I'd appreciate some help.
Your mistake is in the second equality of your first display. In short, the curvature operator is a derivation, so it needs to be applied in both components. Thus, if $T^{ji}$ is a $(2,0)$-tensor field, then \begin{align*} \nabla_i\nabla_j T^{ji} - \nabla_j\nabla_i T^{ji} & = R_{ij}{}^j{}_k T^{ki} + R_{ij}{}^i{}_k T^{jk} \\ & = -R_{ik} T^{ki} + R_{jk} T^{jk} \\ & = 0 . \end{align*} (Note: I am using abstract index notation in this computation.) Applying this with $T^{ji}=\nabla^j X^i$ yields your expected result; in your computation, you missed the first summand, which led to your extra factor of $R_{ij}\nabla^iX^j$.