Analytic methods to quantify power systems stability margins while acknowledging uncertainty will be critical in dynamic security assessment due to the integration of renewable resources and uncontrollable loads. Lyapunov's direct method provides a mathematical framework to assess stability margins. However, analytical construction of Lyapunov functions are limited to very restrictive settings and cannot be extended to accommodate the impact of uncertainty. This paper focuses on the algorithmic construction of Lyapunov functions and the estimation of the robust Region-Of-Attraction (ROA) with sum-of-squares (SOS) optimization programs which can be translated into semidefinite problems and then solved with readily available software. Key to the proposed approach is the formulation of computationally tractable polynomial differential equation models that are affine in parametric/input uncertainty. Numerical case studies for a prototypical power system model validate the accuracy of the polynomial approximations and compare the ROA and the robust ROA with results from repeated time-domain simulations.