Suppression of the J/ψ nuclear modification factor has long been considered a signature of final state effects in large collision systems. In small systems, nuclear modification was assumed to be due to initial state cold nuclear matter effects, until the observation of strong differential suppression of the ψ(2S) state in p+A collisions suggested the presence of final state effects. Here we present PHENIX results for J/ψ nuclear modification as a function of centrality, p T , and rapidity in p+Al, p+Au, and 3 He+Au collisions at 200 GeV. These measurements are compared with various model predictions, including re-weighted nPDFs, a transport model, and a nuclear absorption model. These measurements offer a new perspective into possible cold nuclear matter effects at RHIC energies, including shadowing, nuclear absorption, and the Cronin effect. This talk will also explore recent PHENIX quarkonia measurements of J/ψ production and polarization as functions of p T and rapidity in p+p collisions at 510 GeV, and present data comparisons at RHIC and LHC energies.