The temperature dependencies of NMR shifts in the critical region of two coexisting phases have been simulated using statistical thermodynamics and graph-theory consideration of equilibrium processes of molecular association. Microparameters of magnetic screening of various water and water/pyridine structures used in the statistical averaging have been evaluated by density functional theory calculations (PBE1PBE and B3PW91 functionals in the 6-311++G** basis set). The gauge-including atomic orbital (GIAO) approach has been applied to ensure gauge invariance of the results. Solvent effects were taken into account by a polarized continuum model (PCM). NMR shifts "order parameters" (Deltadelta = |delta+ - delta-|) and "diameters" (phidelta = |(delta+ + delta-)/2 - deltaC|, where delta+, delta-, and deltaC are the chemical shifts of coexisting phases and at the critical point respectively) have been calculated in each case close to the lower critical solution point (TL) and processed using linear regression analysis of Deltadelta approximately |T - TL| and phidelta approximately |T - TL| in the log-log plot. It has been shown that the critical index beta can be evaluated with high precision from the slope of Deltadelta = f(T - TL) at any realistic set of model input parameters. The slope of diameter has been found to depend on both input beta and alpha values. The obtained phidelta slopes (0.58-0.63) are very close to 2beta values. The results are discussed within the concept of complete scaling. Results of simulation are compared and supported by experimental NMR data for water/2,6-lutidine, acetic anhydride/n-heptane, and acetic anhydride/cyclohexane systems.