EFA: N factors are too many for N variables I'm conducting exploratory factor analysis using factanal() function in R.
fit.efa <- factanal(~ x1 + x2 + x3 + x4, 
                factors = 2,  data = df, 
                cor = F,
                na.action = na.omit) 

print(fit.efa, digits=2, cutoff=.3, sort=T)

However, when applying factanal() function, I got the following error:

2 factors are too many for 4 variables

In my case, it is better to use one-factor model for further confirmatory factor analysis, since both Kaiser's criterion and approach based on percentage of explained variance indicated that only one factor is required to represent correctly the latent structure behind my indicators (x1-x4). Nonetheless, according to the chi-squared test, one factor is insufficient (chi = 75.75, df = 2, p = 3.56e-17).
However, I just want to satisfy my interest and to build a two-factor model. So, I don't understand why the factanal() function restricts applying 2 factors for 4 indicators. Every single factor should include at least 2 indicators in order to get an identified model (i.e. df > 0). And in the assumed two-factor model, the number of degrees of freedom is equaled 1. So there are more elements of the input variance-covariance matrix than unknown free parameters, and the unique solution for the model can be found. But I still getting the aforementioned error :(
Could you, please, explain why?
Thanks in advance for your answers and recommendations!

The reason is that factanal implements maximum likelihood estimation, which imposes a constraint.
Given a random sample $X_1,\ldots,X_n$, the factor model assumes that $$X_i-\mu = LF +\Psi,$$
where $L$ is the matrix of factor loadings, $F$ are the factor scores and $\Psi$ is the diagonal matrix of specific variances.
In particular, given an observed iid sample $x_1,\ldots,x_n$ from $N_p(\mu, \Sigma)$, the log-likelihood function is
\begin{eqnarray}
\ell(\theta) & = & \log\prod_{j=1}^n \phi_{p}(x_j;\mu,\Sigma)\\
&= & -\frac{n}{2}\log\vert\Sigma\vert - \sum_{j=1}^n(x_j-\mu)^\top\Sigma^{-1}(x_j-\mu)/2,
\end{eqnarray}
where $\Sigma = LL^\top+\Psi$. It can be shown that the MLE of $\mu$ is $\bar x$. For $L,\Psi$ there is some further work to be done. First, to find a unique solution a constraint has to be placed. This is
$$
L^\top\Psi^{-1} L = \Delta,\quad\text{with}\,\,\,\Delta\,\,\,\text{diagonal matrix}.\quad (*)
$$
The estimation of $L,\Psi$ then proceeds numerically and the solution found has to satisfy $(*)$. It can be shown that the degrees of freedom are
$$\nu = [(p-m)^2 -p-m]/2,$$
where $m$ is the number of factors.
