One of the big puzzles in quarkonia production is how heavy-quark pairs, produced in the initial hard process, hadronize into final quarkonia states. This uncertainty affects the estimation of quarkonia breakup cross sections and color screening in heavy ion collisions. According to Non-Relativistic Quantum Chromodynamics, quarkonia are formed from several intermediate singlet and color octet states. The relative fraction of each contribution can only be determined by experimental data, including yields and angular decay (polarization) measurements in a broad kinematic region. The PHENIX experiment has measured inclusive J/psi polar and azimuthal angular decay coefficients in the mid ($|y|<0.35$) and forward ($1.2<|y|<2.2$) rapidity regions in $p$+$p$ collisions at 200 GeV and 510 GeV. We find that the polarization coefficients are consistent with zero at mid-rapidity and negative at large rapidity. This talk will present the analysis strategy as well as the results in different angles/frames. We will discuss how the results can be explained in terms of singlet and color octet states.