Isostasy is a key concept in geodesy and geophysics. The classical isostatic models of Airy/Heiskanen and Pratt/Hayford imply that the topographic mass surplus and ocean mass deficit are balanced by mountain roots and anti-roots in the former model and by density variations in the topography and the compensation layer below sea bottom in the latter model. In geophysics gravity inversion is an essential topic where isostasy comes to play. The main objective of this study is to compare the prediction of geoid heights from the above isostatic models based on matched asymptotic expansion with geoid heights observed by the Earth Gravitational Model 2008. Numerical computations were carried out both globally and in several regions, showing poor agreements between the theoretical and observed geoid heights. As an alternative, multiple regression analysis including several non-isostatic terms in addition to the isostatic terms was tested providing only slightly better success rates. Our main conclusion is that the geoid height cannot generally be represented by the simple formulas based on matched asymptotic expansions. This is because (a) both the geoid and isostatic compensation of the topography have regional to global contributions in addition to the pure local signal considered in the classical isostatic models, and (b) geodynamic phenomena are still likely to significantly blur the results despite that all spherical harmonic low-degree (below degree 11) gravity signals were excluded from the study.