If $f,g : \mathbb R \to \mathbb R$ are smooth, then $(fg)^{\prime \prime} = f^{\prime \prime} g + 2 f^\prime g^\prime + f g^{\prime \prime}$ and, if we define a scalar product on $1$-forms by $\langle \mathbb d x, \mathbb d x \rangle = 1$, this can be rewritten as $\Delta (fg) = f \Delta g \color{red} + 2 \langle \mathbb d f, \mathbb d g \rangle + g \Delta f$, where $\Delta$ is the usual Laplacian.

When trying to obtain the same result on an arbitrary $n$-dimensional Riemannian manifold, I get the opposite sign in front of the middle term, and I do not understand why.

I make the convention that $\Delta = \mathrm d \mathrm d^* + \mathrm d^* \mathrm d$, with $\mathrm d^* \omega = (-1)^{n(k-1)+1} * \mathrm d * \omega$ for any $k$-form $\omega$, and the Hodge star is defined by $\alpha \wedge * \beta = \langle \alpha, \beta \rangle \mathrm {vol}$ for any two $k$-forms $\alpha, \beta$ (and $\mathrm {vol}$ the Riemannian volume form).

With these conventions, we have that $* \mathrm {vol} = 1$ and

$$\Delta (fg) = \mathrm d^* \mathrm d (fg) = - * \mathrm d * (g \mathrm d f + f \mathrm d g) = \\ -* (\mathrm d g \wedge * \mathrm d f + g \mathrm d * \mathrm d f + \mathrm d f \wedge * \mathrm d g + f \mathrm d * \mathrm d g) = \\ - * \langle \mathrm d g, \mathrm d f \rangle \mathrm {vol} + g \Delta f - * \langle \mathrm d f, \mathrm d g \rangle \mathrm {vol} + f \Delta g = \\ f \Delta g \color{red} - 2 \langle \mathrm d f, \mathrm d g \rangle + g \Delta f \ ,$$

and the red "minus" sign preceding the middle term does not agree with the calculation in the first paragraph. What am I doing wrong?

The geometer's convention is, in $\mathbb{R}$, simply to take second derivative. Its spectrum is nonpositive (as witness by trigonometric functions), and in usual $\mathbb{R}^n$ gives (because of parallelizability no further information is encoded in the laplacian for $p$-forms) $\mathrm{div}\circ\mathrm{grad}$ (which also has nonpositive spectrum).
The second $\Delta=(\mathrm{d}+\mathrm{d}^*)^2$ is the analyst laplacian. It has nonnegative spectrum (as you can see from the square, and this is what analysts like because you can do things like log or arbitrary powers and whatnot unambiguously). In $\mathbb{R}^d$ it is $-\mathrm{div}\circ\mathrm{grad}$.
So $(fg)''=f''g+2f'g'+g''f$ translates to
• $\Delta(fg)=(\Delta f)g+2\nabla f\cdot\nabla g+f(\Delta g)$ for the geometer's laplacian; and
• $\Delta(fg)=(\Delta f)g-2\nabla f\cdot\nabla g+f(\Delta g)$ for the analyst's laplacian.