Recent results from PHENIX experiment are presented especially on collective flow in small systems p+Au, d+Au and $^3$He+Au collisions at $\sqrt{s_{NN}}$ = 200 GeV as well as in beam energy scan from 19.6 GeV, 39 GeV, 62.4 GeV up to 200GeV d+Au collisions. The results include correlations between forward and backward $\eta$ regions from -3.9 to 3.9 with $\eta$-gap from 0.65 up to 6.2 units. Significant long range and near-side azimuthal correlation have been observed even with the largest $\eta$-gap of 6.2 in central d+Au collisions at 200 GeV and this signal has found to be reduced at lower beam energies. The extracted elliptic event anisotropy parameter $v_2$ and the particle species dependence can be interpreted as a signal of collective and elliptic expansion, which is originally driven by the initially small and anisotropic geometry of fluctuating participant density distribution. A similar but increasing collective and elliptic expansion with system size is observed when comparing the central p+Au, d+Au and $^3$He+Au collisions. The full interpretation of these results would require more complete centrality and rapidity dependent studies yet to come in near future. A new way of extracting the small flow signal from relatively large non-flow background is also discussed in order to observe any small change of the azimuthal correlation function with a reference fitting method which would be suitable especially for multiplicity dependent flow analysis in small system.