G03CAF (PDF version)
G03 Chapter Contents
G03 Chapter Introduction
NAG Library Manual

NAG Library Routine Document

G03CAF

Note:  before using this routine, please read the Users' Note for your implementation to check the interpretation of bold italicised terms and other implementation-dependent details.

+ Contents

    1  Purpose
    7  Accuracy

1  Purpose

G03CAF computes the maximum likelihood estimates of the parameters of a factor analysis model. Either the data matrix or a correlation/covariance matrix may be input. Factor loadings, communalities and residual correlations are returned.

2  Specification

SUBROUTINE G03CAF ( MATRIX, WEIGHT, N, M, X, LDX, NVAR, ISX, NFAC, WT, E, STAT, COM, PSI, RES, FL, LDFL, IOP, IWK, WK, LWK, IFAIL)
INTEGER  N, M, LDX, NVAR, ISX(M), NFAC, LDFL, IOP(5), IWK(4*NVAR+2), LWK, IFAIL
REAL (KIND=nag_wp)  X(LDX,M), WT(*), E(NVAR), STAT(4), COM(NVAR), PSI(NVAR), RES(NVAR*(NVAR-1)/2), FL(LDFL,NFAC), WK(LWK)
CHARACTER(1)  MATRIX, WEIGHT

3  Description

Let p variables, x1,x2,,xp, with variance-covariance matrix Σ be observed. The aim of factor analysis is to account for the covariances in these p variables in terms of a smaller number, k, of hypothetical variables, or factors, f1,f2,,fk. These are assumed to be independent and to have unit variance. The relationship between the observed variables and the factors is given by the model:
xi=j=1kλijfj+ei,  i=1,2,,p
where λij, for i=1,2,,p and j=1,2,,k, are the factor loadings and ei, for i=1,2,,p, are independent random variables with variances ψi, for i=1,2,,p. The ψi represent the unique component of the variation of each observed variable. The proportion of variation for each variable accounted for by the factors is known as the communality. For this routine it is assumed that both the k factors and the ei's follow independent Normal distributions.
The model for the variance-covariance matrix, Σ, can be written as:
Σ=ΛΛT+Ψ (1)
where Λ is the matrix of the factor loadings, λij, and Ψ is a diagonal matrix of unique variances, ψi, for i=1,2,,p.
The estimation of the parameters of the model, Λ and Ψ, by maximum likelihood is described by Lawley and Maxwell (1971). The log-likelihood is:
-12n-1logΣ-12n-1traceS,Σ-1+constant,
where n is the number of observations, S is the sample variance-covariance matrix or, if weights are used, S is the weighted sample variance-covariance matrix and n is the effective number of observations, that is, the sum of the weights. The constant is independent of the parameters of the model. A two stage maximization is employed. It makes use of the function FΨ, which is, up to a constant, -2/n-1 times the log-likelihood maximized over Λ. This is then minimized with respect to Ψ to give the estimates, Ψ^, of Ψ. The function FΨ can be written as:
FΨ=j=k+1pθj-logθj-p-k
where values θj, for j=1,2,,p are the eigenvalues of the matrix:
S*=Ψ-1/2SΨ-1/2.
The estimates Λ^, of Λ, are then given by scaling the eigenvectors of S*, which are denoted by V:
Λ^=Ψ1/2VΘ-I1/2.
where Θ is the diagonal matrix with elements θi, and I is the identity matrix.
The minimization of FΨ is performed using E04LBF which uses a modified Newton algorithm. The computation of the Hessian matrix is described by Clark (1970). However, instead of using the eigenvalue decomposition of the matrix S* as described above, the singular value decomposition of the matrix RΨ-1/2 is used, where R is obtained either from the QR decomposition of the (scaled) mean centred data matrix or from the Cholesky decomposition of the correlation/covariance matrix. The routine E04LBF ensures that the values of ψi are greater than a given small positive quantity, δ, so that the communality is always less than one. This avoids the so called Heywood cases.
In addition to the values of Λ, Ψ and the communalities, G03CAF returns the residual correlations, i.e., the off-diagonal elements of C-ΛΛT+Ψ where C is the sample correlation matrix. G03CAF also returns the test statistic:
χ2=n-1-2p+5/6-2k/3FΨ^
which can be used to test the goodness-of-fit of the model (1), see Lawley and Maxwell (1971) and Morrison (1967).

4  References

Clark M R B (1970) A rapidly convergent method for maximum likelihood factor analysis British J. Math. Statist. Psych.
Hammarling S (1985) The singular value decomposition in multivariate statistics SIGNUM Newsl. 20(3) 2–25
Lawley D N and Maxwell A E (1971) Factor Analysis as a Statistical Method (2nd Edition) Butterworths
Morrison D F (1967) Multivariate Statistical Methods McGraw–Hill

5  Parameters

