  We study charmonium production in proton-nucleus ($p$-A) collisions focusing
on final-state effects caused by the formation of an expanding medium. Toward
this end, we utilize a rate equation approach within a fireball model as
previously employed for a wide range of heavy-ion collisions, adapted to the
small systems in $p$-A collisions. The initial geometry of the fireball is
taken from a Monte-Carlo event generator where initial anisotropies are caused
by fluctuations. We calculate the centrality and transverse-momentum dependent
nuclear modification factor ($R_{p{\rm A}}$) as well as elliptic flow ($v_2$)
for both $J/\psi$ and $\psi(2S)$ and compare them to experimental data from
RHIC and the LHC. While the $R_{p{\rm A}}$s show an overall fair agreement with
most of the data, the large $v_2$ values observed in $p$-Pb collisions at the
LHC cannot be accounted for in our approach. While the former finding generally
supports the formation of a near thermalized QCD medium in small systems, the
discrepancy in the $v_2$ suggests that its large observed values are unlikely
to be due to the final-state collectivity of the fireball alone.
