Your confusion seems to be notation related - in conventional abstract index notation, $\nabla_i (\nabla_j f)$ and $\nabla_i \nabla_j f$ are different things! The first means what you think it does - since $f$ is a function, $\nabla_j f = \partial_j f$ is also a function, and thus $\nabla_i(\nabla_j f) = \partial_i \partial_j f$. However, when the parentheses are omitted, all the covariant derivatives are taken before any indices are applied - that is, $$\nabla_i \nabla_j f := \nabla^2 f (\partial_i, \partial_j) = \nabla_i (\nabla_j f) - \nabla_{\nabla_i \partial_j} f.$$
This can be very confusing but it makes computations easier to write down.
Thus the resolution of your problem is that the expression
$$\nabla_j\nabla_k\nabla_if - \nabla_k\nabla_j\nabla_if$$is in fact commuting second covariant derivatives of a one-form $\nabla f$ - it is only at the very end that we take the $i$th component.
If you know how to commute covariant derivatives of an arbitrary tensor, then this is not so hard: since the metric commutes with covariant derivatives, we have $$\begin{align}
\nabla_k \Delta T - \Delta \nabla_k T &= g^{ij}(\nabla_k \nabla_i \nabla_j T - \nabla_i \nabla_j \nabla_k T) \\
&= g^{ij}([\nabla_k, \nabla_i] \nabla_j T +\nabla_i\nabla_k\nabla_j T- \nabla_i \nabla_j \nabla_k T) \\
&= g^{ij}([\nabla_k, \nabla_i] \nabla_j T + \nabla_i([\nabla_k,\nabla_j] T))\tag{1}.\end{align}$$
Here I'm using the abstract index notation, so e.g. $\nabla_i \nabla_j \nabla_k T = (\nabla^3 T)(\partial_i, \partial_j, \partial_k, \ldots)$ is the covariant third derivative.
Since commutators of covariant derivatives always end up just being some contractions with the curvature, we can immediately see from $(1)$ that the formula will have the form $$[\nabla,\Delta]T = R* \nabla T + \nabla( R * T) = R * \nabla T + \nabla R * T \tag{2},$$ where the notation $S * T$ denotes an arbitrary linear combination of contractions of $S \otimes T$. Thus you should expect not only curvature terms, but also derivatives of the curvature. (That is, unless $n = 0$ - meaning $T$ is just a scalar - in which case the second term vanishes because the Hessian of a scalar is symmetric.) For many purposes
$(2)$ is all you need.
To get an explicit expansion of $(1)$ in terms of the curvature tensor, we first need to write down an expression for $[\nabla_k, \nabla_j] T$.$\def\t{\otimes}$ When we take the second covariant derivative of $T = T_{i_1 \ldots i_n}dx^{i_1}\t\cdots\t dx^{i_n}$ and fully expand using the product rule, all the cross terms will vanish thanks to the antisymmetrization; so the derivative commutator also satisfies a Leibniz rule: we find $$[\nabla_k, \nabla_j] T = [\nabla_k,\nabla_j](T_{i_1 \ldots i_n})dx^{i_1}\t\cdots\t dx^{i_n}\\ + T_{i_1 \ldots i_n} ([\nabla_k,\nabla_j]dx^{i_1})\t\cdots \t dx^{i_n} \\ \vdots \\+ T_{i_1 \ldots i_n} dx^{i_1}\t\cdots \t([\nabla_k,\nabla_j]dx^{i_n}).$$
Since the curvature tensor is characterized by $R_{kj}{}^i{}_l dx^l = [\nabla_k,\nabla_j]dx^i$ and the first two derivatives of a scalar commute, this becomes $$[\nabla_k, \nabla_j] T=T_{i_1 \ldots i_n}(R_{kj}{}^{i_1}{}_ldx^l \t \cdots \t dx^{i_n} + \cdots + R_{kj}{}^{i_n}{}_l dx^{i_1} \t \cdots \t dx^l),$$ which (doing some index gymnastics) we can write as
$$[\nabla_k,\nabla_j]T_{i_1\ldots i_n} = R_{kj}{}^l{}_{i_1} T_{l i_2 \ldots i_n}+\cdots + R_{kj}{}^l{}_{i_n}T_{i_1 \ldots i_{n-1} l}.$$
Since you can now commute derivatives of an arbitrary $n$-tensor, you can do it for an $n+1$-tensor too; so if you really want to write out the full mess you can plug this in to $(1)$ for the two commutator terms and then expand everything with the product rule. I don't think there's much insight to be gained by typing it all out - it's just going to be a mess of contractions, some of them of $R$ with $\nabla T$ and some of $\nabla R$ with $T$. Perhaps try it yourself for a very small value of $n$.
Best Answer
You're almost there. Using the Einstein summation convention, so that we sum over repeated indices, we have
$$\Delta \nabla_i f = \nabla^j \nabla_j \nabla_i f = \nabla^j \nabla_i \nabla_j f = - R^j{}_i{}^k{}_j \nabla_k f +\nabla_i \nabla^j \nabla_j f = g^{jk} \operatorname{Ric}_{ij} \nabla_k f + \nabla_i \Delta f.$$
The third equality follows from applying the Ricci identity for $1$-forms that you mention in your edit to $\nabla_j f$, and the fourth from the symmetries of the curvature tensor $R$ and the definitions of $\operatorname{Ric}$ and $\Delta$.