Edit: As mentioned in the comments, the Laplacian on Riemannian manifolds is defined not by $\Delta=*d*d+d*d*$, but $\Delta=\delta d+d\delta$, where the coderivative $\delta$ is defined by $\delta=(-1)^{nk+n+1}*d*$ for a $k$-form in $n$ dimensions. In 3D, this is the same as $\delta=(-1)^k*d*$.
I want to point out an even simpler case which perhaps highlights where you are going wrong.
The laplacian evaluated on $f(x)$ (Note that in the next few calculations the $(-1)^k$ is symbolic and is acting as an operator, i.e. noncommutative, because I am waiting to see what rank the form is before I evaluate it):
$$\Delta f=(-1)^k*d*df + d(-1)^k*d*f=(-1)^k*d*\partial_xf\,dx + d(-1)^k*df\,dx\wedge dy\wedge dz$$
$$=(-1)^k*d\partial_xf\,dy\wedge dz + d(-1)^k*\partial_xf\,dx\wedge dx\wedge dy\wedge dz$$
$$=(-1)^k*\partial^2_xf\,dx\wedge dy\wedge dz=(-1)^0\partial^2_xf=\partial^2_xf$$
Nothing fancy here, just what you'd expect. Now let's try $f(x)\,dy$:
$$\Delta f\,dy=(-1)^k*d*df\,dy + d(-1)^k*d*f\,dy$$
$$=(-1)^k*d*\partial_xf\,dx\wedge dy - d(-1)^k*df\,dx\wedge dz$$
$$=(-1)^k*d\partial_xf\,dz + d(-1)^k*\partial_xf\,dx\wedge dx\wedge dz$$
$$=(-1)^k*\partial^2_xf\,dx\wedge dz=-(-1)^1\partial^2_xf\,dy=\partial^2_xf\,dy$$
So the missing minus sign is coming from this mysterious "codifferential" object. Let's try the full calculation now:
$$\Delta(fg\,dy)=((-1)^k*d*d)(fg\,dy) + (d(-1)^k*d*)(fg\,dy)$$
$$=((-1)^k*d*)(g\partial_xf\,dx\wedge dy +f\partial_zg\,dz\wedge dy) + (d(-1)^k*d)(fg\,dz\wedge dx)$$
$$=((-1)^k*d)(g\partial_xf\,dz-f\partial_zg\,dx) + (d(-1)^k*)(f\partial_yg\,dy\wedge dz\wedge dx)$$
$$=(-1)^1*(g\partial^2_xf\,dx\wedge dz + \partial_xf \partial_yg\,dy\wedge dz-f\partial_y\partial_zg\,dy\wedge dx-f\partial^2_zg\,dz\wedge dx) + d(-1)^0(f\partial_yg)$$
$$=g\partial^2_xf\,dy - \partial_xf \partial_yg\,dx-f\partial_y\partial_zg\,dz+f\partial^2_zg\,dy+\partial_xf\partial_yg\,dx+f\partial^2_yg\,dy+f\partial_z\partial_yg\,dz$$
$$=(\partial^2_xf)g\,dy +f(\partial^2_y+\partial^2_z)g\,dy$$
and all is well.
I think the answers to your first two questions are yes, there is always a codifferential; and no, there is not always a Hodge decomposition. I don't know anything about the Einstein manifold case.
Let's say $E$ is a vector bundle over $M$ with metric and compatible connection $\nabla$. ($E$ could be some tensor bundle with $\nabla$ induced from the Levi-Civita connection, for example.) As you say, $\nabla$ gives an exterior derivative $d^\nabla$ on $E$-valued differential forms that obeys the Leibniz rule
$$ d^\nabla (\omega \otimes e) = d\omega \otimes e + (-1)^{\text{deg} \omega} \omega \wedge \nabla e ,$$
where $\omega$ is a homogeneous form and $e$ is a section of $E$, and $\omega \wedge \nabla e$ means $\sum_i (\omega \wedge dx^i) \otimes \nabla_{\partial_i} e$.
For each degree of forms $k$, you can view $d^\nabla$ as a map $d^\nabla: C^\infty( \Lambda^k \otimes E) \to C^\infty( \Lambda^{k+1} \otimes E)$. Using the metrics, each of those spaces of sections are endowed with $L^2$ inner products. Then $d^\nabla$ always has a formal $L^2$-adjoint $\delta^\nabla$ going the other way: $\delta^\nabla: C^\infty( \Lambda^{k+1} \otimes E) \to C^\infty( \Lambda^k \otimes E)$.
(Note the space of smooth sections can be completed to the space of $L^2$ sections, but $d^\nabla$ and $\delta^\nabla$ are unbounded operators and can generally only be defined on some dense subspace of that $L^2$ space.)
So you can always define $d^\nabla$ and its formal adjoint $\delta^\nabla$. The issue is that $d^\nabla$ squares to zero if and only if the connection $\nabla$ on $E$ is flat, meaning its curvature is zero. (Flat vector bundles $E$, i.e., bundles on which there exists a flat connection, are relatively scarce.) We need $d^\nabla$ to square to zero to even attempt to define the de Rham cohomology of $E$-valued forms, and to have nice properties like the Hodge decomposition. (I want to say that the "Dirac operator" $D = d^\nabla + \delta^\nabla$ and the "Laplacian" $\Delta = D^2$ are not elliptic unless $d^\nabla$ squares to zero, but I need to think about that.)
I don't have a good reference for this at hand, although the Wikipedia article on bundle-valued forms might be useful.
Best Answer
You are mixing the geometer's convention with the analyst's convention.
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