Here's another approach.
We have
$$\begin{eqnarray*}
\int_0^\infty dx\, \left(\frac{\sin x}{x}\right)^n
&=& \lim_{\epsilon\to 0^+}
\frac{1}{2} \int_{-\infty}^\infty dx\,
\left(\frac{\sin x}{x-i\epsilon}\right)^n \\
&=& \lim_{\epsilon\to 0^+}
\frac{1}{2} \int_{-\infty}^\infty dx\,
\frac{1}{(x-i\epsilon)^n}
\left(\frac{e^{i x}-e^{-i x}}{2i}\right)^n \\
&=& \lim_{\epsilon\to 0^+}
\frac{1}{2} \frac{1}{(2i)^n} \int_{-\infty}^\infty dx\,
\frac{1}{(x-i\epsilon)^n}
\sum_{k=0}^n (-1)^k {n \choose k} e^{i x(n-2k)} \\
&=& \lim_{\epsilon\to 0^+}
\frac{1}{2} \frac{1}{(2i)^n}
\sum_{k=0}^n (-1)^k {n \choose k}
\int_{-\infty}^\infty dx\, \frac{e^{i x(n-2k)}}{(x-i\epsilon)^n}.
\end{eqnarray*}$$
If $n-2k \ge 0$ we close the contour in the upper half-plane and pick up the residue at $x=i\epsilon$.
Otherwise we close the contour in the lower half-plane and pick up no residues.
The upper limit of the sum is thus $\lfloor n/2\rfloor$.
Therefore, using the Cauchy differentiation formula, we find
$$\begin{eqnarray*}
\int_0^\infty dx\, \left(\frac{\sin x}{x}\right)^n
&=& \frac{1}{2} \frac{1}{(2i)^n}
\sum_{k=0}^{\lfloor n/2\rfloor} (-1)^k {n \choose k}
\frac{2\pi i}{(n-1)!}
\left.\frac{d^{n-1}}{d x^{n-1}} e^{i x(n-2k)}\right|_{x=0} \\
&=& \frac{1}{2} \frac{1}{(2i)^n}
\sum_{k=0}^{\lfloor n/2\rfloor}
(-1)^k {n \choose k}
\frac{2\pi i}{(n-1)!} (i(n-2k))^{n-1} \\
&=& \frac{\pi}{2^n (n-1)!}
\sum_{k=0}^{\lfloor n/2\rfloor} (-1)^k {n \choose k} (n-2k)^{n-1}.
\end{eqnarray*}$$
The sum can be written in terms of the hypergeometric function but the result is not particularly enlightening.
It can be solved using differentiation under the integral sign. Consider the following integral:
\begin{equation}
I(t)=\int\limits_{-\infty}^{+\infty} \frac{\cos(tx)}{x^{2}+k} \,dx = 2\int\limits_{0}^{+\infty} \frac{\cos(tx)}{x^{2}+k} \,dx
\end{equation}
for any positive real $t$ and $k$. The first derivative with respect to $t$ is:
\begin{equation}
I'(t)= -2\int\limits_{0}^{+\infty} \frac{x\sin(tx)}{x^{2}+k} \,dx
\end{equation}
\begin{equation}
\Leftrightarrow \hspace{.3cm}I'(t)= -2\int\limits_{0}^{+\infty} \frac{x^{2}\sin(tx)}{x(x^{2}+k)} \,dx
\end{equation}
\begin{equation}
\Leftrightarrow \hspace{.3cm}I'(t)= -2\int\limits_{0}^{+\infty} \frac{(x^{2}+k-k)\sin(tx)}{x(x^{2}+k)} \,dx
\end{equation}
\begin{equation}
\Leftrightarrow \hspace{.3cm}I'(t)= -2\int\limits_{0}^{+\infty} \frac{\sin(tx)}{x} \,dx +2k\int\limits_{0}^{+\infty} \frac{\sin(tx)}{x(x^{2}+k)} \,dx
\end{equation}
The first one is just the sine integral as $x\rightarrow \infty$ and it is known to converge to $\frac{\pi}{2}$. Thus:
\begin{equation}
I'(t)= 2k\int\limits_{0}^{+\infty} \frac{\sin(tx)}{x(x^{2}+k)} \,dx -\pi
\end{equation}
Differentiating once more with respect to $t$ yields:
\begin{equation}
I''(t)= 2k\int\limits_{0}^{+\infty} \frac{\cos(tx)}{x^{2}+k} \,dx
\end{equation}
\begin{equation}
\Leftrightarrow \hspace{.3cm}I''(t)-kI(t)=0
\end{equation}
The general solution to the ODE is:
\begin{equation}
I(t)=c_{1}e^{\sqrt{k}t}+c_{2}e^{-\sqrt{k}t}
\end{equation}
Plugging some conditions $\left(I(t=0) \,\,\text{and}\,\, I'(t=0)\right)$ allows you to find that $c_{1}=0$ and that $c_{2}=\frac{\pi}{\sqrt{k}}$. Then:
\begin{equation}
\boxed{\int\limits_{-\infty}^{+\infty} \frac{\cos(tx)}{x^{2}+k} \,dx = \frac{\pi}{\sqrt{k}}e^{-\sqrt{k}t}}
\end{equation}
for positive real values of $t$ and $k$. If you plug $t=2$ and $k=4$, you obtain the desired result.
Best Answer
Here is a partial answer: If $n \geq 2$ is an integer, then
$$ I_n = \frac{1}{(n-1)!2^{n-1}} \sum_{l=0}^{\lfloor n/2 \rfloor} (-1)^l \binom{n}{l} (n-2l)^{n-1} J_{n-2l}, \tag{1} $$
where $J_p$ is defined by
$$ J_{p} = \begin{cases} \displaystyle \frac{4p}{\pi} \sum_{j=1}^{\infty} \frac{\log(2\pi j)}{4j^2-p^2}, & \text{if $p$ is odd}, \\ \displaystyle \frac{\pi}{2}, & \text{if $p$ is even}. \end{cases} $$
Proof of $\text{(1)}$. The case of even $n$ has already been discussed in other postings, so we focus on odd $n$. We first note that, if $n \geq 1$ is an odd integer, then
\begin{align*} \frac{\mathrm{d}^{n-1}}{\mathrm{d}x^{n-1}} \sin^n x &= \frac{1}{(2i)^n} \sum_{l=0}^{\frac{n-1}{2}} (-1)^l \binom{n}{l} \frac{\mathrm{d}^{n-1}}{\mathrm{d}x^{n-1}} (e^{(n-2l)ix} - e^{-(n-2l)ix}) \\ &= \frac{1}{2^{n-1}} \sum_{l=0}^{\frac{n-1}{2}} (-1)^l \binom{n}{l} (n-2l)^{n-1} \sin((n-2l)x). \end{align*}
So by applying integration by parts $(n-1)$-times, we get
\begin{align*} I_n &= \sum_{k=0}^{\infty} \int_{0}^{\pi} \frac{\sin^n x}{(x+k\pi)^n} \, \mathrm{d}x \\ &= \frac{1}{(n-1)!} \sum_{k=0}^{\infty} \int_{0}^{\pi} \biggl( \frac{1}{x+k\pi} - \frac{1}{(k+1)\pi} \biggr) \biggl( \frac{\mathrm{d}^{n-1}}{\mathrm{d}x^{n-1}} \sin^n x \biggr) \, \mathrm{d}x \\ &= \frac{1}{(n-1)!2^{n-1}} \sum_{l=0}^{\frac{n-1}{2}} (-1)^l \binom{n}{l} (n-2l)^{n-1} J_{n-2l}, \end{align*}
where $J_{p}$ is defined by
\begin{align*} J_{p} &= \sum_{k=0}^{\infty} \int_{0}^{\pi} \biggl( \frac{1}{x+k\pi} - \frac{1}{(k+1)\pi} \biggr) \sin(px) \, \mathrm{d}x. \end{align*}
If $p$ is odd, then the above definition is recast as
\begin{align*} J_{p} &= \sum_{k=0}^{\infty} \biggl( \int_{0}^{\pi} \frac{1}{x+k\pi} \sin(px) \, \mathrm{d}x - \frac{2}{p\pi(k+1)} \biggr) \\ &= \lim_{N \to \infty} \biggl( \int_{0}^{N \pi} \frac{\sin(p(x \text{ mod } \pi))}{x} \, \mathrm{d}x - \frac{2}{p\pi} H_N \biggr), \end{align*}
where $H_N = 1 + \frac{1}{2} + \dots + \frac{1}{N}$ is the $N$-th harmonic number. Still assuming that $p$ is an odd integer, Fourier series computation shows that
\begin{align*} \sin(p(x \text{ mod } \pi)) &= \frac{2}{p\pi} - \frac{4p}{\pi} \sum_{n=1}^{\infty} \frac{\cos(2\pi n x)}{4n^2-p^2} \\ &= \frac{4p}{\pi} \sum_{j=1}^{\infty} \frac{1 - \cos(2\pi j x)}{4j^2-p^2}, \end{align*}
and so,
\begin{align*} \int_{0}^{N \pi} \frac{\sin(p(x \mathrm{ mod } \pi))}{x} \, \mathrm{d}x &= \frac{4p}{\pi} \sum_{j=1}^{\infty} \frac{1}{4j^2-p^2} \int_{0}^{N \pi} \frac{1 - \cos(2 j x)}{x} \, \mathrm{d}x \\ &= \frac{4p}{\pi} \sum_{j=1}^{\infty} \frac{1}{4j^2-p^2} (\gamma + \log(2\pi j N) - \operatorname{Ci}(2\pi j N) ). \end{align*}
Plugging this back and using the identity $\frac{4p}{\pi} \sum_{j=1}^{\infty} \frac{1}{4j^2-p^2} = \frac{2}{p\pi}$, which itself follows from the Fourier series of $\sin(p(x \text{ mod } \pi))$, we finally obtain
\begin{align*} J_{p} &= \frac{4p}{\pi} \lim_{N \to \infty} \sum_{j=1}^{\infty} \frac{1}{4j^2-p^2} (\gamma + \log(2\pi j N) - \operatorname{Ci}(2\pi j N) - H_N ) \biggr) \\ &= \frac{4p}{\pi} \sum_{j=1}^{\infty} \frac{\log(2\pi j)}{4j^2-p^2} \end{align*}
as desired. $\square$
Addendum. Here is a Mathematica code for numerical verification of $\text{(1)}$: