Three hundred and nine carbon-carbon, carbon-nitrogen, and carbon-oxygen pi-bond lengths in high precision crystal structures of 31 purine and pyrimidine nucleobases were related to the Pauling pi-bond order, its analogues corrected to crystal packing effects, the numbers of non-hydrogen atoms around the bond, and the sum of atomic numbers of the bond atoms. Principal Component Analysis (PCA) and Hierarchical Cluster Analysis (HCA) demonstrated that the bond lengths in the nucleobases are three-dimensional phenomenon, characterized by nine distinct classes of bonds. Bond lengths predicted by Linear Regression models, Pauling Harmonic Potential Curves, Multiple Linear Regression, Principal Component, and Partial Least Squares Regression were compared to those calculated by molecular mechanics, semiempirical, and ab initio methods using PCA-HCA procedure on the calculated bond lengths, statistical parameters, and structural aromaticity indices. Incorporation of crystal packing effects into bond orders makes multivariate models to be competitive to semiempirical results, while further improvement of quantum chemical calculations can be achieved by geometry optimization of molecular clusters.43787-809