1:     MATRIX – CHARACTER(1)Input
On entry: selects the type of matrix on which factor analysis is to be performed.
MATRIX='D'
The data matrix will be input in X and factor analysis will be computed for the correlation matrix.
MATRIX='S'
The data matrix will be input in X and factor analysis will be computed for the covariance matrix, i.e., the results are scaled as described in Section 8.
MATRIX='C'
The correlation/variance-covariance matrix will be input in X and factor analysis computed for this matrix.
Constraint: MATRIX='D', 'S' or 'C'.
2:     WEIGHT – CHARACTER(1)Input
On entry: if MATRIX='D' or 'S', WEIGHT indicates if weights are to be used.
WEIGHT='U'
No weights are used.
WEIGHT='W'
Weights are used and must be supplied in WT.
Note:  if MATRIX='C', WEIGHT is not referenced.
Constraint: if MATRIX='D' or 'S', WEIGHT='U' or 'W'.
3:     N – INTEGERInput
On entry: if MATRIX='D' or 'S' the number of observations in the data array X.
If MATRIX='C' the (effective) number of observations used in computing the (possibly weighted) correlation/variance-covariance matrix input in X.
Constraint: N>NVAR.
4:     M – INTEGERInput
On entry: the number of variables in the data/correlation/variance-covariance matrix.
Constraint: MNVAR.
5:     X(LDX,M) – REAL (KIND=nag_wp) arrayInput
On entry: the input matrix.
If MATRIX='D' or 'S', X must contain the data matrix, i.e., Xij must contain the ith observation for the jth variable, for i=1,2,,n and j=1,2,,M.
If MATRIX='C', X must contain the correlation or variance-covariance matrix. Only the upper triangular part is required.
6:     LDX – INTEGERInput
On entry: the first dimension of the array X as declared in the (sub)program from which G03CAF is called.
Constraints:
  • if MATRIX='D' or 'S', LDXN;
  • if MATRIX='C', LDXM.
7:     NVAR – INTEGERInput
On entry: p, the number of variables in the factor analysis.
Constraint: NVAR2.
8:     ISX(M) – INTEGER arrayInput
On entry: ISXj indicates whether or not the jth variable is included in the factor analysis. If ISXj1, the variable represented by the jth column of X is included in the analysis; otherwise it is excluded, for j=1,2,,M.
Constraint: ISXj>0 for NVAR values of j.
9:     NFAC – INTEGERInput
On entry: k, the number of factors.
Constraint: 1NFACNVAR.
10:   WT(*) – REAL (KIND=nag_wp) arrayInput
Note: the dimension of the array WT must be at least N if WEIGHT='W' and MATRIX='D' or 'S', and at least 1 otherwise.
On entry: if WEIGHT='W' and MATRIX='D' or 'S', WT must contain the weights to be used in the factor analysis. The effective number of observations in the analysis will then be the sum of weights. If WTi=0.0, the ith observation is not included in the analysis.
If WEIGHT='U' or MATRIX='C', WT is not referenced and the effective number of observations is n.
Constraint: if WEIGHT='W', the sum of weights >NVAR, WTi0.0, for i=1,2,,n.
11:   E(NVAR) – REAL (KIND=nag_wp) arrayOutput
On exit: the eigenvalues θi, for i=1,2,,p.
12:   STAT(4) – REAL (KIND=nag_wp) arrayOutput
On exit: the test statistics.
STAT1
Contains the value FΨ^.
STAT2
Contains the test statistic, χ2.
STAT3
Contains the degrees of freedom associated with the test statistic.
STAT4
Contains the significance level.
13:   COM(NVAR) – REAL (KIND=nag_wp) arrayOutput
On exit: the communalities.
14:   PSI(NVAR) – REAL (KIND=nag_wp) arrayOutput
On exit: the estimates of ψi, for i=1,2,,p.
15:   RES(NVAR×NVAR-1/2) – REAL (KIND=nag_wp) arrayOutput
On exit: the residual correlations. The residual correlation for the ith and jth variables is stored in RESj-1j-2/2+i, i<j.
16:   FL(LDFL,NFAC) – REAL (KIND=nag_wp) arrayOutput
On exit: the factor loadings. FLij contains λij, for i=1,2,,p and j=1,2,,k.
17:   LDFL – INTEGERInput
On entry: the first dimension of the array FL as declared in the (sub)program from which G03CAF is called.
Constraint: LDFLNVAR.
18:   IOP(5) – INTEGER arrayInput
On entry: options for the optimization. There are four options to be set:
iprint controls iteration monitoring;
  if iprint0, then there is no printing of information else if iprint>0, then information is printed at every iprint iterations. The information printed consists of the value of FΨ at that iteration, the number of evaluations of FΨ, the current estimates of the communalities and an indication of whether or not they are at the boundary.
maxfun the maximum number of function evaluations.
acc the required accuracy for the estimates of ψi.
eps a lower bound for the values of ψ, see Section 3.
Let ε=machine precision then if IOP1=0, then the following default values are used:
  • iprint=-1
  • maxfun=100p
  • acc=10ε
  • eps=ε
