Neutral current quasielastic (anti)neutrino scattering cross sections on a $^{12}$C target are analyzed using a realistic spectral function $S(p,E)$ that gives a scaling function in accordance with the ($e,e$) scattering data. The spectral function accounts for the nucleon-nucleon (NN) correlations by using natural orbitals (NOs) from the Jastrow correlation method and has a realistic energy dependence. The standard value of the axial mass $M_A= 1.032$ GeV is used in all calculations. The role of the final-state interaction (FSI) on the spectral and scaling functions, as well as on the cross sections is accounted for. A comparison of the calculations with the empirical data of the MiniBooNE and BNL experiments is performed. Our results are analyzed in comparison with those when NN correlations are not included, and also with results from other theoretical approaches, such as the relativistic Fermi gas (RFG), the relativistic mean field (RMF), the relativistic Greens function (RGF), as well as with the SuperScaling Approach (SuSA) based on the analysis of quasielastic electron scattering.