As usual, the trick is to use representation in local coordinates. Note first that using the linearity it is sufficient to show the statement holds for $\omega := f \mathrm d x^I$ for any increasing multi index $I$. Then we have to show that
$$(\mathrm d\delta + \delta \mathrm d)\omega = (\Delta f) \mathrm dx^I.$$
We get
\begin{align*}
\mathrm d \delta \omega &= \mathrm d \left( -\sum_{k=1}^p (-1)^{k-1} \partial_{jk} f \mathrm dx^{i_1} \wedge \dots \wedge \widehat{\mathrm dx^{i_1}} \wedge \dots \mathrm dx^{i_p} \right)\\
&= -\sum_{k=1}^p (-1)^{k-1} \sum_{l=1}^n \partial_l \partial_{jk} f \mathrm dx^{l} \wedge \mathrm dx^{i_1} \wedge \dots \wedge \widehat{\mathrm dx^{i_k}} \wedge \dots \wedge \mathrm dx^{i_p}\\
&= -\sum_{k=1}^p \partial_{jk}^2 f \mathrm dx^I - \sum_{k=1}^p (-1)^{k-1} \sum_{\substack{l=1\\ l \not\in \{i_1, \dots, i_p\}}}^n \partial_l \partial_{jk} f \mathrm dx^{l} \wedge \mathrm dx^{i_1} \wedge \dots \wedge \widehat{\mathrm dx^{i_k}} \wedge \dots \wedge \mathrm dx^{i_p}
\end{align*}
Analogously, but easier, one can show that
\begin{align*}
\delta\mathrm d\omega = -\sum_{k=1}^p \partial_{jk}^2 f \mathrm dx^I + \sum_{k=1}^p (-1)^{k-1} \sum_{\substack{l=1\\ l \not\in \{i_1, \dots, i_p\}}}^n \partial_l \partial_{jk} f \mathrm dx^{l} \wedge \mathrm dx^{i_1} \wedge \dots \wedge \widehat{\mathrm dx^{i_k}} \wedge \dots \wedge \mathrm dx^{i_p}
\end{align*}
Thus, we get
$$ (\mathrm d \delta + \delta\mathrm d) \omega = \left(-\sum_k \partial_k^2 f\right)\mathrm dx^I = (-\Delta f) \mathrm dx^I.$$
PS: Maybe check if a don't mess up with the sign convention, since I usually use the definition $\delta := (-1)^{n(p+1)} * \mathrm d *$.