The study is designed to demonstrate and test the utility of the proposed two-stage matching method compared to other analytic methods traditionally employed for multisite observational studies. More specifically, the study addresses the following research questions: (1) How do different specifications of the matching method influence covariate balance? (2) How do different specifications in the matching method influence inferences about treatment effect and effect heterogeneity? The different matching method specifications include differences in the propensity score model and whether a between-site match, within-site match, or two-stage matching process is used. The simulation results indicate that the two-stage matching method balances the desire for within-site covariate balance and the desire to retain as many treatment units in the analysis as possible. Relative to more straightforward matching methods, however, the two-stage matching method does not result in greater covariate balance nor less biased effect estimation. As a result, more straightforward methods that address the nested data structure--such as within-site matching or pooled matching with a random-intercept-and-slope propensity score model--might be preferable to the more complex two-stage matching method. These conclusions are based on a finite set of data generating conditions, with a small set of important confounders at both the unit and site level and a reasonable within-site sample size for matching. Future research should examine the performance of various propensity score model and matching methods under more extreme data conditions. (Contains 2 tables and 5 figures.)