Factor analysis is a popular tool for identifying and summarizing associations among multiple measures. When measures on organizations, areas, or similar higher-level units are obtained by summarizing data from groups of individuals, associations at the group level are often of primary interest while those at the individual level might not even be meaningfully defined. These data structures can be described with multivariate multilevel models, where the key estimand is the group-level covariance matrix. O’Mallley and Zaslavsky (2008) extended the Fay-Herriott (1979) modeling approach, in which individual-level data for each group are summarized to sufficient statistics (an approximately unbiased point estimate and an estimate of its sampling variance), which become the inputs to the multilevel model. The group-level model is multivariate normal, sometimes with box constraints corresponding to limited-range group measures, making the prior distributions truncated multivariate normal. Sampling from the posterior distribution is facilitated by exploiting a property of the SOUP prior for binomial data (Meng and Zaslavsky 2002). We then apply familiar factor-analysis tools to draws of the covariance matrix, combining Bayesian inference and exploratory data analysis. We illustrate the methods with two motivating examples from health services research, one concerning regional variation in service utilization and the other an analysis of a health care quality survey with structured nonresponse. We conclude by comparing our methods with other approaches to analysis of similarly structured data.