In a previous paper (I) we carried out a model calculation for nematic liquid crystals which treats spatial correlations between molecules as for a classical isotropic liquid and orientational correlations in the meanfield approximation. Generalized Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) equations were derived for the density function and an orientationally averaged pair-correlation function. They were solved for a simple pairwise potential chosen to fit certain experimental data on $p$-azoxyanisole (PAA). A free-energy expression was obtained by identifying its Euler-Lagrange equation to the first BBGKY equation. In this paper we reformulated the molecular theory of nematic liquid crystals rigorously on the basis of a general probability distribution function ${P}_{N}(1,\dots{},N)$. It was shown that both the mean-field theory (MFT) and the above-stated orientation-averaged pair-correlation theory (OAPC) arise from applying low-order approximations to ${P}_{N}$. In the case of MFT, ${P}_{N}$ was approximated as a product of single-particle functions. In the case of OAPC, spatial correlations were introduced, but decoupling of spatial and orientational dependences were assumed. It is now clear how higher-order-approximation schemes can be systematically constructed. Furthermore, for each scheme chosen, the present formulation offers a uniquely consistent formula for the free energy. This places our previous work on a firm theoretical foundation. The theory in OAPC approximation was applied to model systems described by a more general pairwise potential than that employed in I. Potential parameters were chosen to mimic PAA and methoxybenzylidene butylaniline. For each substance we calculated the isotropic-nematic coexistence curve, temperature dependence of the order parameter, latent heat, and volume change at transition. The results were compared to experimental data.