A popular method for estimating a causal treatment effect with observational data is the difference-in-differences (DiD) model. In this work, we extend the classical DiD model to the hierarchical context in which data cannot be matched at the most granular level. Our motivating example is an application to assess the impact of primary care redesign policy on diabetes outcomes in Minnesota, in which patient-level outcomes are not matched longitudinally, and thus the mean change in outcome, treatment(s), and covariates are measured at the clinic level. We propose a Bayesian hierarchical difference-in-differences (HDiD) model which estimates the treatment effect by regressing the treatment on a latent variable representing the mean change in group-level outcome. We first apply the HDiD model to measure the impact of primary care redesign on clinics in Minnesota from 2008 to 2012. We go on to present theoretical and empirical results showing that an HDiD model that fails to adjust for a particular class of confounding variables biases the treatment effect estimate. Motivated by the need for covariate adjustment, we propose and implement various approaches to perform variable selection using a structured Bayesian spike-and-slab model in the HDiD context, evaluating their properties with theoretical results and through simulation. We then conduct a cross-sectional analysis to measure the specific primary care services and resources which most strongly relate with improved diabetes outcomes in 2017. We conclude by introducing an additional variable selection approach which leverages the temporality of the HDiD context, as well as any structure in the predictors.