Using the next-to-leading order QCD-corrected effective Hamiltonian, charmless exclusive nonleptonic decays of the B s meson into η or η ′ are calculated within the generalized factorization approach. Nonfactorizable contributions, which can be parametrized in terms of the effective number of colors N eff c for P P and V P decay modes, are studied in two different schemes: (i) the one with the "homogeneous" structure in which (is assumed, and (ii) the "heterogeneous" one in which the possibility ofis considered, where N eff c (V + A) denotes the effective value of colors for the (V − A)(V + A) penguin operators and N eff c (V − A) for the (V − A)(V − A) ones. For processes depending on the N eff c -stable a i such as B s → (π, ρ)(η ′ , η), the predicted branching ratios are not sensitive to the factorization approach we choose. While for the processes depending on the N eff c -sensitive a i such asB s → ωη ( ′ ) , there is a wide range for the branching ratios depending on the choice of the N eff c involved. We have included the QCD anomaly effect in our calculations and found that it is important forB s → η ( ′ ) η ( ′ ) . The effect of the (cc) → η ′ mechanism is found to be tiny due to a possible CKM-suppression and the suppression in the decay constants except for theB s → φη decay within the "naive" factorization approach, where the internal W diagram is CKM-suppressed and the penguin contributions are compensated.