If IOP10, then
  • iprint=IOP2
  • maxfun=IOP3
  • acc=10-l where l=IOP4
  • eps=10-l where l=IOP5
Constraint: if IOP10, IOPi must be such that maxfun1, εacc<1 and εeps<1, for i=3,4,5.
19:   IWK(4×NVAR+2) – INTEGER arrayWorkspace
20:   WK(LWK) – REAL (KIND=nag_wp) arrayWorkspace
21:   LWK – INTEGERInput
On entry: the dimension of the array WK as declared in the (sub)program from which G03CAF is called. The length of the workspace.
Constraints:
  • if MATRIX='D' or 'S', LWKmax5×NVAR×NVAR+33×NVAR-4/2,N×NVAR+7×NVAR+NVAR×NVAR-1/2;
  • if MATRIX='C', LWK5×NVAR×NVAR+33×NVAR-4/2.
22:   IFAIL – INTEGERInput/Output
On entry: IFAIL must be set to 0, -1​ or ​1. If you are unfamiliar with this parameter you should refer to Section 3.3 in the Essential Introduction for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value -1​ or ​1 is recommended. If the output of error messages is undesirable, then the value 1 is recommended. Otherwise, because for this routine the values of the output parameters may be useful even if IFAIL0 on exit, the recommended value is -1. When the value -1​ or ​1 is used it is essential to test the value of IFAIL on exit.
On exit: IFAIL=0 unless the routine detects an error or a warning has been flagged (see Section 6).

6  Error Indicators and Warnings

If on entry IFAIL=0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF).
Note: G03CAF may return useful information for one or more of the following detected errors or warnings.
Errors or warnings detected by the routine:
IFAIL=1
On entry,LDFL<NVAR,
orNVAR<2,
orNNVAR,
orNFAC<1,
orNVAR<NFAC,
orM<NVAR,
orMATRIX='D' or 'S' and LDX<N,
orMATRIX='C' and LDX<M,
orMATRIX'D', 'S' or 'C',
orMATRIX='D' or 'S' and WEIGHT'U' or 'W',
orIOP10 and IOP3 is such that maxfun<1,
orIOP10 and IOP4 is such that acc1.0,
orIOP10 and IOP4 is such that acc<machine precision,
orIOP10 and IOP5 is such that eps1.0,
orIOP10 and IOP5 is such that eps<machine precision,
orMATRIX='C' and LWK<5×NVAR×NVAR+33×NVAR-4/2,
orMATRIX='D' or 'S' and
LWK<max5×NVAR×NVAR+33×NVAR-4/2,N×NVAR+7×NVAR+ NVAR×NVAR-1/2.
IFAIL=2
On entry,WEIGHT='W' and a value of WT<0.0.
IFAIL=3
On entry, there are not exactly NVAR elements of ISX>0, or the effective number of observations NVAR.
IFAIL=4
On entry, MATRIX='D' or 'S' and the data matrix is not of full column rank, or MATRIX='C' and the input correlation/variance-covariance matrix is not positive definite.
This exit may also be caused by two of the eigenvalues of S* being equal; this is rare (see Lawley and Maxwell (1971)), and may be due to the data/correlation matrix being almost singular.
IFAIL=5
A singular value decomposition has failed to converge. This is a very unlikely error exit.
IFAIL=6
The estimation procedure has failed to converge in the given number of iterations. Change IOP to either increase number of iterations maxfun or increase the value of acc.
IFAIL=7
The convergence is not certain but a lower point could not be found. See E04LBF for further details. In this case all results are computed.

7  Accuracy

The accuracy achieved is discussed in E04LBF with the value of the parameter XTOL given by acc as described in parameter IOP.

8  Further Comments

The factor loadings may be orthogonally rotated by using G03BAF and factor score coefficients can be computed using G03CCF. The maximum likelihood estimators are invariant to a change in scale. This means that the results obtained will be the same (up to a scaling factor) if either the correlation matrix or the variance-covariance matrix is used. As the correlation matrix ensures that all values of ψi are between 0 and 1 it will lead to a more efficient optimization. In the situation when the data matrix is input the results are always computed for the correlation matrix and then scaled if the results for the covariance matrix are required. When you input the covariance/correlation matrix the input matrix itself is used and you are advised to input the correlation matrix rather than the covariance matrix.

9  Example

This example is taken from Lawley and Maxwell (1971). The correlation matrix for nine variables is input and the parameters of a factor analysis model with three factors are estimated and printed.

9.1  Program Text

Program Text (g03cafe.f90)

9.2  Program Data

Program Data (g03cafe.d)

9.3  Program Results

Program Results (g03cafe.r)


G03CAF (PDF version)
G03 Chapter Contents
G03 Chapter Introduction
NAG Library Manual

© The Numerical Algorithms Group Ltd, Oxford, UK. 2012