  Ultra-high energy cosmic rays (UHECRs) are atomic nuclei with energies over
ten million times energies accessible to human-made particle accelerators.
Evidence suggests that they originate from relatively nearby extragalactic
sources, but the nature of the sources is unknown. We develop a multilevel
Bayesian framework for assessing association of UHECRs and candidate source
populations, and Markov chain Monte Carlo algorithms for estimating model
parameters and comparing models by computing, via Chib's method, marginal
likelihoods and Bayes factors. We demonstrate the framework by analyzing
measurements of 69 UHECRs observed by the Pierre Auger Observatory (PAO) from
2004-2009, using a volume-complete catalog of 17 local active galactic nuclei
(AGN) out to 15 megaparsecs as candidate sources. An early portion of the data
("period 1," with 14 events) was used by PAO to set an energy cut maximizing
the anisotropy in period 1; the 69 measurements include this "tuned" subset,
and subsequent "untuned" events with energies above the same cutoff. Also,
measurement errors are approximately summarized. These factors are problematic
for independent analyses of PAO data. Within the context of "standard candle"
source models (i.e., with a common isotropic emission rate), and considering
only the 55 untuned events, there is no significant evidence favoring
association of UHECRs with local AGN vs. an isotropic background. The
highest-probability associations are with the two nearest, adjacent AGN,
Centaurus A and NGC 4945. If the association model is adopted, the fraction of
UHECRs that may be associated is likely nonzero but is well below 50%. Our
framework enables estimation of the angular scale for deflection of cosmic rays
by cosmic magnetic fields; relatively modest scales of $\approx\!3^{\circ}$ to
$30^{\circ}$ are favored. Models that assign a large fraction of UHECRs to a
single nearby source (e.g., Centaurus A) are ruled out unless very large
deflection scales are specified a priori, and even then they are disfavored.
However, including the period 1 data alters the conclusions significantly, and
a simulation study supports the idea that the period 1 data are anomalous,
presumably due to the tuning. Accurate and optimal analysis of future data will
likely require more complete disclosure of the data.
