Categories &

Functions List

Function Reference: factoran

statistics: lambda = factoran (X, m)
statistics: [lambda, psi] = factoran (X, m)
statistics: [lambda, psi, T] = factoran (X, m)
statistics: [lambda, psi, T, stats] = factoran (X, m)
statistics: [lambda, psi, T, stats, F] = factoran (X, m)
statistics: […] = factoran (…, Name, Value)

Common factor analysis.

lambda = factoran (X, m) fits the common factor model with m common factors to the N-by-P data matrix X, whose rows are observations and columns are variables, and returns the P-by-m matrix lambda of factor loadings. The model is

 
 x = mu + lambda * f + e

where f are the common factors and e the variable-specific errors, uncorrelated with each other and with f. The analysis is carried out on the correlation matrix, so the loadings are in standardized units and a variable’s communality sum (lambda(i,:) .^ 2) and its specific variance psi(i) sum to one.

[lambda, psi] = factoran (…) also returns the P-by-1 vector psi of specific variances.

[lambda, psi, T] = factoran (…) also returns the m-by-m rotation matrix T that was applied to the loadings. It is the identity when 'Rotate' is 'none'.

[lambda, psi, T, stats] = factoran (…) also returns a structure stats with the fields

FieldDescription
loglikethe maximized log-likelihood, up to a constant.
dfethe error degrees of freedom, ((p - m)^2 - p - m) / 2.
chisqthe likelihood ratio statistic testing m common factors against an unrestricted covariance.
pthe significance of chisq.

Source Code: factoran

The last two are present only when they can be computed: they are omitted when the degrees of freedom are not positive, when a specific variance has reached its lower bound (a Heywood case, where the likelihood is on the boundary), and for the 'paf' extraction, which does not maximize a likelihood. Test for them with isfield before using them.

[lambda, psi, T, stats, F] = factoran (…) also returns the N-by-m matrix F of predicted factor scores. Scores are not available from a covariance matrix, only from data.

Name-Value pairs

NameValue
'Extraction'How the loadings are estimated, either 'ml' (default) for maximum likelihood or 'paf' for principal axis factoring. An Octave extension; MATLAB fits by maximum likelihood only. See the note below.
'Xtype'Whether X holds 'data' (default) or a 'covariance' (or correlation) matrix.
'Nobs'The number of observations behind a covariance matrix. Required for stats when 'Xtype' is 'covariance'.
'Delta'The lower bound on the specific variances, a scalar in [0, 1) (default 0.005). Bounding them away from zero keeps the likelihood finite; a solution that reaches the bound is a Heywood case and is reported by a warning.
'Rotate'The rotation applied to the loadings, passed to rotatefactors: 'varimax' (default), 'none', 'quartimax', 'equamax', 'parsimax', 'orthomax', or 'promax'.
'Normalize'Whether the rotation normalizes the rows of the loadings (Kaiser normalization), 'on' (default) or 'off'.
'Power'The exponent of the 'promax' target, a scalar not less than 1 (default 4).
'Scores'How F is predicted, either 'wls' (default, also named 'Bartlett') or 'regression' (also named 'Thomson').
'Maxit'The iteration limit of the extraction (default 500).
'Tolerance'The convergence tolerance of the extraction (default 1e-8).

Source Code: factoran

Choosing the extraction

The two extractions fit the same model but estimate it differently, and they answer to different circumstances.

'ml' maximizes the likelihood of a multivariate normal, and is what MATLAB’s factoran does. Use it when you want the likelihood ratio test in stats to decide how many factors the data support, and when the data are plausibly normal. It can fail to converge, or push a specific variance to zero, when the model asks for more factors than the data hold.

'paf' iterates communalities on the reduced correlation matrix. It makes no distributional assumption and is stable where maximum likelihood struggles, which is why it remains available here, but it provides no likelihood and therefore no test: stats carries only dfe.

The two agree closely when the model fits the data well and diverge when it does not, so a large difference between them is itself informative.

Number of factors

m must leave the model identified, that is (p - m)^2 >= p + m. With six variables at most three factors can be fitted, and only the smaller counts leave degrees of freedom to test.

References

  1. Lawley, D. N., and Maxwell, A. E., Factor Analysis as a Statistical Method, 2nd Edition, Butterworths, 1971.
  2. Joreskog, K. G., "Some contributions to maximum likelihood factor analysis", Psychometrika 32(4), 443-482, 1967.
  3. Harman, H. H., Modern Factor Analysis, 3rd Edition, University of Chicago Press, 1976.

See also: rotatefactors, pca, pcacov, princomp, barttest

Source Code: factoran

Six measured variables built from two underlying factors, plus noise. Factor analysis recovers the structure without being told it: the first three variables load on one factor and the last three on the other.

 rng (42);
 F = randn (300, 2);
 X = F * [0.8 0.1; 0.7 0.2; 0.75 0.15; 0.15 0.8; 0.2 0.7; 0.1 0.75]' ...
     + 0.6 * randn (300, 6);
 lambda = factoran (X, 2);
 printf ("loadings on the two rotated factors:\n");
loadings on the two rotated factors:
 disp (round (lambda * 1000) / 1000);
   0.863000   0.077000
   0.699000   0.264000
   0.749000   0.148000
   0.189000   0.798000
   0.290000   0.701000
   0.040000   0.819000

