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.