This is an excellent question, which predates the theory of General Relativity by a few years, and has only been astrophysically verified after 1974 (Hulse-Taylor) or even later. As far as I am aware, a complete theoretical answer in General Relativity itself is still elusive.

I'll largely stick with the theory, which I guess is what you are interested in.[1]

The topic is in Part VIII (Chapter 35) of Misner, Thorne & Wheeler's Gravitation ("MTW"), which is the gold standard reference/textbook for General Relativity.

I will later try to briefly summarize the section.

Instead, first, at the root of my answer is the non-linear nature of the Einstein Field Equations, which imply that gravity self-gravitates. This is usually side-stepped by textbooks, which instead proceed to linearize the Einstein Field Equations, that is, they consider the weak limit of gravitation. This is usually to split a metric with a gravitational wave into some minimally-or-even-non-dynamical background and the wave-part, using the former to define the speed of the latter.

For strong gravitational waves, or a very dynamical background, we cannot do this. I think this regime is best studied in a vacuum solution, i.e., where there is only gravitation that self-interacts, although most of the work in this regime appears to have the aim of resolving questions about the distribution of matter in the very early universe (e.g. Misner's mixmaster). An interesting exact solution of the Einstein Field Equations is the https://en.wikipedia.org/wiki/Kasner_metric which can generate singularities and other features formed by gravitational-wave interactionsn. That is, there is manifestly non-linear gravitational self-interaction in the highly-dynamical Kasner chaos. (This is the worst case for the linearized treatments in textbooks).

The Kasner metric can be applied to a Lorentzian manifold (3 space, 1 time dimension), and so is consistent with our universe's causal structure. An interesting feature of this solution is that it engages only two constants in the Einstein Field Equations: c and G. The only velocity scale that we can construct from any combination of these two constants is c itself.

This is highly suggestive that in a universe like ours gravitational radiation must propagate at c.

Again, this is just an argument that there may be an answer within the theory of General Relativity itself, without treating the speed of gravitational waves as a postulate.

This argument is made without regard to the obvious craziness in a Kasner chaos universe compared to our own. General Relativity admits complete gravitational solutions for all sorts of craziness, including universes with any number of space and time dimensions as long as there are at least two total, stress-energy which is distributed very differently from ours (including negative energy, or energy that pops in or out of existence without a cause), and so forth. It is very General. (Special Relativity is very Special: it's defined on -- and only on -- a gravity-free (flat, no gravitational waves) spacetime of exactly 3 spatial and 1 timelike dimension.)

Returning to MTW, they conclude that the propagation speed instead can be no greater than c.

The authors develop an exact vacuum plane-wave solution in §35.9, where the only thing in the spacetime is single large pulse of gravitational radiation in otherwise totally empty flat spacetime, and a set of "test particles", which are well defined probes in General Relativity defined so as to not perturb the solution. They proceed to compare this solution to that of an electromagnetic plane wave in Special Relativity, and arrive at a more physical viewpoint where the gravitational plane wave, if it has anything like a physical source (e.g. a pair of masses in mutual orbit; they return to this in Ch. 36) should be more like a set of "ripples in the spacetime curvature ... propagating on a very slightly curved background spacetime ... The most striking difference between the background and the ripples is not in the magnitude of their spacetime curvatures, but in their characteristic lengths". There is a characteristic length of this background spacetime, determined by its (much much larger) radius of curvature.

They then grind out effective stress-energy tensors, which would couple with any matter in a non-vacuum environment. The argument is that anything in the stress-energy tensor must in certain causal structures (like the one in our universe) must propagate at no more than c.

Their treatment is the basis for a particular type of graviton (Ex 35.16), the scattering and redshift of gravitational waves (.17, .18), and various ways to express them without resort to an effective stress-energy tensor.

In these approaches, c is the limiting speed of gravitational radiation, but gravitational radiation may move slower than c in non-vacuum. In vacuum, the pure pulse may in some circumstances develop a trailing edge that propagates at less than c, but when this can even happen the effect is weak when the wavelength is short compared to the background length scale, or when the amplitude of the pulse is small. Perhaps gravitational astronomy can hope to find someday a high-enough amplitude wave to put this to the test.

For clarity, though, the linearization-is-good-theory claim is on solid footing. LIGO uses the linearized equations and have found agreement between the speed of gravitational waves they've detected to their theory to something around nineteen decimal places. A good multimessenger signal, which seems inevitable, will almost certainly improve that. There is no good reason to expect the speed of gravity in a full solution to the Einstein Field Equations to differ. One can make the same argument about the march of results from post-Newtonian expansions (as below) too. The only "wiggle" room is that cosmic inflation is probably much more dynamical, and the gravitational waves are of much greater amplitude.

- --

[1] The non-theory answer is that the behaviour of orbits in known astrophysical systems which are very post-Newtonian (think black holes, or extremely fast-moving galaxies containing predictable spectra from hydrogen or light curves from supernovae) are consistent with a speed limit on any gravitational interaction, and that the speed limit is very close to c. It turns out to be hard to measure the speed limit exactly. See e.g. Will @ https://arxiv.org/abs/astro-ph/0301145 which discusses light from quasars being gravitationally lensed by Jupiter and how it would look different under theories that admit a propagation speed c_{gravity} different from c_{electromagnetism}. One could also compare bimetric theories of gravitation in which in the early universe gravitational self-interaction propagates differently from electromagnetism, with a view to resolving some questions in the distribution of galaxies in our sky. These generally have to decay the additional metric (meaning gravitational radiation propagates like electromagnetism) in the very early part of the universe or we get something very different from the cosmic web of galaxies that we observe.