Anharmonic quantum chemical calculations were employed to simulate and interpret a near-infrared (NIR) spectrum of caffeine. First and second overtones, as well as binary and ternary combination bands, were obtained, accurately reproducing the lineshape of the experimental spectrum in the region of 10,000-4000 cm-1 (1000-2500 nm). The calculations enabled performing a detailed analysis of NIR spectra of caffeine, including weak bands due to the second overtones and ternary combinations. A highly convoluted nature of NIR spectrum of caffeine was unveiled, with numerous overlapping bands found beneath the observed spectral lineshape. To properly reflect that intrinsic complexity, the band assignments were provided in the form of heat maps presenting the contributions to the NIR spectrum from various kinds of vibrational transitions. These contributions were also quantitatively assessed in terms of the integral intensities. It was found that the combination bands provide the decisively dominant contributions to the NIR spectrum of caffeine. The first overtones gain significant importance between 6500-5500 cm-1, while the second overtones are meaningful in the higher wavenumber regions, particularly in the 10,000-7000 cm-1 region. The obtained detailed band assignments enabled deep interpretation of the absorption regions of caffeine identified in the literature as meaningful for analytical applications of NIR spectroscopy focused on quantitative analysis of caffeine content in drugs and natural products.