Why do the Einstein field equations (EFE) involve the Ricci curvature tensor instead of Riemann curvature tensor?

I think this question is more trivial than you think.

You should ask yourself why should the full Riemann tensor appear. I'll sketch a heuristic derivation of the field equations.

We know that with small velocities and a static field, the Poisson equation $$\Delta\phi=4\pi G\rho$$ is approximately satisfied. From special relativity we know that the mass/energy density $\rho$ must change with two Lorentz factors under a Lorentz transformation. Thus it is the time-time component of a rank two tensor $T_{\mu\nu}$. Using the equivalence principle, we promote this to a curved spacetime tensor. When we look for field equations, we demand that they be tensor equations. For one thing, this means we must have the same number of indices on both sides. We posit $$D_{\mu\nu}=\kappa T_{\mu\nu}$$ with $\kappa$ some constant. We don't know what $D_{\mu\nu}$ is, but the principle of covariant conservation fixes it to be the Einstein tensor. Note that the general form is the natural generalization of the Poisson equation.

You might propose an equation with more indices, such as $$R_{\mu\nu\rho\sigma}=\kappa'T_{\mu\nu}T_{\rho\sigma}$$ with some appropriate antisymmetrization scheme. What are the vacuum equations? They would be $$R_{\mu\nu\rho\sigma}=0$$ But this just says spacetime is flat! We know this is incorrect. Black holes are certainly vacuum solutions but are also certainly not flat spacetime solutions.

In summary, the Ricci tensor has the ability to vanish without the full Riemann tensor vanishing. The general form of the equations is determined by the Poisson equation to be a rank two equation. In my mind, these two facts are the most effective argument.

Regarding why the Ricci tensor, not the Riemann, appears in the EFEs the answer is in the Newtonian theory.Very heuristically consider Newton's law

$\ddot{r}=-\frac{GM}{r²} $.

Now let's try to make this a more local statement. In order to do so I'll divide both sides by $r$ and some numerical factors in order for the volume $V=\frac{4\pi r³}{3}$ of some hypothetical sphere enclosing the mass to appear in the denominator

$\frac{\ddot{r}}{r}= -\frac{4\pi G}{3} \left(\frac{M}{\frac{4\pi r³}{3}}\right)=-\frac{8\pi G}{3}\left(\frac{1}{2}\rho\right)$,

where $\rho$ is the matter density. Now use the fact that

$\frac{\ddot{V}}{V}=3\frac{\ddot{r}}{r}+O(V^{-2})$,

to write Newton's law, in a first approximation, as

$\frac{\ddot{V}}{V}=-8\pi G\left(\frac{1}{2}\rho\right)$,

so that Newton's law can be interpreted as relating the fractional volume change with the matter distribution. Now from the point of view of General Relativity gravity is curvature, which is described by the Riemann tensor. In order to obtain Newton's law you could ask for the part of the Riemann curvature that describes volume changes, and will not be surprised to find that it is the Ricci tensor (wikipedia has a section on this interpretation of Ricci curvature). The other part, the Weyl tensor, describes shape deforming, volume preserving curvature, and therefore cannot be directly related to Newton's law.

When writing the EFEs as

$R_{\mu\nu}=8\pi G (T_{\mu\nu}-1/2g_{\mu\nu}T^\mu_{\;\mu})$,

and taking the time-time component you get on the left hand side the fractional change in volume (plus the $O(V^{-2})$ term that I neglected), and the right hand side provides the one-half density (plus terms of pressure divided by $c²$ which don't appear in the newtonian limit).

Regarding the Riemann tensor, it is true that it contains more information than the Ricci tensor, but this is not irrelevant. In fact the Riemann curvature satisfies the second Bianchi identity $R_{\alpha\beta[\mu\nu;\lambda]}=0$, and if you're masochistic enough to decompose the Riemann tensor in Ricci and Weyl parts (here in wikipedia for the expression) you'll see that one can write the derivatives of the Weyl part in terms of the Ricci tensor and its derivative. Now use the EFEs to substitute the Ricci tensor for the energy-momentum and the second Bianchi indentity will give you a system of partial differential equations relating the Weyl components and the energy-momentum (and it's derivatives). So, in a certain sense, the curvature equations of general relativity are the Einstein Equations plus the Bianchi indentity, that determines the whole curvature. In the absence of matter the manifold will be Ricci flat (by EFEs) but the Weyl part need not be null, the solution to the system of PDEs will be determined by the boundary conditions. In relativity this boundary conditions are things like the spacetime being asymptotically flat, symmetries, etc.

This situation is analogous to electromagnetism. One can write Maxwell's equations covariantly as $\partial_\alpha F^{\alpha\beta}=J^\beta$ and $\partial_{[\gamma}F_{\alpha\beta]}=0$. But you could write only the first equation (the one with the sources) and just define the Faraday tensor to be a closed form, which gives you the second equation, Gauss' law and Faraday's law. Since the Bianchi identity is, well, an identity, the equations for the Weyl part are already determined by the definition of Riemann curvature.

I hope this clarifies the presence of the Ricci tensor in the EFEs and how the rest of the Riemann curvature is obtained.

EDIT: Per Rod Vance's suggestion I'll put the Bianchi identities explicitly in terms of Weyl and Ricci parts. Straight out of the wikipedia link he furnished below, the Bianchi identities read (for our four dimensional case):

$\nabla_\mu C^\mu_{\; \nu\lambda\sigma}=\nabla_{[\lambda}( R_{\sigma]\mu}-\frac{R}{3}g_{\sigma]\mu})$,

where brackets denote antisymmetrization as usual, and one should note that although the derivative of the metric is zero, one must apply Leibniz rule here to take in account derivatives of the scalar curvature as well. Therefore in the presence of matter the Weyl tensor will describe the tidal forces (differential gravity) sourced by this objects, but even in the absence of matter one will still have to solve a PDE to find $C^\mu_{\;\nu\lambda\sigma}$, the solution of which is subject to the aforementioned boundary conditions. Gravitational waves in particular are vacuum solutions where there are gravitational tidal forces in a given region of the spacetime even in the absence of matter.

The Einstein equations are equivalent to an relation on the Riemann tensor. Given $a, b$ that are linearly independent vectors,

$$R(a \wedge b) = C(a \wedge b) + 4\pi [a \wedge T(b) - b \wedge T(a) - \frac{2}{3} T a \wedge b]$$

where $T(a)$ is the stress energy tensor acting on $a$ and $T$ is its trace, and $C$ is the Weyl ("conformal") tensor.

This equation subsumes the Einstein equations, and it forms part of the basis for numerical tetrad approaches to gravity.

I say this equation subsumes the Einstein equations because it requires strictly more information than the stress-energy tensor to characterize the Riemann this way. The Weyl tensor describes the state of gravitational radiation in the system, and so you can see the Riemann has two distinct contributions: one from stress-energy, and one from gravitational radiation.

The symmetries of the Weyl tensor mean that you can contract both sides and eliminate it, yielding the familiar Einstein equations in terms of the Ricci tensor. This can be convenient, as often we might treat different systems with the same distribution of stress energy--but different distributions of gravitational radiation--as "equivalent" in some sense.

Compare with common electrodynamics problems: e.g. solving for the electric field from a spherical charge density. Technically, you could add any divergence- and curl-free electric field and still satisfy the Maxwell equations, but everybody knows that and it's tacitly not considered except when that's really relevant.