Are the matrix elements of $S$-matrix Lorentz invariant?

The whole point of QFT is that it is a framework that allows you to define Lorentz (co)variant scattering amplitudes. In fact, under some general hypothesis it is the only framework with that property. The expression in the OP is not manifestly Lorentz covariant, although it turns out to be, after a very cumbersome analysis. See ref.1 for a detailed discussion.

Alternatively, the Lorentz covariance of the scattering amplitudes follows immediately from the Lorentz invariance of transition amplitudes, say, as given by the LSZ formula, e.g. $$ \langle p_{1},\ldots ,p_{n}\ {\mathrm {out}}|q_{1},\ldots ,q_{m}\ {\mathrm {in}}\rangle =\int \prod _{{i=1}}^{{m}}\left\{{\mathrm {d}}^{4}x_{i}{\frac {ie^{{-iq_{i}\cdot x_{i}}}\left(\Box _{{x_{i}}}+m^{2}\right)}{(2\pi )^{{{\frac {3}{2}}}}Z^{{{\frac {1}{2}}}}}}\right\}\prod _{{j=1}}^{{n}}\left\{{\mathrm {d}}^{4}y_{j}{\frac {ie^{{ip_{j}\cdot y_{j}}}\left(\Box _{{y_{j}}}+m^{2}\right)}{(2\pi )^{{{\frac {3}{2}}}}Z^{{{\frac {1}{2}}}}}}\right\}\langle 0|{\mathrm {T}}\varphi (x_{1})\ldots \varphi (x_{m})\varphi (y_{1})\ldots \varphi (y_{n})|0\rangle $$ for a scalar field theory in four spacetime dimensions. In turns, the correlation function $\langle 0|{\mathrm {T}}\varphi (x_{1})\ldots \varphi (y_{n})|0\rangle$ is Lorentz invariant, which is manifest when written as a path integral, $$ \langle 0|{\mathrm {T}}\varphi (x_{1})\ldots \varphi (y_{n})|0\rangle\sim \int \varphi (x_{1})\ldots \varphi (y_{n})\ e^{-I[\phi]}\ \mathrm d\varphi $$ where $I$ is the action, a Lorentz scalar.

All the elements in the above formulas are manifestly Lorentz covariant, and so it immediately follows that so are the $S$-matrix elements. The case of arbitrary spin is handled similarly, although the details are more involved. Again, see ref.1 for the explicit construction.

FWIW, the $S$-matrix can also be written directly as a path integral, sidestepping the need of ever introducing the LSZ formula. See e.g. ref.2, where we can read, for example, (p.11)

Comparing [...] with [...] we can say that to get the $S$-matrix in the momentum representation one is to calculate the Feynman functional integral [...]. This observation shows that the $S$-matrix is defined by means of a functional integral in a more elegant manner than the evolution operator. [...] This is specially attractive in the case of field theory because [it] becomes manifestly Lorentz covariant [...].

References.

1. Weinberg - Quantum theory of fields, Vol.1. Foundations.

2. Faddeev - Introduction to Functional Methods, from Methods in Field Theory, edited by North Holland. See also this PSE post.

1)
A normalization condition for one-particle states $| \vec p \rangle$, which is Lorentz invariant is
$\langle \vec p | \vec q \rangle = 2 E_p (2 \pi)^3 \delta^{(3)} (\vec p - \vec q)$
(The factor $2$ is conventional)

However, if we write
$\delta^{(3)} (\vec k) = \int \frac{d^3x}{(2 \pi)^3} exp (i \vec k \cdot \vec x)$
we can state
$\delta^{(3)} (\vec 0) = \frac{V}{(2 \pi)^3}$
Hence
$\frac{\delta^{(3)} (\vec 0')}{V'} = \frac{\delta^{(3)} (\vec 0)}{V}$

So, the normalization in your post is invariant as well.

2)
The S-matrix element $\langle f|S| i \rangle$ for $n$ asymptotic momentum eigenstates is given by the LSZ (Lehmann-Symanzik-Zimmermann) reduction formula which for scalar quantum fields $\phi (x)$ states
$\langle f|S| i \rangle = [i \int d^4 x_1 exp (-i p_1 x_1) (\Box_1 + m^2)] \cdot \cdot \cdot [i \int d^4 x_n exp (i p_n x_n) (\Box_n + m^2)] \langle \Omega | T \{ \phi (x_1) \cdot \cdot \cdot \phi (x_n)\} | \Omega \rangle$
where $-i$ in the exponent applies to initial states and $+i$ to final states

The LSZ formula is constructed from Lorentz covariant fields, hence the S-matrix element is an invariant.