Option to compute corrected error estimate

[marklhc]

For a unidimensional test with factor score scoring vector $\boldsymbol{a}^\top$, when the factor scoring weights are determined by sample estimates $\hat{\boldsymbol{a}}$ with asymptotic covariance matrix $V{\hat{\boldsymbol{a}}}$, the expected conditional variance of the factor scores is (using law of total variance): $$\mathrm{Var}(\hat{\boldsymbol{a}}^\top \boldsymbol{Y} \mid \theta) = \mathrm{E}(\hat{\boldsymbol{a}})^\top \boldsymbol{\Theta} \mathrm{E}(\hat{\boldsymbol{a}}) + \mathrm{Tr}(V{\hat{\boldsymbol{a}}} \boldsymbol{\Theta})$$

The first component is approximated by substituting $\mathrm{E}(\hat{\boldsymbol{a}})$ by a consistent estimate, as is currently done. The second component is additional error due to sampling error in $\hat{\boldsymbol{a}}$, and will be close to zero in large sample. However, in small samples, it is worth adding the correction factor.

We can use the delta method to approximate $V_{\hat{\boldsymbol{a}}}$, with the gradients computed using finite difference (as in grand standardization) with the internal lavaan functions.

This can be extended to multidimensional tests, and probably require vectorizing the scoring matrix.

[marklhc]

Need to accommodate multiple groups

[marklhc]

Need to add a unit test.