The function returns the empirical reliability of factor scores. If
\(\sigma^2\) is the sample variance of the estimated scores and \(\bar{sigma}^2_{error}\)
is the average if the squared scores, that is
$$\bar{sigma}^2_{error}=\frac{1}{N}\sum^N_i=1 se_{scores}^2$$
then the subtraction method, a classical reliability estimate similar to
classical test theory is returned using 'sub' yields
$$\frac{\sigma^2-\bar{sigma}^2_{error}}{\sigma^2}$$
for the reliability of the scores. If 'div' is chosen, and alternative
division based approach is used.
$$\frac{\sigma^2}{\sigma^2+\bar{sigma}^2_{error}}$$.
If 'irt' is chosen, a plot returning the standard error of the scores with
the scores is returned per factor.