Anharmonicity of molecular vibrational motions is closely associated with the thermal property of crystals. However, the origin of anharmonicity is still not fully understood. Low-frequency vibrations, which are usually defined in the terahertz (THz) range, show excellent sensitivity to anharmonicity. In this work, anharmonicity of isonicotinamide in the form I was investigated by using temperature-dependent terahertz time-domain spectroscopy and the quasi-harmonic approximation (QHA) approach at PBE-D3 and PBE-MBD levels. Both DFT calculations suggest the variation of π-π stacking conformation dominates in the thermal expansion of the unit cell. Frequency shifts of the modes in THz range obtained by QHA approach are found to be qualitatively consistent with experimental observations, demonstrating QHA approach is a useful tool for the interpretation of frequency shifts of modes induced by temperature.