The noncentrality parameter is $\delta^{2}$, the projection for the restricted model is $P_{r}$, $\beta$ is the vector of true parameters, $X$ is the design matrix for the unrestricted (true) model, $|| x ||$ is the norm:
$$
\delta^{2} = \frac{|| X \beta - P_{r} X \beta ||^{2}}{\sigma^{2}}
$$
You can read the formula like this: $E(y | X) = X \beta$ is the vector of expected values conditional on the design matrix $X$. If you treat $X \beta$ as an empirical data vector $y$, then its projection onto the restricted model subspace is $P_{r} X \beta$, which gives you the prediction $\hat{y}$ from the restricted model for that "data". Consequently, $X \beta - P_{r} X \beta$ is analogous to $y - \hat{y}$ and gives you the error of that prediction. Hence $|| X \beta - P_{r} X \beta ||^{2}$ gives the sum of squares of that error. If the restricted model is true, then $X \beta$ already is within the subspace defined by $X_{r}$, and $P_{r} X \beta = X \beta$, such that the noncentrality parameter is $0$.
You should find this in Mardia, Kent & Bibby. (1980). Multivariate Analysis.
For the first question, the default method in SAS to find the df is not very smart; it looks for terms in the random effect that syntactically include the fixed effect, and uses that. In this case, since trt is not found in ind, it's not doing the right thing. I've never tried BETWITHIN and don't know the details, but either the Satterthwaite option (satterth) or using ind*trt as the random effect give correct results.
PROC MIXED data=Data;
CLASS ind fac trt;
MODEL y = trt /s ddfm=satterth;
RANDOM ind /s;
run;
PROC MIXED data=Data;
CLASS ind fac trt;
MODEL y = trt /s;
RANDOM ind*trt /s;
run;
As for the second question, your SAS code doesn't quite match your R code; it only has a term for fac*ind, while the R code has a term for both ind and fac*ind. (See the Variance Components output to see this.) Adding this gives the same SE for trt in all models in both Q1 and Q2 (0.1892).
As you note, this is an odd model to fit as the fac*ind term has one observation for each level, so is equivalent to the error term. This is reflected in the SAS output, where the fac*ind term has zero variance. This is also what the error message from lme4 is telling you; the reason for the error is that you most likely misspecified something as you're including the error term in the model in two different ways. Interestingly, there is one slight difference in the nlme model; it's somehow finding a variance term for the fac*ind term in addition to the error term, but you will notice that the sum of these two variances equal the error term from both SAS and nlme without the fac*ind term. However, the SE for trt remains the same (0.1892) as trt is nested in ind, so these lower variance terms don't affect it.
Finally, a general note about degrees of freedom in these models: They are computed after the model is fit, and so differences in the degrees of freedom between different programs or options of a program do not necessarily mean that the model is being fit differently. For that, one must look at the estimates of the parameters, both fixed effect parameters and covariance parameters.
Also, using the t and F approximations with a given number of degrees of freedom is fairly controversial. Not only are there several ways to approximate the df, some believe the practice of doing so is not a good idea anyway. A couple words of advice:
If everything is balanced, compare the results with the
traditional least squares method, as they should agree. If it's
close to balanced, compute them yourself (assuming balance) so that
you can make sure the ones you're using are in the right ballpark.
If you have a large sample size, the degrees of freedom don't
matter very much as the distributions get close to a normal and
chi-squared.
Check out Doug Bates's methods for inference. His older method is
based on MCMC simulation; his newer method is based on profiling the
likelihood.
Best Answer
In the Anderson and Rubin Statistic case you compute the test statistic directly from the model without using an estimator: your assumption is that at $\beta=\beta_0$ you have $\mathbb{E}\left( Z^\prime[Y-X\beta]\right)=0$. For each instrument $z_k$ you have a central limit theorem (CLT) $\sqrt{n}(\frac{1}{n}\sum_iz_{k,i}(y_i-x_i^\prime\beta_0))\overset{d}{\to} \mathcal{N}(0,\sigma^2_k)$. Stacking all the moments together you get a vector CLT: $\sqrt{n}(Z^\prime(Y-X\beta_0)/N)\overset{d}{\to} \mathcal{N}(0,V)$ where $V=var(\varepsilon_i Z_i)$. You can transform the normal vector into a $\chi^2$ distribution as follow: $AR(\beta_0)=n(Z^\prime(Y-X\beta_0)/N)^\prime V^{-1} (Z^\prime(Y-X\beta_0)/N) \overset{d}{\to} \chi^2_{K}$. This is the Anderson-Rubin test where $K$ the degrees of freedom is also the number of instruments. The more instruments you have, the larger $K$ and the larger the critical values for the $\chi^2$ distribution will be. Hence it becomes harder to reject when K is large. If your data is homoskedastic, then you can do a susbset test which is dominated by a $\chi^2$ with fewer degrees of freedom but this does not work in the heteroskedastic case if the parameters not in the subset are weakly identified... There is a long litterature on this, I think this summarizes the issue well: http://economics.mit.edu/files/9890