How many factors do the data support? The likelihood ratio test in stats answers it. These data were built from two factors, and the test rejects one factor while accepting two. Three factors leave no degrees of freedom, so there is nothing left to test with.

 rng (42);
 F = randn (300, 2);
 X = F * [0.8 0.1; 0.7 0.2; 0.75 0.15; 0.15 0.8; 0.2 0.7; 0.1 0.75]' ...
     + 0.6 * randn (300, 6);
 for m = 1:3
   [~, ~, ~, stats] = factoran (X, m);
   if (isfield (stats, "p"))
     printf ("%d factor(s): chisq = %8.3f, dfe = %d, p = %.4f\n", ...
             m, stats.chisq, stats.dfe, stats.p);
   else
     printf ("%d factor(s): nothing to test against (dfe = %d)\n", ...
             m, stats.dfe);
   endif
 endfor
1 factor(s): chisq =  287.486, dfe = 9, p = 0.0000
2 factor(s): chisq =    3.645, dfe = 4, p = 0.4561
warning: factoran: some specific variances are at their lower bound; the fit is a Heywood case.
warning: called from
    factoran at line 309 column 5
    __eval_demo__ at line 84 column 9
    __demo_notebook__ at line 43 column 3
    __build_demos__ at line 79 column 7
    function_texi2html at line 135 column 5
    package_texi2html at line 336 column 9

3 factor(s): nothing to test against (dfe = 0)

Rotation decides how a fit is presented, not how good it is. The unrotated solution puts most of the variance on a general first factor; varimax turns it so each variable loads mainly on one factor, which is easier to read. The specific variances are untouched either way.

 rng (42);
 F = randn (300, 2);
 X = F * [0.8 0.1; 0.7 0.2; 0.75 0.15; 0.15 0.8; 0.2 0.7; 0.1 0.75]' ...
     + 0.6 * randn (300, 6);
 [Lnone, psi_none] = factoran (X, 2, "Rotate", "none");
 [Lvari, psi_vari, T] = factoran (X, 2, "Rotate", "varimax");
 printf ("unrotated:\n");  disp (round (Lnone * 100) / 100);
unrotated:
   0.6900  -0.5300
   0.6900  -0.2800
   0.6500  -0.4000
   0.6800   0.4600
   0.6900   0.3200
   0.5800   0.5800
 printf ("varimax:\n");    disp (round (Lvari * 100) / 100);
varimax:
   0.860000   0.080000
   0.700000   0.260000
   0.750000   0.150000
   0.190000   0.800000
   0.290000   0.700000
   0.040000   0.820000
 printf ("the rotation matrix takes one to the other: %d\n", ...
         max (max (abs (Lnone * T - Lvari))) < 1e-10);
the rotation matrix takes one to the other: 1
 printf ("specific variances unchanged: %d\n", ...
         max (abs (psi_none - psi_vari)) < 1e-10);
specific variances unchanged: 1

Factor scores place each observation on the factors, so they can be plotted or used as inputs downstream. The two predictors optimise different things and are not equal, but they agree closely on the ordering of observations.

 rng (42);
 F = randn (300, 2);
 X = F * [0.8 0.1; 0.7 0.2; 0.75 0.15; 0.15 0.8; 0.2 0.7; 0.1 0.75]' ...
     + 0.6 * randn (300, 6);
 [~, ~, ~, ~, Fwls] = factoran (X, 2, "Scores", "wls");
 [~, ~, ~, ~, Freg] = factoran (X, 2, "Scores", "regression");
 printf ("first three observations, weighted least squares:\n");
first three observations, weighted least squares:
 disp (round (Fwls(1:3,:) * 1000) / 1000);
  -0.603000  -0.996000
  -0.319000   0.097000
  -1.115000   1.357000
 printf ("first three observations, regression:\n");
first three observations, regression:
 disp (round (Freg(1:3,:) * 1000) / 1000);
  -0.555000  -0.855000
  -0.262000   0.064000
  -0.864000   1.066000
 printf ("the two agree on factor 1 to a correlation of %.4f\n", ...
         corr (Fwls(:,1), Freg(:,1)));
the two agree on factor 1 to a correlation of 0.9981

Two ways to estimate the same model. Maximum likelihood is the default and is what MATLAB does; principal axis factoring is an Octave extension that assumes no distribution. They agree closely when the model fits, so a large gap between them is a warning about the fit. Only maximum likelihood carries a test.

 rng (42);
 F = randn (300, 2);
 X = F * [0.8 0.1; 0.7 0.2; 0.75 0.15; 0.15 0.8; 0.2 0.7; 0.1 0.75]' ...
     + 0.6 * randn (300, 6);
 Lml = factoran (X, 2, "Extraction", "ml");
 Lpaf = factoran (X, 2, "Extraction", "paf");
 printf ("largest loading difference between the extractions: %.4f\n", ...
         max (abs (abs (Lml(:)) - abs (Lpaf(:)))));
largest loading difference between the extractions: 0.0047
 [~, ~, ~, sml] = factoran (X, 2, "Extraction", "ml");
 [~, ~, ~, spaf] = factoran (X, 2, "Extraction", "paf");
 printf ("ml  reports: %s\n", strjoin (fieldnames (sml)', ", "));
ml  reports: loglike, dfe, chisq, p
 printf ("paf reports: %s\n", strjoin (fieldnames (spaf)', ", "));
paf reports: loglike, dfe