A calculation method for higher-order moments of physical quantities, including magnetization and energy, based on the higher-order tensor renormalization group is proposed. The physical observables are represented by impurity tensors. A systematic summation scheme provides coarse-grained tensors including multiple impurities. Our method is compared with the Monte Carlo method on the two-dimensional Potts model. While the nature of the transition of the $q$-state Potts model has been known for a long time owing to the analytical arguments, a clear numerical confirmation has been difficult due to extremely long correlation length in the weakly first-order transitions, e.g., for $q=5$. A jump of the Binder ratio precisely determines the transition temperature. The finite-size scaling analysis provides critical exponents and distinguishes the weakly first-order and the continuous transitions.