Using a diagrammatic approach to Eulerian perturbation theory, we calculate analytically the variance and skewness of the density and velocity divergence induced by gravitational evolution from Gaussian initial conditions, including corrections beyond leading order. Except for the power spectrum, previous calculations in cosmological perturbation theory have been confined to leading order (tree level): we extend these to include loop corrections. For scale-free initial power spectra, P(k) similar to k(n) with -2 less than or equal to n less than or equal to 2, the one-loop variance sigma(2) = [delta(2)] = sigma(l)(2) + 1.82 sigma(l)(4), and the skewness S-3 = [delta(3)]/sigma(4) = 34/7 + 9.8 sigma(l)(2), where sigma(l) is the rms fluctuation of the density field to linear order. (These results depend weakly on the spectral index n, due to the nonlocality of the nonlinear solutions to the equations of motion.) Thus, loop corrections for the (unsmoothed) density held begin to dominate over tree-level contributions (and perturbation theory presumably begins to break down) when sigma(l)(2) similar or equal to 1/2. For the divergence of the velocity field, loop dominance does not occur until sigma(l)(2) approximate to 1. We also compute loop corrections to the variance, skewness, and kurtosis for several nonlinear approximation schemes, where the calculation can be easily generalized to one-point cumulants of higher order and arbitrary number of loops. We find that the Zeldovich approximation gives the best approximation to the loop corrections of exact perturbation theory, followed by the linear potential approximation (LPA) and the frozen flow approximation (FFA), in qualitative agreement with the relative behavior of tree-level results. In LPA and FFA, loop corrections are infrared divergent for spectral indices n less than or equal to -1; this is related to the breaking of Galilean invariance in these schemes.