  A common feature of methods for analyzing samples of probability density
functions is that they respect the geometry inherent to the space of densities.
Once a metric is specified for this space, the Fr\'echet mean is typically used
to quantify and visualize the average density from the sample. For
one-dimensional densities, the Wasserstein metric is popular due to its
theoretical appeal and interpretive value as an optimal transport metric,
leading to the Wasserstein-Fr\'echet mean or barycenter as the mean density. We
extend the existing methodology for samples of densities in two key directions.
First, motivated by applications in neuroimaging, we consider dependent density
data, where a $p$-vector of univariate random densities is observed for each
sampling unit. Second, we introduce a Wasserstein covariance measure and
propose intuitively appealing estimators for both fixed and diverging $p$,
where the latter corresponds to continuously-indexed densities. We also give
theory demonstrating consistency and asymptotic normality, while accounting for
errors introduced in the unavoidable preparatory density estimation step. The
utility of the Wasserstein covariance matrix is demonstrated through
applications to functional connectivity in the brain using functional magnetic
resonance imaging data and to the secular evolution of mortality for various
countries.
