Introduction
This page makes no attempt to explain Factor Analysis comprehensively. It is assumed that users are already familiar with the concepts and purposes of Factor Analysis, and the explanations provided are to help users to make decisions on program parameters.
Only Exploratory Factor Analysis is covered. Confirmatory Factor Analysis is a separate subject and is not covered by StatsToDo
The page provides Factor Analysis using the following procedures
- Principal Components are calculated from Eigen Values and Vectors, using the Jacobi method
- Factor rotation allows the options of
- Orthogonal rotation (Varimax) produces factors that are uncorrelated with each other
- Oblique rotation (Oblimin, and in the R codes also Promax) produces a close fit, but allows the factors to be correlated
- Conversion of the rotated factors into the W matrix, which are coefficients to calculate factor scores
- Calculate factor scores if data is available
In addition, this page provides a discussion on sample size, and a program on Parallel Analysis. Paralle Analysis is a big subject and disccused in its own page in
ParallelAnalysis.php
Options and Decisions
The following options are available in using the programs on this page
Sample size is discussed in the sample size panel
Principal components based on the Covariance or the Correlation matrix.
- The Covariance matrix was used when Factor Analysis was first described, and it results in coefficients (W matrix) that can translate measurements directly into factor scores. The disadvantage is that the factor matrix is intuitively difficult to interpret, and it reflects the numerical values, so that length using inches or cms will result in very different factor values
- The Correlation matrix is currently the most common option, and conceptually it turns all measurements into standard deviate measurements (z = (value-mean)/SD). The disadvantage is that, to calculate factor scores, the primary measurements need to be converted to the z value, thus requiring an additional calculation.
- The program on this page uses the Correlation matrix by default when raw data is presented. However, there is an option to calculate starting with a matrix, and users can enter either the Covariance or the Correlation matrix
The number of factors to retain. Two options are available to the user
- Specify the number of factors to retain. This must be a number from 1 to the number ov variables
- Use the K1 rule where all factors with Eigen value >=1 are retained. This is done by specify the number of factors to retain as <1 or > number of variables
- Discussions on how to decide on which option is discussed in more details in ParallelAnalysis.php
Orthogonal (Varimax) or
Oblique (Oblimin) Rotation
- This is usually determined by the theoretical construct of the purpose of the exercise, whether the factors are correlated or not. For example, in the measurement of capability, the physical strength may be considered uncorrelated to mental acuity, but in measuring intelligence, the ability to comprehend and that of mathematics may be correlated.
- Without a preconceived framework, the oblique rotation can be tried first to obtained the best fit, and if the results show no meaningful correlation between the factors, then the orthogonal used to produce the final results.
Other Panels on This Page
The
Complete Program performs Factor Analysis in one step, but requires the user to specify tyoe of data entry, number of factors to retain, and type of rotation
The Subprograms performs the same functions, but in individual steps, allowing users to modify intermediary results
R codes contains the whole program in R
Parallel Analysis, a duplication of the Javascript program from ParallelAnalysis.php
Sample size, a discussion and some tables on how to determine the sample size necessary for Factor Analysis
References
Algorithms : It is difficult to find the algorithms for calculations
associated with Factor Analysis, as most modern text book and technical manuals
advise users to use one of the commercial packages. I have eventually found some
useful algorithms in old text books, and they are as follows
Press WH, Flannery VP, Teukolsky SA, Vetterling WT (1989). Numerical Recipes in Pascal.
Cambridge University Press IBSN 0-521-37516-9 p.395-396 and p.402-404.
Jacobi method for finding Eigen values and Eigen vectors
Norusis MJ (1979) SPSS Statistical Algorithms Release 8. SPSS Inc Chicago
- p. 86 for converting Eigen values and vectors to Principal Components
- p. 91-93 for Varimax Rotation
- p. 94-97 for Oblimin Rotation
- p. 97-98 for Factor scores
Text books I learnt Factor Analysis some time ago, so all my text books are old,
but they are adequate in explaining the basic concepts and provide the calculations used in these pages.
Users should search for better and newer text books.
- Thurston LL (1937) Multiple Factor Analysis. University of Chicago Press. I have not read this
book, but this is quoted almost universally as it is the original Factor
Analysis text book which set out the Principles and procedures
- Gorsuch RL (1974) Factor Analysis. W. B. Saunders Company London 1974 ISBN 0-7216-4170-9
A standard text book teaching Factor Analysis at the Master's level. This is my
copy and I believe there are later additions of this book available
Parallel Analysis
Sample Size
Mundfrom DJ, Shaw DG, Tian LK (2005) Minimum sample size recommendations for conducting factor analysis.
International Journal of Testing 5:2:p 159-168
Free Factor Analysis software and its user manual can be downloaded from
http://psico.fcep.urv.es/utilitats/factor/Download.html. This is a package for Windows written by Drs. Lorezo-Seva and Ferrendo from Universitat Rovira i Virgili in Terragona in Spain. The presentation and options are very similar to that from SPSS, and the manual is excellent. The best part is that it is free and yes it is in English.
Teaching and discussion papers on the www There is an enormous list of
discussion papers, technical notes and tutorials on the www that can be easily found by Google search.
The following is a small sample of this.
This panel provides explanations for each section of the program for Factor Analysis. The program is divided into 2 sections. Section 1 contains all the subroutine functions, and section 2 the main program.
Please note that user input is required in 2 places. Firstly the data input, and secondly the choice of method to determine the number of factors to retain.
Section 1. Functions
Four (4) functions are presented here
# create eigen value, eigen vector and principal components from correlation matrix
CorToPComp <- function(corMx)
{
eigenResults <- eigen(corMx) # Calculate eigen value and vector
evl <- eigenResults$values # evl = eigen value
evc <- eigenResults$vectors # evc = eigen vector
nv = length(evl)
evm <- matrix(data=0, nrow = nv, ncol=nv) # evm = principal component
for(j in 1:nv)
{
x = sqrt(evl[j])
evm[,j] = evc[,j] * x
}
list("evl"= evl, "evc"= evc, "evm"= evm) # returns eigen values, eigen vectors and prinComp
}
The
CorToPComp function takes a correlation matrix (corMx), and returns a class containing 3 structures. These are
- evl, the eigen value array
- evc, the eigen vector matrix
- evm, the principal component matrix
# Create eigen value, eigen vector and principal components from data matrix
DataFrameToPCom <- function(df)
{
cMx <- cor(df)
print(cMx) # correlation matrix
res <- CorToPComp(cMx) # returns eigen values, eigen vectors and prinComp
}
The
DataFrameToPCom function takes the dataframe containing a matrix of measurements (df), with rows representing cases and column representing variables. From this, the correlation matrix (cMx) is calculated, and presents cMx to the
CorToPComp function to receive the eigen value, vector, and principle component. It then returns these 3 structures.
# Calculate the number of factors to retain by K1 rule
NumberOfFactorsByK1 <- function(arEigenValue)
{
f = 0
for(i in 1:length(arEigenValue))
{
if(arEigenValue[i]>=1)
{
f = f + 1
}
}
f # return number of factors
}
The
NumberOfFactorsByK1 function takes the eigen value array, and counts the number of values that are >=1, and returns that value as the number of factors to retain
# Calculate the number of factors to retain by Parallel Analysis
NumberOfFactorsByParallel <- function(arEigenValue, ssiz, ite) # ite=number of iterations
{
nc = ssiz #sample size, number of rows of data
nv = length(arEigenValue) #number of variables
# parallel analysis begin
ex <- rep(0, nv)
exx <- rep(0, nv)
for(i in 1:ite)
{
dMx <- replicate(nv, rnorm(nc)) # random data matrix
pComp <- eigen(cor(dMx))$values # array of eigen values
for(j in 1:nv)
{
ex[j] = ex[j] + pComp[j]
exx[j] = exx[j] + pComp[j]^2
}
}
resAr <- rep(0,nv)
for(j in 1:nv)
{
mean = ex[j] / ite # mean
sd = sqrt((exx[j] - ex[j]^2/ite)/(ite-1)) # SD
lim = mean + qnorm(0.95) * sd # 95 percentile
resAr[j] = lim
}
print(resAr) # array of 95 percentile eigen values
print(arEigenValue)
f = 0
for(i in 1:nv)
{
if(arEigenValue[i]>=resAr[i])
{
f = f + 1
}
}
f
}
The
NumberOfFactorsByParallel function accepts 3 parameters
- arEigenValue is the array of eigen value calculated from the data or correlation matrix
- ssiz is the sample size (number of rows) of the data
- ite is the number of iterations to be used to calculate the 95the percentile. This should be 100 or more, usually 1000 is sufficient
The program then iterates, calculating eigen values using normally distributed random data and establishes the 95th percentile values by simulation
The simulated values are then compared with the input eigen value array. The number of eigen value from the input array that is >= to that in the simulated array is counted, and returned as the number of factors to retain.
Section 2. Main Program
The code in section 2 are executed in the order they present
# Step 1. Data input
myDat = ("
0.178 -0.338 0.470 1.211 0.508 0.354
0.304 -0.230 -0.981 -1.218 -0.337 -0.780
1.155 0.350 -0.310 -0.083 -2.224 -1.502
0.562 0.148 1.376 0.616 0.887 1.342
......
")
myDataFrame <- read.table(textConnection(myDat),header=FALSE)
summary(myDataFrame)
The data is a table (matrix), without header. The rows representing cases and column variables. A summary of the data is presented as follows
V1 V2 V3 V4 V5 V6
Min. :-1.5520 Min. :-0.88300 Min. :-1.6040 Min. :-2.06200 Min. :-2.2240 Min. :-1.75400
1st Qu.:-0.1710 1st Qu.:-0.33800 1st Qu.:-0.9810 1st Qu.:-0.42100 1st Qu.:-1.1100 1st Qu.:-0.61600
Median : 0.1780 Median : 0.14100 Median : 0.2950 Median : 0.15700 Median :-0.4100 Median :-0.17100
Mean : 0.1406 Mean : 0.07696 Mean : 0.1205 Mean : 0.01056 Mean :-0.3199 Mean :-0.09372
3rd Qu.: 0.6360 3rd Qu.: 0.39100 3rd Qu.: 1.2760 3rd Qu.: 0.39000 3rd Qu.: 0.5080 3rd Qu.: 0.66200
Max. : 1.4290 Max. : 1.35600 Max. : 1.7680 Max. : 1.26700 Max. : 2.2750 Max. : 1.43000
# Step 2. Calculate Principal Component Matrix
nc = nrow(myDataFrame) # number of cases
nv = ncol(myDataFrame) # number of variables
res <- DataFrameToPCom(myDataFrame) # eigen values, eigen vector, Principal components
res$evl # eigen values
#res$evc # eigen vector
pComp <- res$evm # All Principal components
pComp # All Principal components
Step 2 calculate the eigen value, vector, and the Principal Component. It prints out the correlation matrix, eigen value array and the Principal Component matrix, as follows
Correlation matrix
V1 V2 V3 V4 V5 V6
V1 1.0000000 0.65568487 0.2313851 0.27914974 0.1181262 0.2401638
V2 0.6556849 1.00000000 0.2329578 0.05241047 0.2292866 0.3559332
V3 0.2313851 0.23295783 1.0000000 0.65596377 0.2608266 0.2263553
V4 0.2791497 0.05241047 0.6559638 1.00000000 0.3849751 0.3061149
V5 0.1181262 0.22928662 0.2608266 0.38497508 1.0000000 0.7628703
V6 0.2401638 0.35593319 0.2263553 0.30611491 0.7628703 1.0000000
eigen value array
[1] 2.6752733 1.3097426 1.1512484 0.4477224 0.2252526 0.1907607
Principal Components Matrix
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] -0.5959219 0.6554206 0.22250426 -0.3505435042 0.06255386 0.19748157
[2,] -0.6026316 0.6998645 -0.05526649 0.2888334979 -0.14161525 -0.20122309
[3,] -0.6440947 -0.2740522 0.57047741 0.3895179069 0.08849390 0.15823263
[4,] -0.6738459 -0.4272552 0.46373143 -0.2988512091 -0.07734747 -0.23031085
[5,] -0.7223140 -0.3366373 -0.50819274 -0.0194946053 -0.28227127 0.16315911
[6,] -0.7525434 -0.1392112 -0.54766295 -0.0007828776 0.32831983 -0.08105191
# USER INPUT REQUIRED
# Step 3. User to determine number of factors to retain (one of 3 options)
#nf = 3 # choice a specify a number of factors
nf = NumberOfFactorsByK1(res$evl) # choice b using k1 rule
#nf = NumberOfFactorsByParallel(res$evl, nc, 1000) # choice c using parallel analysis
nf # number of factors to retain
# END USER INPUT
Step 3 is for the user to choose the method of determining the number of factors to retain. Three options are available
- nf=n, the user can arbitrarily stipulate the number of factor to retain (e.g. 3)
- nf = NumberOfFactorsByK1(res$evl), calls the k1 function with the eigen value array
- nf = NumberOfFactorsByParallel(res$evl, nc, ite), calls the Parallel Analysis function with the eigen value array, the number of cases (row, sample size), and the number of iterations to use (e.g. 1000)
The user can choose by commenting (place # as first character) or uncommenting (remove #) on the relevant row of command. In the default code, the K1 rule is chosen and it returns 3 factors to retain, as follows
> nf # number of factors to retain
[1] 3
# Step 4. Create z scores required for calculating factor values for each case
arMean <- rep(0,nv)
arSD <- rep(0,nv)
zMx <- matrix(0,nrow=nc, ncol=nv)
for(i in 1:nv)
{
arMean[i] = mean(myDataFrame[ , i])
arSD[i] = sd(myDataFrame[ , i])
zMx[, i] <- (myDataFrame[ , i] - arMean[i]) / arSD[i]
}
arMean # array of means
arSD # array of SDs
zMx # matrix of z values
Step 4 calculates the mean and Standard Deviation for each variable (column). It then transform the values into standardized z values (z=(value-mean) / SD). These are used later for calculating factor scores. The results are as follows
> arMean # array of means
[1] 0.14064 0.07696 0.12048 0.01056 -0.31988 -0.09372
> arSD # array of SDs
[1] 0.7599281 0.5908474 1.1104260 0.7961356 1.0809935 0.9073096
> zMx # matrix of z values
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0.04916255 -0.70231326 0.314762090 1.507833555 0.76585104 0.493458903
[2,] 0.21496772 -0.51952497 -0.991943656 -1.543154170 -0.01583728 -0.756390101
[3,] 1.33481056 0.46211590 -0.387671047 -0.117517666 -1.76145370 -1.552149344
[4,] 0.55447354 0.12023408 1.130665195 0.760473449 1.11645445 1.582392604
......
# Step 5. Calculate Factor Score if there is only 1 Principal Component
if(nf<2)
{
pCompMx <- matrix(pComp[,1:nf])
print("Principal Component")
print(pCompMx) # principal component matrix to be used for subsequent rotation and processing
# Calculate coefficient matrix for scores
coeffMx <- pCompMx %*% solve(t(pCompMx) %*% pCompMx)
print("Coefficients")
print(coeffMx) # coefficient matrix
# Calculate Factor Scores
scoreMx <- zMx %*% coeffMx
print("Factor Scores")
print(scoreMx) # factor scores for varimax
} else
Step 5 is performed if there is only one (1) Principal component retained. The principal component is translated into the coefficient matrix, and the factor scores are calculated. As this is not the option in the example code, no output is shown. However, if the number of factor was chosen as 1 or if the Parallel Analysis was used, then the following results is shown
[1] "Principal Component"
[,1]
[1,] -0.5959219
[2,] -0.6026316
[3,] -0.6440947
[4,] -0.6738459
[5,] -0.7223140
[6,] -0.7525434
[1] "Coefficients"
[,1]
[1,] -0.2227518
[2,] -0.2252598
[3,] -0.2407585
[4,] -0.2518793
[5,] -0.2699963
[6,] -0.2812959
[1] "Factor Scores"
[,1]
[1,] -0.65390672
[2,] 0.91369654
[3,] 0.63370706
......
# Step 6 Performs Factor rotation if there is more than 1 principal component
{
pCompMx <- pComp[,1:nf]
print("Principal Components")
print(pCompMx) # principal component matrix to be used for subsequent rotation and processing
# perform Varimax rotation
vMax <- varimax(pCompMx)
#vMax
vMx <- pCompMx %*% vMax$rotmat
print("Varimax Factor Loadings")
print(vMx) #varimax loading matrix
# Calculate coefficient matrix for scores
vCoeffMx <- vMx %*% solve(t(vMx) %*% vMx)
print("Coefficient Matrix")
print(vCoeffMx) # varimax coefficient matrix
# Calculate Factor Scores
vScoreMx <- zMx %*% vCoeffMx
print("Factor Scores")
print(vScoreMx) # factor scores for varimax
# Perform Promax rotation
pMax <- promax(pCompMx)
#pMax
pMx <- pCompMx %*% pMax$rotmat
print("Promax Factor Loadings")
print(pMx) # promax loading matrix
# Calculate coefficient matrix for scores
pCoeffMx <- pMx %*% solve(t(pMx) %*% pMx)
print("Coefficient Matrix")
print(pCoeffMx) # promax coefficient matrix
# Calculate Factor Scores
pScoreMx <- zMx %*% pCoeffMx
print("Factor Scores")
print(pScoreMx) # factor scores for promax
# Performs Oblimin Rotation
#install.packages("GPArotation") # if not already installed
library(GPArotation)
obmn <- oblimin(pCompMx, normalize=TRUE)
print("Oblimin Factor Loadings and Factor Correlations")
print(obmn)
oMx <- obmn$loadings
#oMx # oblimin loading matrix (factor pattern)
# Calculate coefficient matrix for scores
oCoeffMx <- oMx %*% solve(t(oMx) %*% oMx)
print("Coefficient Matrix")
print(oCoeffMx) # oblimin coefficient matrix
# Calculate Factor Scores
oScoreMx <- zMx %*% oCoeffMx
print("Factor Scores")
print(oScoreMx) # factor scores for oblimin
}
Section 6 is executed if more than 1 Principal component is retained. The retained factors are subjected to rotation. The rotated factors are then converted to the coefficient matrix, and the factor scores are calculated using the z matrix and the coefficient matrix.
Three types of rotations are carried out, Varimax, Promax, and Oblimin. The results have the same structure but different values, and are shown separately as follows
Three 3 factors have eigen values >=1, and these are retained, as follows.
[1] "Principal Components"
[,1] [,2] [,3]
[1,] -0.5959219 0.6554206 0.22250426
[2,] -0.6026316 0.6998645 -0.05526649
[3,] -0.6440947 -0.2740522 0.57047741
[4,] -0.6738459 -0.4272552 0.46373143
[5,] -0.7223140 -0.3366373 -0.50819274
[6,] -0.7525434 -0.1392112 -0.54766295
The first rotation is Orthogonal (Varimax) and the results are as follows
[1] "Varimax Factor Loadings"
[,1] [,2] [,3]
[1,] -0.02453600 0.88913790 0.207458613
[2,] -0.21308829 0.90034535 -0.001106254
[3,] -0.08583428 0.16383896 0.883853119
[4,] -0.22653935 0.03792675 0.893814143
[5,] -0.92160147 0.03512097 0.206734722
[6,] -0.90981342 0.21391356 0.110075392
[1] "Coefficient Matrix"
[,1] [,2] [,3]
[1,] 0.14559564 0.560732513 0.04210757
[2,] -0.01747185 0.564782857 -0.13890447
[3,] 0.13675107 0.001150852 0.57323095
[4,] 0.02871355 -0.102217672 0.56642661
[5,] -0.56263156 -0.125188522 -0.03989932
[6,] -0.54909048 0.004753418 -0.12330412
[1] "Factor Scores"
[,1] [,2] [,3]
[1,] -0.59607769 -0.616382588 1.04273054
[2,] 0.28465380 -0.017896153 -1.26758238
[3,] 1.97319968 1.234168820 -0.03510730
[4,] -1.24194545 0.170139383 0.84586987
.......
The second rotation is Oblique (Promax) and the results are as follows
[1] "Promax Factor Loadings"
[,1] [,2] [,3]
[1,] 0.10895336 0.90058699 0.13062313
[2,] -0.12048997 0.91037793 -0.11768389
[3,] 0.06214795 0.07796628 0.90044284
[4,] -0.09963805 -0.06950278 0.90178181
[5,] -0.93819973 -0.08681418 0.07649840
[6,] -0.91956666 0.10906335 -0.04013362
[1] "Coefficient Matrix"
[,1] [,2] [,3]
[1,] 0.06917274 0.53965839 0.07797890
[2,] -0.06418369 0.54467281 -0.07464688
[3,] 0.04806237 0.04494003 0.54263501
[4,] -0.04451634 -0.04438377 0.54157032
[5,] -0.53110712 -0.05804648 0.03337254
[6,] -0.52092238 0.05964160 -0.03703614
[1] "Factor Scores"
[,1] [,2] [,3]
[1,] -0.66731997 -0.42380216 1.05094056
[2,] 0.47166717 -0.18724198 -1.29096075
[3,] 1.79334119 0.96951114 -0.20571498
[4,] -1.36613427 0.41134414 1.03830371
.......
The third rotation is Oblique (Oblimin) and the results are as follows. Please note that the matrix Phi is the correlation matrix between the rotated factors. Rotation matrix is an intermediary result produced by R and can be ignored.
[1] "Oblimin Factor Loadings and Factor Correlations"
Oblique rotation method Oblimin Quartimin converged.
Loadings:
[,1] [,2] [,3]
[1,] 0.0946 0.8959 0.1375
[2,] -0.1316 0.9056 -0.1080
[3,] 0.0514 0.0822 0.8968
[4,] -0.1078 -0.0635 0.8980
[5,] -0.9345 -0.0804 0.0820
[6,] -0.9174 0.1137 -0.0324
Rotating matrix:
[,1] [,2] [,3]
[1,] 0.534 -0.418 -0.461
[2,] 0.339 0.944 -0.487
[3,] 0.856 0.126 0.817
Phi:
[,1] [,2] [,3]
[1,] 1.000 -0.234 -0.292
[2,] -0.234 1.000 0.207
[3,] -0.292 0.207 1.000
[1] "Coefficient Matrix"
[,1] [,2] [,3]
[1,] 0.07318699 0.54287847 0.07602563
[2,] -0.06171397 0.54743794 -0.07871153
[3,] 0.05236099 0.04092966 0.54530484
[4,] -0.04109614 -0.05017189 0.54376227
[5,] -0.53316646 -0.06611237 0.02825672
[6,] -0.52274532 0.05300352 -0.04301873
[1] "Factor Scores"
[,1] [,2] [,3]
[1,] -0.66482370 -0.44502856 1.05097451
[2,] 0.46311661 -0.16992759 -1.29069359
[3,] 1.80422927 1.00183304 -0.19319602
[4,] -1.36133404 0.38501672 1.02623918
...........