Analyzing large power grids directly is computationally expensive. To tackle this issue, various power grid reduction methods, which aim to reduce original power grid to a smaller one while preserving voltage responses at port nodes as much as possible, have been proposed. In the past decade, node elimination based methods incorporated with graph sparsification have shown good performance in the reduction of modern large-scale power grids. The resulting sparsified reduced power grids enable more efficient transient analysis than simulating original power grids directly, with tolerable simulation accuracy. However, due to the large number of ports, the efficiency gained from these elimination-plus-sparsification methods can be still limited. In this work, we propose an effective node aggregation scheme to further reduce the sparsified reduced grids generated by elimination-plus-sparsification methods without increasing simulation errors. We first propose an accuracy-aware multilevel node aggregation framework with an effective scheme to determine aggregation levels automatically, which ensures simulation accuracy does not deteriorate. Then, an aggregation error induced spectral node similarity metric is proposed based on rigorous proof. Experimental results show that the proposed method can reduce the sparsified reduced grids to much smaller sizes with simulation accuracy even improved. It enables a further 3.0X speedup of transient simulation on average. Compared with the transient analysis of original power grids, the proposed method achieves an average 12.4X speedup with average voltage error below 1 mV and maximum error below 6 mV for all test cases with 1.8 V supply voltage. Finally, we demonstrate the advantages of the proposed method over other power grid reduction methods [1], [2] in terms of accuracy, model size and time for reduction.