The vibrational properties of solids are fundamental to a large number of physical phenomena, including phase stability and thermal conduction. The canonical approach to modeling these properties requires knowledge of the interatomic forces constants (FCs). The problem of extracting the parameters in the FC expansion from a set of reference forces can be cast in linear form making it amenable to linear regression techniques. Here, we consider the efficiency of various common regression methods for FC extraction and the efficacy of the resulting models for predicting various thermodynamic properties. The regression approach drastically reduces the required number of reference calculations, which constitute the computationally most demanding task in FC extraction, compared to explicit enumeration techniques. It thereby becomes possible to extract both harmonic and high-order anharmonic FCs for large systems with low symmetry, including defects and surfaces. It is shown that ordinary least-squares, especially in connection with feature elimination, often yields the best performance in terms of convergence with respect to training set size and sparsity of the solution. Regression based on the least absolute shrinkage and selection operator (LASSO) on the other hand, while useful in some cases, tends to yield a larger number of features, with a noise level that has a detrimental effect on the prediction of e.g., the thermal conductivity. Finally, we also consider methods for the prediction of the temperature dependence of vibrational spectra from high-order FC expansions via molecular dynamics simulations as well as self-consistent phonons.