In special relativity the 4-tuple of numbers, energy and 3 components of momentum, form a 4-vector with respect to Lorentz transformations.
Is it true that the analogous 4-tuple of energy and momentum of electromagnetic field in vacuum is also a 4-vector with respect to Lorentz transformation?
The answer depends on whether the electromagnetic field in question is interacting with other things. Here is the sequence of logic.
energy-momentum of a single particle is a 4-vector. It transforms as a 4-vector.
the total energy-momentum of a collection of particles in a collision all at the same place is a 4-vector
the total energy-momentum of a set of particles located at different places is a 4-vector as long as the set of particles is isolated. But if they are not isolated then the total energy-momentum may be changing, and the timing of this change is such that the 4-vector sum will depend on which set of events are deemed simultaneous. This sum is then not guaranteed to transform as the Lorentz transformation, so it is not a 4-vector.
an electromagnetic field is an extended object, like a collection of particles at different places. Its total energy-momentum is a 4-vector if the field is not interacting with other things. For example a light pulse in empty space is isolated in the right way and so its total energy-momentum is a 4-vector. But the total energy-momentum of an e.m. field interacting with charges, such as the static field around a charged sphere for example, is not a 4-vector, because such a field is not an isolated system.
The analysis of this aspect of electromagnetism is associated with a problem called the "four thirds problem", in which different ways to find the energy or momentum of the field of a charged object disagree by a factor 4/3. The reason is that one has to take into account the forces that are keeping the charged object in equilibrium, whatever they may be, and the meeting of these forces at the charge is providing a continuous exchange of momentum between the e.m. field and these other forces. The other forces may, for example, be internal forces in a charged sphere which keep it in mechanical equilibrium, or the forces in the separator which holds the plates of a capacitor apart. They are called Poincare stresses.
Yes, it is true. If $\Sigma$ is any Minkowskian time slice and $T_{EM}^{\mu\nu}$ is the stress-energy tensor of the EM field (assumed to be an isolated system), then $$P^\mu = \int_\Sigma T_{EM}^{\mu0} dx^1dx^2dx^3$$ (assuming the integral converges) turns out to be a four vector and has the meaning you said.
Notice that the four vector $P^\mu$ does not depend on the Minkowski reference frame whose $\Sigma$ is a spatial section. Actually all that is valid for every isolated field or isolated continuous system that vanishes sufficiently fast at spatial infinity.
Here is the shortest proof I am able to produce.
So assume that, in Minkowski spacetime, $$\nabla_\nu T^{\mu \nu}=0\tag{1}$$ (isolated system!).
Consider two time slices $\Sigma$ and $\Sigma'$ associated to a pair of generally different future-oriented Minkowskian reference frames $x^0,x^1,x^2,x^3$ and $x'^0,x'^1,x'^2,x'^3$ respectively.
Finally, consider a $4D$ cylinder $C$ with bases in $\Sigma$ and $\Sigma'$ and lateral surface of timelike type.
We henceforth consider the portion of the system contained in the cylinder and we suppose that this lateral surface of $C$ is sufficiently close to the spatial infinity that the system vanishes outside $C$ between $\Sigma$ and $\Sigma'$.
If $\omega$ is a constant co-vector, (1) implies $$\nabla_\nu (\omega_\mu T^{\mu\nu})=0\:.$$ The divergence theorem applied to $C$, dropping the contribution of the lateral surface leads to $$\int_{\Sigma} \omega_\mu T^{\mu0} dx^1dx^2dx^3= \int_{\Sigma'} \omega'_\alpha T'^{\alpha 0} dx'^1dx'^2dx'^3\:.$$ This means that the map $$\omega \mapsto \int_{\Sigma} \omega_\mu T^{\mu0} dx^1dx^2dx^3 = \int_{\Sigma'} \omega'_\alpha T'^{\alpha 0} dx'^1dx'^2dx'^3$$ does not depend of the used Minkowski reference frame nor on the instants of time used to define $\Sigma$ and $\Sigma'$ in each reference frame.
The crucial observation is now that the above map
(a) is linear
(b) produces real numbers.
A linear map form the space of co-vectors to $\mathbb{R}$ is an element to the dual space of co-vectors, i.e., it is a vector!
We conclude that ther must exist a constant four vector $P$, such that $$\omega_\mu P^\mu =\int_{\Sigma} \omega_\mu T^{\mu0} dx^1dx^2dx^3 = \int_{\Sigma'} \omega'_\alpha T'^{\alpha 0} dx'^1dx'^2dx'^3 = \omega'_\alpha P'^\alpha\:,$$ for every co-vector $\omega$.
Since $\omega$ is constant, it can be carried out of the integration symbol. Its arbitrariness eventually implies that $$P^\mu = \int_{\Sigma} T^{\mu0} dx^1dx^2dx^3$$ defines a four-vector (the same independently of the choice of the reference frame). In particular, $$P'^\alpha = \Lambda^\alpha_\mu P^\mu \tag{2}$$ when $$x'^\alpha = \Lambda^\alpha_\mu x^\mu + b^\alpha\:.$$
As a final remark, I stress that if we use the same reference frame so that $\Sigma$ and $\Sigma'$ are simply the rest spaces at different times, independence of $P^\mu$ from the time labeling the time slice just says that the component of the four-momentum are separately constant in time in that (generic) reference system. Energy and momentum are separately conserved in every fixed inertial reference frame.
ADDENDUM: the simplest example of a consistent EM field + charges model.
Coming back to the case of the EM field, when we add the source of the field to the picture and we consider the total system
EM field + charges,
the dynamically conserved stress energy tensor is the total stress-energy tensor of the system. The conserved total quantities should be defined accordingly.
As discussed in the last item of @Andrew Steane's excellent answer, a classical model of a finite-radius elementary electric charge gives rise to the insurmountable problem of 4/3 factor. This happens when one tries to interpret the mass and the momentum of an elementary charge as of completely electromagnetic nature. One could try to add unknown Poincaré's stresses without solving the physical problem (and the zero-radius model leads to the even worst issue of an infinite self-force). A posteriori, all that seems to be related to the physical fact that an electron is a quantum particle and any classical description is hopeless from scratch. However the quantum problem is affected by similar issues, those of the renormalization
In any cases, effective macroscopic models, for charged fluids or gases, can be constructed an used for the total system EM field + charges, where the charges are those of a macroscopic system. The simplest one is the model of rarefied gas where the internal stresses of the gas can be approximatively disregarded. This is nothing but the perfect-fluid stress-energy tensor with zero pressure. (The model can be extended to the presence of pressure and macroscopic stresses for a concrete macroscopic fluid.) Its stress energy tensor reads $$T_M^{ab} := \mu V^aV^b\:.$$
Above, (a) $V$ is the field of four-velocities,
(b) $\mu \geq 0$ is the density of mass (mesured at rest with an integral line of $V$.
(c) We also assume that the gas is charged with a density of electric charge $$\rho = k \mu$$ for some constant $k$ (it means that the microsopical charge carriers are identical).
(d) We eventaully associate a current $J^a = \rho V^a$ to the gas.
In this situation, the total stress-energy tensor is $$T^{ab} = T_M^{ab}+ T_{EM}^{ab}$$ (in this concrete case no third interaction term shows up on the right-hand side in view of the physical simplicity of the system).
We should also assume that the Maxwell equations are satisfied and that the field is completely generated by the said charges. $$\nabla_a F^{ab} = -J^b\:; \quad \epsilon^{abcd}\nabla_b F_{cd}=0\:.$$ Notice that from the former, since $F$ is antysimmetric, we infer both the conservation of the charge and the one of the mass ($k$ is constant!), simply taking the divergence of both sides.
A quite lengthy computation proves that $$\nabla_a T^{ab} =0$$ is equivalent to $$\nabla_a T_M^{ab} = F_L^b\:, \quad F_L^b= - J_a F^{ab}$$ the latter density of four-force being just the Lorentz one.
In turn, the former equation, taking the conservation of the mass into account, is equivalent to $$\mu V^a\nabla_a V^b = F^b_L\:.$$ This identity asserts that the integral lines of the field of four-velocities evolves according to the presence of the Lorentz force. Notice that, in the absence of the EM field, we would instead have the geodesic equation $$V^a\nabla_a V^b =0\:.$$
