Principal Components Analysis with Rotation
A Python Program
Classical theory of principal components
analysis (PCA) explains that principal components are constructed one by one
using the criterion of maximum variance of each component, so rotation for the
set of principal components is not considered. However, when principal
components analysis is considered from the point of view of orthogonal
projection of the data (Okamoto, 2006),
rotation is a natural method of analysis. The orthogonal projection P projects
the data Z in p-dimensional space 
, where p is the number of the variables in the data, onto
q-dimensional subspace 
of so that the projected
data ZP has maximum variance in q-dimensional space. This q-dimensional
subspace, onto which the data Z is projected, can has any coordinate basis,
that is, a coordinate basis can be chosen without any mathematical constraints.
We can choose the basis that may give us easy interpretation of the results
from the PCA. Popular methods for choice of a coordinate basis are varimax and
promax rotations. The Python program for PCA in this website uses varimax and
promax rotations.
This multi-dimensional projective view of
PCA is also presented by Deisenroth, Faisal, and Ong (2020), “Mathematics for Machine
Learning.” Cambridge
University Press, chapter 10, in which a projection matrix P is denoted by 
. Another notational difference is that although Deisenroth et al.
denotes the data matrix by the variable-case form, Okamoto (2006) denotes the data
matrix by the case-variable form.
The program files and a sample data file are archived in the file pcaprg.zip, which can be downloaded by clicking the file name pcaprg.zip. The program and the data in the file can be freely used on the user’s responsibility, although all rights are reserved. Listing 1 and Listing 2 of the main program and the class for PCA are shown at the end of this website.
It is easy to install Python. How to install Python is explained in the website
http://y-okamoto-psy1949.la.coocan.jp/Python/en/BeginningPython/
When you run the program MainPrg.py of Listing 1, the name of the input data file will be required (Figure 1).
Figure 1
An input data file should be made as a CSV file, and follow the format shown in Figure 2.
Figure 2 DataV6.csv (Opened by Excel)
The data file DataV6.csv shown in Figure 1 was created by Excel choosing the file extension “.csv”. In the first row, names of variables are written in the order of case identification variable, then the variables to be analyzed. From the second row, values of the variables of each case are written in one row in the columns under the respective variable names. A value of the case identification variable is treated as string type, and values of variables to be analyzed by PCA are treated as numerical type (float type).
When the Enter key is pressed after setting the name of the input data file, a name for the output file is required (Figure 3).。
Figure 3
Set any name of a text file for the output. Press the Enter key, then calculation starts. Eigen values as lambdas are displayed (Figure 4).
Figure 4
For the meaning of lambda values, see Okamoto (2006). Following display of lambda values, the number of principal components is asked. When specific reasons for the number of components are not considered, the number of lambda values larger than 1 is usually chosen. In Figure 5, 2 is set.
Figure 5
After setting the number of components, a value for xdim will be required (Figure 5). The value for xdim denotes the component for the abscissa in a two-dimensional map of the pattern matrix. Set the integer number, which you desire. If the integer value less than 1 is set, drawing of a map is passed and the program goes to the next step. When the integer equal to or larger than 1 is set, the value for ydim, which denotes the component for the ordinate, will be required (Figure 6).
Figure 6
In Figure 6, 2 is set for ydim. After setting the value for ydim, the map corresponding to the selected components will be displayed (Figure 7).
Figure 7
In Figure 7, variables are gathered in two places, which are far from the axes. This configuration is typical for the initial principal components. The image of the map can be saved by clicking the floppy-disk icon at the bottom of the window. The window is closed by clicking the X icon of the window.
After the window is closed, the next value for xdim is required. If the value less than 1 is set, display of the map for the initial principal components ends and the next step of showing the map for the varimax rotated pattern starts.
The value for xdim for the map of the pattern after the varimax rotation is required as is shown in Figure 8.
Figure 8
The values for xdim and ydim are required in the same way as in the case of the initial principal components. In Figure 9, xdim = 1 and ydim = 2 are set.
Figure 9
The map in Figure 10 is the corresponding map.
Figure 10
Variables are gathered in two groups, each of which is near the abscissa or the ordinate,
When the value less than 1 is set for xdim, display of the map for the varimax rotated pattern ends and calculation of the promax rotation starts. After the calculation ends and the results of the analysis are saved, a message saying that the output file was saved is displayed and the program ends (Figure 11).
Figure 11
The contents of the output file for the input data file DataV6.txt (Figure 2) is as follows.
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ Results.txt ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
Input data file =
DataV6.csv
n_vars = 6
Var_names...
CaseID P-1 P-2 P-3 P-4 P-5 P-6
1 46.000 57.000 53.000 65.000 54.000 67.000
2 19.000 16.000 38.000 34.000 26.000 40.000
3 28.000 30.000 23.000 20.000 35.000 18.000
4 40.000 44.000 55.000 48.000 43.000 45.000
5 64.000 62.000 67.000 57.000 71.000 70.000
6 49.000 43.000 61.000 61.000 52.000 59.000
7 27.000 36.000 57.000 64.000 30.000 50.000
8 52.000 48.000 46.000 52.000 44.000 50.000
9 66.000 65.000 57.000 57.000 63.000 48.000
10 54.000 52.000 26.000 34.000 55.000 26.000
11 62.000 58.000 55.000 65.000 54.000 65.000
12 81.000 83.000 54.000 51.000 84.000 54.000
13 32.000 36.000 19.000 33.000 36.000 26.000
14 81.000 80.000 34.000 27.000 69.000 31.000
15 32.000 41.000 71.000 56.000 41.000 58.000
16 70.000 70.000 29.000 24.000 66.000 30.000
17 48.000 45.000 35.000 34.000 46.000 38.000
18 56.000 49.000 46.000 46.000 49.000 35.000
19 66.000 65.000 71.000 72.000 66.000 57.000
20 55.000 64.000 60.000 69.000 63.000 47.000
Lambda...
Lambd cum.sqr
%
Lambda-1 =
1.85921 3.45665 57.61
Lambda-2 =
1.47763 5.64003 94.00
Lambda-3 =
0.39244 5.79404 96.57
Lambda-4 =
0.34309 5.91175 98.53
Lambda-5 =
0.22609 5.96287 99.38
Lambda-6 =
0.19269 6.00000 100.00
Pattern matrix for the
principal components =
comp.1 comp.2
P-1 0.76686 0.61656
P-2 0.81292 0.55673
P-3 0.72398 -0.63147
P-4 0.68197 -0.66907
P-5 0.83447 0.51828
P-6 0.72260 -0.61503
Principal
components...
CaseID comp.1 comp.2
1 0.61232 -0.82129
2 -1.76651 -0.94212
3 -1.99807 0.54979
4 -0.33724 -0.57566
5 1.32511 -0.33666
6 0.35849 -0.93765
7 -0.53145 -1.57079
8 -0.10764 -0.30908
9 0.81074 0.22764
10 -0.68697 1.15484
11 0.83043 -0.54224
12 1.57715 1.14908
13 -1.61827 0.38890
14 0.38392 2.15026
15 -0.01195 -1.49276
16 -0.07016 1.93720
17 -0.71796 0.40174
18 -0.25206 0.26182
19 1.36769 -0.45348
20 0.83244 -0.23956
Varimax rotation was
applied...
Pattern for the
Varimax rotated components =
comp.1 comp.2
P-1 0.98285 -0.04715
P-2 0.97765 -0.12248
P-3 0.12304 -0.95277
P-4 0.06666 -0.95305
P-5 0.96827 -0.16555
P-6 0.13291 -0.93955
Sorted Pattern for
Varimax rotation =
comp.1 comp.2
P-1 0.98285 -0.04715
P-2 0.97765 -0.12248
P-5 0.96827 -0.16555
P-4 0.06666 -0.95305
P-3 0.12304 -0.95277
P-6 0.13291 -0.93955
Components for
Varimax...
CaseID comp.1 comp.2
1 -0.08642 -1.02077
2 -1.94693 0.46648
3 -1.13077 1.73664
4 -0.63419 -0.20717
5 0.76846 -1.13080
6 -0.35357 -0.93952
7 -1.43953 -0.82315
8 -0.28555 -0.15994
9 0.75775 -0.36733
10 0.25177 1.31993
11 0.26190 -0.95658
12 1.94247 -0.18600
13 -0.95323 1.36433
14 1.71343 1.35468
15 -0.99897 -1.10929
16 1.23228 1.49638
17 -0.27090 0.77684
18 -0.01500 0.36313
19 0.72285 -1.24647
20 0.46413 -0.73138
Promax rotation was
applied...
Pattern for the Promax
rotated components =
comp.1 comp.2
P-1 0.99665 0.06582
P-2 0.98263 -0.01159
P-3 0.01527 -0.95714
P-4 -0.04225 -0.96396
P-5 0.96809 -0.05659
P-6 0.02686 -0.94252
Correlation matrix of
the components...
comp.1 comp.2
comp.1 1.00000 -0.22419
comp.2 -0.22419 1.00000
Sorted Pattern for
Promax rotation =
P-1 0.99665 0.06582
P-2 0.98263 -0.01159
P-5 0.96809 -0.05659
P-4 -0.04225 -0.96396
P-3 0.01527 -0.95714
P-6 0.02686 -0.94252
Components for
Promax...
CaseID comp.1 comp.2
1 0.02941 -1.00453
2 -1.98715 0.68293
3 -1.31965 1.85301
4 -0.60674 -0.13438
5 0.89125 -1.21021
6 -0.24521 -0.89369
7 -1.33736 -0.65567
8 -0.26566 -0.12673
9 0.79439 -0.45039
10 0.10110 1.28314
11 0.36824 -0.98001
12 1.95105 -0.40373
13 -1.10120 1.46307
14 1.54949 1.15295
15 -0.86731 -0.98964
16 1.05542 1.34797
17 -0.35689 0.80242
18 -0.05591 0.36250
19 0.85899 -1.32000
20 0.54376 -0.77903
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
If the data contains categorical variables, these categorical data may be quantified, and then can be analyzed by multivariate analysis of quantitative variables. A method for quantification is explained in the website, in that the associated Python program is also uploaded.
The following listing (Listing 1) is the main program file MainPrg.py for PCA discussed in this website. The class InitPCA, by which PCA is done, is shown in Listing 2.
import csv
import numpy
import math
import mymodule as mm
from init_pca import
InitPCA
from rotate import
Varimax, Promax
from PlotComp import
PlotPattern
from ArrangePat import
ArrangePtn
"""
Principal Components Analysis with Rotation of Components
Yasuharu Okamoto, 2016.10, 2021.01
"""
#
# Prepare
the input data file (CSV format)
#
s = input("Input
data file (*.csv) = ")
with open(s, 'r') as
f:
data = [v for v in
csv.reader(f)]
#
# Prepare
the output file
#
s_out =
input("Output trxt file (*.txt) = ")
f_out = open(s_out,
"w")
f_out.write("Input
data file = " + s + "\n\n")
#print('data:\n',
data)
#
#
Calculate the principal components
#
calc_init_pca =
InitPCA(data, f_out)
var_names, n_data,
n_vars, CaseID, n_comp, Pat, PComp = calc_init_pca.calc_pca()
"""
var_names:
A list of variable names
n_data: The
number of cases in data
n_vars: The
number of variables
CaseID: The
variable for case idenfication (string type)
n_com:
The number of componets
Pat:
The pattern matrix of the principal components
PComp:
The matrix of the principal components
"""
f_out.write("\n\nPrincipal
components...\n")
f_out.write("%10s"%var_names[0])
for j in range(0,
n_comp):
f_out.write("%10s"%("comp.%d"%(j + 1)))
f_out.write("\n")
for i in range(0,
n_data):
f_out.write("%10s"%CaseID[i])
for j in range(0, n_comp):
f_out.write("%10.5f"%PComp[i, j])
f_out.write("\n")
PlotPtn = PlotPattern(
Pat, var_names, n_vars, n_comp, " ...principal components" )
while True:
print("\nSet the
component number (from 1 to %d) to plot the pattern."%n_comp)
print("If you do not
want to plot the pattern, set the number less than 1.")
xdim = int(input("xdim
= "))
if xdim < 1:
break
ydim = int(input("ydim
= "))
PlotPtn.Plot(xdim, ydim)
rot = Varimax(Pat,
n_vars, n_comp)
PatV = rot.rotateV()
f_out.write("\n\nVarimax
rotation was applied...\n\n")
f_out.write("\nPattern
for the Varimax rotated components =\n")
f_out.write("%10s"%"
")
for j in range(0,
n_comp):
f_out.write("%10s"%("comp.%d"%(j + 1)))
f_out.write("\n")
for i in range(0,
n_vars):
f_out.write("%10s"%var_names[i + 1])
for j in range(0, n_comp):
f_out.write("%10.5f"%PatV[i, j])
f_out.write("\n")
ArgPtn =
ArrangePtn(var_names, PatV, n_vars, n_comp)
PtnArgd =
ArgPtn.arrange()
f_out.write("\n\nSorted
Pattern for Varimax rotation =\n")
f_out.write("%10s"%"
")
for j in range(0,
n_comp):
f_out.write("%10s"%("comp.%d"%(j + 1)))
f_out.write("\n")
for i in range(0,
n_vars):
f_out.write("%10s"%PtnArgd[i][0])
for j in range(0, n_comp):
f_out.write("%10.5f"%PtnArgd[i][1][j])
f_out.write("\n")
PlotPtn = PlotPattern(
PatV, var_names, n_vars, n_comp, " ...Varimax rotated components" )
while True:
print("\nSet the
component number (from 1 to %d) to plot the pattern for the Varimax
criterion."%n_comp)
print("If you do not
want to plot the pattern, set the number less than 1.")
xdim = int(input("xdim
= "))
if xdim < 1:
break
ydim = int(input("ydim
= "))
PlotPtn.Plot(xdim, ydim)
PPP = PatV *
numpy.linalg.inv(numpy.transpose(PatV) * PatV)
PCompV = PComp *
numpy.transpose(Pat) * PPP
f_out.write("\n\nComponents
for Varimax...\n")
f_out.write("%10s"%var_names[0])
for j in range(0,
n_comp):
f_out.write("%10s"%("comp.%d"%(j
+ 1)))
f_out.write("\n")
for i in range(0,
n_data):
f_out.write("%10s"%CaseID[i])
for j in range(0, n_comp):
f_out.write("%10.5f"%PCompV[i, j])
f_out.write("\n")
obrot = Promax(PatV,
n_vars, n_comp, 4.0)
Patob, Phi =
obrot.rotateP();
f_out.write("\n\nPromax
rotation was applied...\n\n")
f_out.write("\nPattern
for the Promax rotated components =\n")
f_out.write("%10s"%"
")
for j in range(0,
n_comp):
f_out.write("%10s"%("comp.%d"%(j + 1)))
f_out.write("\n")
for i in range(0,
n_vars):
f_out.write("%10s"%var_names[i + 1])
for j in range(0, n_comp):
f_out.write("%10.5f"%Patob[i, j])
f_out.write("\n")
f_out.write("\n\nCorrelation
matrix of the components...\n")
f_out.write("%10s"%"
")
for j in range(0,
n_comp):
f_out.write("%10s"%("comp.%d"%(j + 1)))
f_out.write("\n")
for i in range(0,
n_comp):
f_out.write("%10s"%("comp.%d"%(i + 1)))
for j in range(0, n_comp):
f_out.write("%10.5f"%Phi[i, j])
f_out.write("\n")
ArgPtn =
ArrangePtn(var_names, Patob, n_vars, n_comp)
PtnArgd =
ArgPtn.arrange()
f_out.write("\n\nSorted
Pattern for Promax rotation =\n")
for i in range(0,
n_vars):
f_out.write("%10s"%PtnArgd[i][0])
for j in range(0, n_comp):
f_out.write("%10.5f"%PtnArgd[i][1][j])
f_out.write("\n")
PPP = Patob *
numpy.linalg.inv(numpy.transpose(Patob) * Patob)
PCompob = PComp *
numpy.transpose(Pat) * PPP
f_out.write("\n\nComponents
for Promax...\n")
f_out.write("%10s"%var_names[0])
for j in range(0,
n_comp):
f_out.write("%10s"%("comp.%d"%(j + 1)))
f_out.write("\n")
for i in range(0,
n_data):
f_out.write("%10s"%CaseID[i])
for j in range(0, n_comp):
f_out.write("%10.5f"%PCompob[i, j])
f_out.write("\n")
f_out.close()
print(s_out + "
was saved.")
Listing
2. The class InitPCA for principal components analysis
import numpy
import math
import scipy.stats as ss
import mymodule as mm
class InitPCA:
"""
Calculation of the Principal Components
"""
def __init__(self, data,
f_out):
self.data = data
self.f_out = f_out
def calc_pca(self):
data
= numpy.array(self.data)
n_vars = len(data[0]) - 1
#
#
var_names: List of the variable names
#
var_names
= data[0]
n_data = len(data) - 1
CaseID = data.T[0][1:]
X =
[]
for
i in range(n_vars):
X.append([float(v) for v in data.T[i + 1][1:]])
self.f_out.write("\n")
self.f_out.write("n_vars = %d\n"%n_vars)
self.f_out.write("Var_names...\n")
for
vn in var_names:
self.f_out.write("%10s"%vn)
self.f_out.write("\n")
self.f_out.write("\n")
for
i in range(0, n_data):
self.f_out.write("%10s"%CaseID[i])
for j in range(0, n_vars):
self.f_out.write("%10.3f"%X[j][i])
self.f_out.write("\n")
R =
numpy.corrcoef(X)
Lmbd2, V = numpy.linalg.eig(R)
Lmbd
= []
for i
in range(0, len(Lmbd2)):
if Lmbd2[i] > 0.0:
Lmbd.append(math.sqrt(Lmbd2[i]))
else:
Lmbd.append(0.0)
sort_index = numpy.argsort(numpy.array(Lmbd) * (-1))
Lmbd
= numpy.array(Lmbd)[sort_index]
n_eign = 0
for
v in Lmbd:
if v > 0.0:
n_eign += 1
ListLV = V.T[sort_index].T
self.f_out.write("\nLambda...\n")
sum_L2 = 0.0
for
i in range(0, n_eign):
sum_L2 += mm.sqr(Lmbd[i])
cum_sum = 0.0
self.f_out.write("%30s"%"Lambd ")
self.f_out.write("%10s"%"cum.sqr")
self.f_out.write("%10s\n"%"%")
print("\n")
print(("%30s"%"Lambda") +
("%10s"%"cum.sqr") + ("%10s"%"%") +
"\n")
for
i in range(0, n_eign):
self.f_out.write("%20s"%("Lambda-%d = "%(i + 1)))
self.f_out.write("%10.5f"%Lmbd[i])
cum_sum += mm.sqr(Lmbd[i])
self.f_out.write("%10.5f"%cum_sum)
self.f_out.write("%10.2f\n"%(100.0 * cum_sum / sum_L2))
print("%20s"%("Lambda-%d = "%(i + 1)),
end="")
print("%10.5f"%Lmbd[i], end = "")
print("%10.5f"%cum_sum, end = "")
print("%10.2f"%(100.0 * cum_sum / sum_L2))
s0 =
input("\nNumber of components = ")
n_comp = int(s0)
if
(n_comp > n_vars):
n_comp = n_vars
elif
(n_comp < 1):
n_comp = 1
if
(n_comp > n_eign):
n_comp = n_eign
Z =
ss.zscore(X, axis = 1).T
L1 =
numpy.diag(Lmbd[:n_comp])
V1 =
ListLV[:, :n_comp]
Pat
= V1 @ L1
self.f_out.write("\n\nPattern matrix for the principal components
=\n")
self.f_out.write("%10s"%" ")
for
i in range(0, n_comp):
self.f_out.write("%10s"%("comp.%d"%(i + 1)))
self.f_out.write("\n")
for
i in range(0, n_vars):
self.f_out.write("%10s"%var_names[i + 1])
for j in range(0, n_comp):
self.f_out.write("%10.5f"%Pat[i,
j])
self.f_out.write("\n")
PComp = numpy.matrix(Z @ V1 @ numpy.linalg.inv(L1))
return var_names, n_data, n_vars, CaseID, n_comp, Pat, PComp
Reference