Using the Lagrange multiplier statistic in regression I don't understand how the LM statistic works because if we run the regression and then record the residual, and then run the regression using the residual as the dependent variable and then using the unrestricted model, how can we run a regression of a constant on 5 different variables? For example, I'm using Stata and here's my commands:
reg bwght cigs parity faminc

where bwght is my dependent variable and cigs, etc. are the dependent variables
scalar epsilon=e(rss)

here is where I saved my error term (I think)
 reg epsilon cigs parity faminc motheduc fatheduc

I tried to run this as a regression but I was told that there was an error because epsilon isn't a variable. Which I can understand why. 
Which leads to my question. How can we conduct a Lagrange Multiplier test?
 A: The LM test is a principle for constructing tests in a variety of situations. In your case, I am assuming that you are interested in the LM test for linear regression specification, in particular for testing for omitted variables in your model. An auxiliary regression of the form you are attempting is a convenient way of computing the LM test statistic.

LM test for omitted variables
Suppose that the linear regression model that you have is composed of two sets of two sets of regressors -- those collected in the vector $\boldsymbol{X}_{1i}$ and those in the vector $\boldsymbol{X}_{2i}$. 
$$
Y_i = \boldsymbol{X}_{1i}'\boldsymbol{\beta}_1 + \boldsymbol{X}_{2i}'\boldsymbol{\beta}_2+\varepsilon_i
$$
Suppose you are interested in testing whether the regressors $\boldsymbol{X}_{2i}$ belong to the model or not (that $\boldsymbol{\beta}_2=\boldsymbol{0}$). The LM test in this context involves estimating the restricted (under the hypothesis) model
$$
Y_i = \boldsymbol{X}_{1i}'\widetilde{\boldsymbol{\beta}}_1 +\widetilde{\varepsilon}_i
$$
where the $\widetilde{\cdot}$ indicate estimates. You then use the constructed residuals, $\widetilde{\varepsilon}_i$ in an auxiliary regression on the full set of regressors, 
$$
\widetilde{\varepsilon}_i =\boldsymbol{X}_{1i}'\boldsymbol{\gamma}_1 + \boldsymbol{X}_{2i}'\boldsymbol{\gamma}_2+\nu_i
$$ 
Estimate this regression, and compute the statistic $nR^2$, where $n$ is the sample size, and $R^2$ is the uncenterd coefficient of determination. Under the null, this is distributed $\chi^2(k_2)$, where $k_2$ is the number of regressors omitted from the main model. The decision rule is $\mathbf{1}_{[nR^2 > c_{1-\alpha}]}$, where $c_{1-\alpha}$ is the $1-\alpha$-th quantile from the $\chi^2(k_2)$ distribution.

Computing the LM test in Stata
Here is an example taken from the Wooldridge book (Example 5.3).
use "http://fmwww.bc.edu/ec-p/data/wooldridge/crime1", clear
reg narr86 pcnv ptime86 qemp86  // restricted regression
predict epstilde, resid  // predict the residuals
reg epstilde pcnv ptime86 qemp86 avgsen tottime  // auxiliary regression
predict nuhat, resid

mat accum mResid = epstilde nuhat, noconstant
di "The LM test statistic is: " e(N)*(1 - mResid[2,2]/mResid[1,1]) ///
    " and the 10% critical value is: " invchi2tail(2, 0.1)

which gives the output
. di "The LM test statistic is: " e(N)*(1 - mResid[2,2]/mResid[1,1]) ///
>     " and the 10% critical value is: " invchi2tail(2, 0.1)
The LM test statistic is: 4.0707559 and the 10% critical value is: 4.6051702

This demonstrates that you are unable to reject the null hypothesis that the omitted regressors do not belong to the model.
A: The Wooldridge example from Fg Nu can be improved upon in a couple of ways. First, to get the exact p value for test statistic, we can change the final line to:
scalar LM = e(N)*(1 - mResid[2,2]/mResid[1,1])
di "The LM test statistic is: " LM " and the associated p value is: " chi2tail(2, LM)

Which gives the output:
The LM test statistic is: 4.0707072 and the associated p value is: .13063428

Also, we can use Stata's test command to simplify the code, at the pedagogic cost of hiding some math. And we can use the constraint facility to demonstrate how to generalize to hypotheses that are not simply that certain coefficients are zero:
use "http://fmwww.bc.edu/ec-p/data/wooldridge/crime1", clear
constraint 1 avgsen = 0
constraint 2 tottime= 0
cnsreg narr86 pcnv ptime86 qemp86 avgsen tottime, c(1 2)  // restricted regression
predict double epstilde, resid  // predict the residuals
reg epstilde pcnv ptime86 qemp86 avgsen tottime  // auxiliary regression
test (avgsen = 0) (tottime = 0)

Which gives the output:
 ( 1)  avgsen = 0
 ( 2)  tottime = 0

       F(  2,  2719) =    2.03
            Prob > F =    0.1310

The result looks superficially different because the regress and test commands incorporate small-sample degrees-of-freedom corrections. So the previous chi2 test statistic is divided by 2 (since there are 2 constraints) and viewed as an F statistic. The resulting p value is nearly the same.
