High entropy alloys (HEAs) are a series of novel materials that demonstrate many exceptional mechanical properties. To understand the origin of these attractive properties, it is important to investigate the thermodynamics and elucidate the evolution of various chemical phases. In this work, we introduce a data-driven approach to construct the effective Hamiltonian and study the thermodynamics of HEAs through canonical Monte Carlo simulation. The main characteristic of our method is to use pairwise interactions between atoms as features and systematically improve the representativeness of the dataset using samples from Monte Carlo simulation. We find this method produces highly robust and accurate effective Hamiltonians that give less than 0.1 mRy test error for all the three refractory HEAs: MoNbTaW, MoNbTaVW, and MoNbTaTiW. Using replica exchange to speed up the MC simulation, we calculated the specific heats and short-range order parameters in a wide range of temperatures. For all the studied materials, we find there are two major order-disorder transitions occurring respectively at $T_1$ and $T_2$, where $T_1$ is near room temperature but $T_2$ is much higher. We further demonstrate that the transition at $T_1$ is caused by W and Nb while the one at $T_2$ is caused by the other elements. By comparing with experiments, {color{black} the results provide insight into the role of chemical ordering in the strength and ductility of HEAs.