There are two main contributions that can be included in the DFT context when calculating phonons at finite temperature:
- Thermal expansion. As temperature increases, the volume of materials changes (typically increases). If you have the volume as a function of temperature, then you can calculate the phonons at each of these volumes to obtain the thermal expansion contribution to finite temperature phonons. To calculate the volumes at each temperature, a commonly used approach is the quasi-harmonic approximation. In this method, you minimize the free energy at a given temperature with respect to the volume of the system, where the free energy is calculated using the harmonic approximation at each volume. To do these calculations in VASP, you simply need to do a series of harmonic phonon calculations at different volumes.
- Phonon-phonon interactions. As temperature increases, the harmonic approximation (in which phonons are non-interacting) becomes increasingly worse, and anharmonic phonon-phonon interactions become important. These interactions change with temperature, and lead to a renormalization of phonon frequencies with temperature. There are multiple ways in which these corrections can be incorporated, for example using the self-consistent harmonic approximation, which, roughly speaking, finds the "best" harmonic fit to the anharmonic energy. I do not think that any of these methods are natively available in VASP, and instead you will probably have to use another code (potentially interfaced with VASP), such as SCAILD or SSCHA.
Overall, both effects can be important. I am not familiar with the structures you are interested in, but in the oxide perovskties (e.g. BaTiO$_3$) which also exhibit temperature driven phase transitions between different structures, phonon-phonon interactions are the dominant mechanism.