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.