  We propose a hierarchical approach to testing general relativity with
multiple gravitational wave detections. Unlike existing strategies, our method
does not assume that parameters quantifying deviations from general relativity
are either common or completely unrelated across all sources. We instead assume
that these parameters follow some underlying distribution, which we parametrize
and constrain. This can be then compared to the distribution expected from
general relativity, i.e. no deviation in any of the events. We demonstrate that
our method is robust to measurement uncertainties and can be applied to
theories of gravity where the parameters beyond general relativity are related
to each other, as generally expected. Our method contains the two extremes of
common and unrelated parameters as limiting cases. We apply the hierarchical
model to the population of 10 binary black hole systems so far detected by LIGO
and Virgo. We do this for a parametrized test of gravitational wave generation,
by modeling the population distribution of beyond-general-relativity parameters
with a Gaussian distribution. We compute the mean and the variance of the
population and show that both are consistent with general relativity for all
parameters we consider. In the best case, we find that the population
properties of the existing binary signals are consistent with general
relativity at the ~1% level. This hierarchical approach subsumes and extends
existing methodologies, and is more robust at revealing potential subtle
deviations from general relativity with increasing number of detections.
