In a previous work we developed a convex infinite dimensional linear programming (LP) approach to approximating the region of attraction (ROA) of polynomial dynamical systems subject to compact basic semialgebraic state constraints. Finite dimensional relaxations to the infinite-dimensional LP lead to a truncated moment problem in the primal and a polynomial sum-of-squares problem in the dual. This primal-dual linear matrix inequality (LMI) problem can be solved numerically with standard semidefinite programming solvers, producing a hierarchy of outer (i.e. exterior) approximations of the ROA by polynomial sublevel sets, with a guarantee of almost uniform and set-wise convergence. In this companion paper, we show that our approach is flexible enough to be modified so as to generate a hierarchy of polynomial inner (i.e.,interior) approximations of the ROA with similar convergence guarantees.