Short answer: the formula is true. Long answer: see below for the calculation.
Choose local coordinates $(x^1, \ldots, x^n)$ and a local frame field $(e_1, \ldots, e_r)$ for $E$, both of these over an open neighbourhood $U$, and let $\newcommand{\Tud}[3]{{#1}^{#2}_{\phantom{#2}{#3}}}$ $\nabla e_\mu = \Tud{\omega}{\nu}{\mu} \otimes e_\nu$ for some $\Tud{\omega}{\nu}{\mu} \in \Omega^1(U)$. (I am using the summation convention for repeated indices.) Let $\Tud{\Omega}{\nu}{\mu} \in \Omega^2(U)$ be such that $\nabla \nabla e_\mu = \Tud{\Omega}{\nu}{\mu} \otimes e_\nu$ Now, by the Leibniz rule, we have $\nabla \nabla e_\mu = \mathrm{d}\Tud{\omega}{\nu}{\mu} \otimes e_\nu - \Tud{\omega}{\nu}{\mu} \wedge \Tud{\omega}{\rho}{\nu} \otimes e_\rho$ thus, $\Tud{\Omega}{\nu}{\mu} = \mathrm{d}\Tud{\omega}{\nu}{\mu} + \Tud{\omega}{\nu}{\rho} \wedge \Tud{\omega}{\rho}{\mu}$
Since $\nabla (W^\mu e_\mu) = \mathrm{d} W^\mu \otimes e_\mu + W^\mu \Tud{\omega}{\nu}{\mu} \otimes e_\nu$ we have $\nabla \nabla (W^\mu e_\mu) = \mathrm{d} W^\mu \wedge \Tud{\omega}{\nu}{\mu} \otimes e_\nu + W^\mu \, \mathrm{d} \Tud{\omega}{\nu}{\mu} \otimes e_\nu - \mathrm{d}W^\mu \wedge \Tud{\omega}{\nu}{\mu} \otimes e_\nu - W^\mu \Tud{\omega}{\nu}{\mu} \wedge \Tud{\omega}{\rho}{\nu} \otimes e_\rho$ thus, $\nabla \nabla$ is indeed $C^\infty(U)$-linear, with $\nabla \nabla (W^\mu e_\mu) = W^\mu \Tud{\Omega}{\nu}{\mu} \otimes e_\nu = W^\mu \nabla \nabla e_\mu$ This means we only need to check the formula for the frame field instead of arbitrary sections of $E$.
Now, observe that, if $X$ and $Y$ are vector fields, $\begin{align} \nabla_X \nabla_Y e_\mu = \nabla_X ( \langle \Tud{\omega}{\nu}{\mu} , Y \rangle e_\nu) & = \langle \mathrm{d} \langle \Tud{\omega}{\nu}{\mu} , Y \rangle, X \rangle e_\nu + \langle \Tud{\omega}{\nu}{\mu} , Y \rangle \langle \Tud{\omega}{\rho}{\nu} , X \rangle e_\rho \\ \nabla_Y \nabla_X e_\mu = \nabla_Y ( \langle \Tud{\omega}{\nu}{\mu} , Y \rangle e_\nu) & = \langle \mathrm{d} \langle \Tud{\omega}{\nu}{\mu} , X \rangle, Y \rangle e_\nu + \langle \Tud{\omega}{\nu}{\mu} , X \rangle \langle \Tud{\omega}{\rho}{\nu} , Y \rangle e_\rho \\ \nabla_{[X, Y]} e_\mu & = \langle \Tud{\omega}{\nu}{\mu} , [X, Y] \rangle e_\nu \end{align}$ where I have written $\langle - , - \rangle$ for the canonical pairing of a 1-form and a vector field. So we have $\begin{align} (\nabla_X \nabla_Y - \nabla_Y \nabla_X - \nabla_{[X, Y]}) e_\mu & = (\langle \mathrm{d} \langle \Tud{\omega}{\nu}{\mu} , Y \rangle, X \rangle - \langle \mathrm{d} \langle \Tud{\omega}{\nu}{\mu} , X \rangle, Y \rangle - \langle \Tud{\omega}{\nu}{\mu} , [X, Y] \rangle) \, e_\nu \\ & \phantom{=} + (\langle \Tud{\omega}{\nu}{\rho} , X \rangle \langle \Tud{\omega}{\rho}{\mu} , Y \rangle - \langle \Tud{\omega}{\nu}{\rho} , Y \rangle \langle \Tud{\omega}{\rho}{\mu} , X \rangle) \, e_\nu \end{align}$ while on the other hand $\Tud{\Omega}{\nu}{\mu} (X, Y) = \langle \Tud{\Omega}{\nu}{\mu} , X \wedge Y \rangle = \langle \mathrm{d}\Tud{\omega}{\nu}{\mu} , X \wedge Y \rangle + \langle \Tud{\omega}{\nu}{\rho} \wedge \Tud{\omega}{\rho}{\mu} , X \wedge Y \rangle$ but by Cartan's formula for the exterior derivative $\langle \mathrm{d}\Tud{\omega}{\nu}{\mu} , X \wedge Y \rangle = \langle \mathrm{d} \langle \Tud{\omega}{\nu}{\mu} , Y \rangle , X \rangle - \langle \mathrm{d} \langle \Tud{\omega}{\nu}{\mu} , X \rangle , Y \rangle - \langle \Tud{\omega}{\nu}{\mu} , [X, Y] \rangle$ Now, by definition $\langle \Tud{\omega}{\nu}{\rho} \wedge \Tud{\omega}{\rho}{\mu} , X \wedge Y \rangle = \langle \Tud{\omega}{\nu}{\rho} , X \rangle \langle \Tud{\omega}{\rho}{\mu}, Y \rangle - \langle \Tud{\omega}{\nu}{\rho} , Y \rangle \langle \Tud{\omega}{\rho}{\mu}, X \rangle $ and so $\Tud{\Omega}{\nu}{\mu} (X, Y) e_\nu = (\nabla_X \nabla_Y - \nabla_Y \nabla_X - \nabla_{[X, Y]}) e_\mu $ which is exactly what we want.