因子分析法

  1. 因子分析の考え方
  2. 因子負荷量の求め方
  3. 因子得点の求め方
  4. 因子分析のプログラム

1.因子分析の考え方

トップに戻る


              

2.因子負荷量の求め方

トップに戻る


3.因子得点の求め方

トップに戻る


4.因子分析のプログラム

(1)因子負荷量を求めるプログラムfactor1.R
factor1<-function(r){
n<-ncol(r)
evalue<-eigen(r)$values
evector<-eigen(r)$vectors
evalue1<-evalue[round(evalue,6)>0]
nc<-length(evalue1)
cont.evalue1<-evalue1/n
d1<-matrix(0,ncol=nc,nrow=nc)
d1[col(d1)>=row(d1)]<-1
cum.evalue<-cont.evalue1%*%d1
lamda<-diag(evalue1,nc)
evector<-evector[,1:nc]
b<-apply(evector^2,2,sum)
for ( j in 1:nc) evector[,j]<-evector[,j]/sqrt(b[j])
a<-evector%*%sqrt(lamda)
list(r=round(r,3),evalue=round(evalue1,3),
     cum=round(cum.evalue,3),loading=round(a,3),
     evector=round(evector,3))
}

実行例

x<-matrix(c(5,4,3,2,1,4,4,3,4,3,
            5,5,4,4,3,5,5,3,4,3,
            3,3,2,2,1,3,3,2,2,2,
            4,3,4,2,2,2,3,2,3,1,
            5,5,4,3,2,2,3,2,2,1),ncol=10,byrow=T)
out<-factor1(cor(x))
print(out)

出力 

$r
       [,1]  [,2]   [,3]  [,4]  [,5]   [,6]  [,7]  [,8]  [,9]  [,10]
 [1,] 1.000 0.839  0.688 0.562 0.468  0.343 0.562 0.612 0.559  0.280
 [2,] 0.839 1.000  0.559 0.839 0.598  0.383 0.559 0.456 0.250  0.250
 [3,] 0.688 0.559  1.000 0.562 0.802 -0.086 0.250 0.102 0.280 -0.280
 [4,] 0.562 0.839  0.562 1.000 0.869  0.514 0.688 0.408 0.280  0.280
 [5,] 0.468 0.598  0.802 0.869 1.000  0.275 0.535 0.218 0.299  0.000
 [6,] 0.343 0.383 -0.086 0.514 0.275  1.000 0.943 0.910 0.767  0.959
 [7,] 0.562 0.559  0.250 0.688 0.535  0.943 1.000 0.919 0.839  0.839
 [8,] 0.612 0.456  0.102 0.408 0.218  0.910 0.919 1.000 0.913  0.913
 [9,] 0.559 0.250  0.280 0.280 0.299  0.767 0.839 0.913 1.000  0.750
[10,] 0.280 0.250 -0.280 0.280 0.000  0.959 0.839 0.913 0.750  1.000

$evalue
[1] 5.861 2.671 0.816 0.651

$cum
      [,1]  [,2]  [,3] [,4]
[1,] 0.586 0.853 0.935    1

$loading
       [,1]   [,2]   [,3]   [,4]
 [1,] 0.753 -0.368 -0.454  0.303
 [2,] 0.727 -0.478  0.071  0.487
 [3,] 0.438 -0.807 -0.304 -0.254
 [4,] 0.762 -0.469  0.443  0.059
 [5,] 0.614 -0.647  0.268 -0.364
 [6,] 0.843  0.486  0.230 -0.028
 [7,] 0.965  0.203  0.122 -0.112
 [8,] 0.892  0.413 -0.179  0.037
 [9,] 0.809  0.331 -0.368 -0.316
[10,] 0.724  0.678  0.092  0.084

$evector
       [,1]   [,2]   [,3]   [,4]
 [1,] 0.311 -0.225 -0.503  0.376
 [2,] 0.300 -0.293  0.079  0.604
 [3,] 0.181 -0.494 -0.337 -0.315
 [4,] 0.315 -0.287  0.490  0.073
 [5,] 0.253 -0.396  0.297 -0.451
 [6,] 0.348  0.297  0.255 -0.035
 [7,] 0.399  0.124  0.135 -0.138
 [8,] 0.368  0.253 -0.199  0.046
 [9,] 0.334  0.202 -0.407 -0.392
[10,] 0.299  0.415  0.102  0.104


      
(2)因子得点を求めるプログラムfscore.R
fscore<-function(x,a){
   n<-ncol(x)
   m<-ncol(a)
   z<-scale(x,center=T,scale=T)
   a2<-t(a)%*%a
   f<-z%*%a%*%solve(a2)
   list(rawdata=round(x,3),loading=round(a,3),fscore=round(f,3))
}
実行例
x<-matrix(c(5,4,3,2,1,4,4,3,4,3,
            5,5,4,4,3,5,5,3,4,3,
            3,3,2,2,1,3,3,2,2,2,
            4,3,4,2,2,2,3,2,3,1,
            5,5,4,3,2,2,3,2,2,1),ncol=10,byrow=T)
a<-factor1(cor(x))$loading
out<-fscore(x,a)
print(out)

出力
$rawdata
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    5    4    3    2    1    4    4    3    4     3
[2,]    5    5    4    4    3    5    5    3    4     3
[3,]    3    3    2    2    1    3    3    2    2     2
[4,]    4    3    4    2    2    2    3    2    3     1
[5,]    5    5    4    3    2    2    3    2    2     1

$loading
       [,1]   [,2]   [,3]   [,4]
 [1,] 0.753 -0.368 -0.454  0.303
 [2,] 0.727 -0.478  0.071  0.487
 [3,] 0.438 -0.807 -0.304 -0.254
 [4,] 0.762 -0.469  0.443  0.059
 [5,] 0.614 -0.647  0.268 -0.364
 [6,] 0.843  0.486  0.230 -0.028
 [7,] 0.965  0.203  0.122 -0.112
 [8,] 0.892  0.413 -0.179  0.037
 [9,] 0.809  0.331 -0.368 -0.316
[10,] 0.724  0.678  0.092  0.084

$fscore
       [,1]   [,2]   [,3]   [,4]
[1,]  0.456  1.085 -1.225  0.564
[2,]  1.499 -0.179  0.831 -0.482
[3,] -1.011  0.901  1.156  0.172
[4,] -0.671 -0.487 -0.686 -1.429
[5,] -0.272 -1.320 -0.076  1.174
 
(3)varimax回転のプログラムvarimax2
varimax2<-function(am)
{
m <- ncol(am)
n <- nrow(am)
commu <- apply(am^2,1,sum)
icomm <- 1/sqrt(commu)
h <- diag(icomm,nrow=length(commu))
aini <- h%*%am
b <- aini
m1 <- m-1
kk <- 1000
for ( k in 1:kk){b2 <- b;
for ( p in 1:m1) {p1<- p+1;for ( q in p1:m){ 
                bb <- b^2; bpq1 <-  bb[,p] - bb[,q];
                bpq2 <- b[,p]*b[,q];
                a1 <- sum(bpq1); b1 <- 2*sum(bpq2);
                c1 <- sum(bpq1^2)-4*sum(bpq2^2);
                d1 <- 4*sum(bpq1*bpq2);
                v1 <- d1 - 2*a1*b1/n;
                v2 <- c1 - (a1^2-b1^2)/n;
                if ( v1 > 0 && v2 > 0) theta4 <- atan(v1/v2);
                if ( v1 > 0 && v2 < 0) theta4 <- atan(v1/v2) + 3.14;
                if ( v1 < 0 && v2 > 0) theta4 <- atan(v1/v2);
                if ( v1 < 0 && v2 < 0) theta4 <- atan(v1/v2) - 3.14;
                theta <- theta4/4;
                costhe <- cos(theta);
                sinthe <- sin(theta);
                t1 <- diag(m);
                t1[p,p] <- costhe;
                t1[p,q] <- -sinthe;
                t1[q,q] <- costhe;
                t1[q,p] <- sinthe;
                b <- b%*%t1;}}
                bidf<- (b2 - b)^2;
                dif <- sum(bidf)/(m*n);
                if(dif < 0.001) {
                       kk<- k;
                       break}}
                b<-solve(h)%*%b ;
                eigen1<- apply(b^2,2,sum);
                list(data=round(am,3),commu=round(commu,3),
                     loading=round(b,3),evalue=round(eigen1,2),iter=kk)
}


(4)主成分分析との比較
 営業社員のデータを標準得点に変換して後、主成分 分析を行う
prcomp2(scale(x))
$data
      [,1] [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8] [,9] [,10]
[1,]  0.67    0 -0.45 -0.67 -0.96  0.61  0.45  1.10    1     1
[2,]  0.67    1  0.67  1.57  1.43  1.38  1.57  1.10    1     1
[3,] -1.57   -1 -1.57 -0.67 -0.96 -0.15 -0.67 -0.73   -1     0
[4,] -0.45   -1  0.67 -0.67  0.24 -0.92 -0.67 -0.73    0    -1
[5,]  0.67    1  0.67  0.45  0.24 -0.92 -0.67 -0.73   -1    -1

$cov
      [,1] [,2]  [,3] [,4] [,5]  [,6] [,7] [,8] [,9] [,10]
 [1,] 1.00 0.84  0.69 0.56 0.47  0.34 0.56 0.61 0.56  0.28
 [2,] 0.84 1.00  0.56 0.84 0.60  0.38 0.56 0.46 0.25  0.25
 [3,] 0.69 0.56  1.00 0.56 0.80 -0.09 0.25 0.10 0.28 -0.28
 [4,] 0.56 0.84  0.56 1.00 0.87  0.51 0.69 0.41 0.28  0.28
 [5,] 0.47 0.60  0.80 0.87 1.00  0.28 0.53 0.22 0.30  0.00
 [6,] 0.34 0.38 -0.09 0.51 0.28  1.00 0.94 0.91 0.77  0.96
 [7,] 0.56 0.56  0.25 0.69 0.53  0.94 1.00 0.92 0.84  0.84
 [8,] 0.61 0.46  0.10 0.41 0.22  0.91 0.92 1.00 0.91  0.91
 [9,] 0.56 0.25  0.28 0.28 0.30  0.77 0.84 0.91 1.00  0.75
[10,] 0.28 0.25 -0.28 0.28 0.00  0.96 0.84 0.91 0.75  1.00

$cor
      [,1] [,2]  [,3] [,4] [,5]  [,6] [,7] [,8] [,9] [,10]
 [1,] 1.00 0.84  0.69 0.56 0.47  0.34 0.56 0.61 0.56  0.28
 [2,] 0.84 1.00  0.56 0.84 0.60  0.38 0.56 0.46 0.25  0.25
 [3,] 0.69 0.56  1.00 0.56 0.80 -0.09 0.25 0.10 0.28 -0.28
 [4,] 0.56 0.84  0.56 1.00 0.87  0.51 0.69 0.41 0.28  0.28
 [5,] 0.47 0.60  0.80 0.87 1.00  0.28 0.53 0.22 0.30  0.00
 [6,] 0.34 0.38 -0.09 0.51 0.28  1.00 0.94 0.91 0.77  0.96
 [7,] 0.56 0.56  0.25 0.69 0.53  0.94 1.00 0.92 0.84  0.84
 [8,] 0.61 0.46  0.10 0.41 0.22  0.91 0.92 1.00 0.91  0.91
 [9,] 0.56 0.25  0.28 0.28 0.30  0.77 0.84 0.91 1.00  0.75
[10,] 0.28 0.25 -0.28 0.28 0.00  0.96 0.84 0.91 0.75  1.00

$evalue
 [1] 5.86 2.67 0.82 0.65 0.00 0.00 0.00 0.00 0.00 0.00

$cum.cont
 [1] 0.59 0.85 0.93 1.00 1.00 1.00 1.00 1.00 1.00 1.00

$evector
      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10]
 [1,] 0.31 -0.22 -0.50  0.38 -0.43  0.00  0.05 -0.42  0.20 -0.24
 [2,] 0.30 -0.29  0.08  0.60  0.27  0.06  0.00  0.37  0.17  0.46
 [3,] 0.18 -0.49 -0.34 -0.31  0.02  0.39 -0.16  0.29 -0.50 -0.07
 [4,] 0.31 -0.29  0.49  0.07 -0.15 -0.49  0.25 -0.07 -0.46 -0.19
 [5,] 0.25 -0.40  0.30 -0.45  0.03  0.08  0.04  0.03  0.67 -0.19
 [6,] 0.35  0.30  0.25 -0.03 -0.12  0.64  0.44 -0.23 -0.13  0.18
 [7,] 0.40  0.12  0.14 -0.14  0.05 -0.07 -0.72 -0.40 -0.08  0.31
 [8,] 0.37  0.25 -0.20  0.05  0.72 -0.06  0.09 -0.06 -0.03 -0.47
 [9,] 0.33  0.20 -0.41 -0.39 -0.10 -0.42  0.33  0.19  0.06  0.44
[10,] 0.30  0.41  0.10  0.10 -0.41  0.04 -0.27  0.59  0.06 -0.35

$norm.evector
 [1] 1 1 1 1 1 1 1 1 1 1

$prcompscore
      [,1]  [,2]  [,3]  [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]  1.10  1.77 -1.11  0.46    0    0    0    0    0     0
[2,]  3.63 -0.29  0.75 -0.39    0    0    0    0    0     0
[3,] -2.45  1.47  1.04  0.14    0    0    0    0    0     0
[4,] -1.62 -0.80 -0.62 -1.15    0    0    0    0    0     0
[5,] -0.66 -2.16 -0.07  0.95    0    0    0    0    0     0

$loading
      [,1]  [,2]  [,3]  [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,] 0.75 -0.37 -0.45  0.30    0    0  NaN  NaN  NaN   NaN
 [2,] 0.73 -0.48  0.07  0.49    0    0  NaN  NaN  NaN   NaN
 [3,] 0.44 -0.81 -0.30 -0.25    0    0  NaN  NaN  NaN   NaN
 [4,] 0.76 -0.47  0.44  0.06    0    0  NaN  NaN  NaN   NaN
 [5,] 0.61 -0.65  0.27 -0.36    0    0  NaN  NaN  NaN   NaN
 [6,] 0.84  0.49  0.23 -0.03    0    0  NaN  NaN  NaN   NaN
 [7,] 0.97  0.20  0.12 -0.11    0    0  NaN  NaN  NaN   NaN
 [8,] 0.89  0.41 -0.18  0.04    0    0  NaN  NaN  NaN   NaN
 [9,] 0.81  0.33 -0.37 -0.32    0    0  NaN  NaN  NaN   NaN
[10,] 0.72  0.68  0.09  0.08    0    0  NaN  NaN  NaN   NaN

$evector2
      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10]
 [1,] 0.31 -0.22 -0.50  0.38 -0.43  0.00  0.05 -0.42  0.20 -0.24
 [2,] 0.30 -0.29  0.08  0.60  0.27  0.06  0.00  0.37  0.17  0.46
 [3,] 0.18 -0.49 -0.34 -0.31  0.02  0.39 -0.16  0.29 -0.50 -0.07
 [4,] 0.31 -0.29  0.49  0.07 -0.15 -0.49  0.25 -0.07 -0.46 -0.19
 [5,] 0.25 -0.40  0.30 -0.45  0.03  0.08  0.04  0.03  0.67 -0.19
 [6,] 0.35  0.30  0.25 -0.03 -0.12  0.64  0.44 -0.23 -0.13  0.18
 [7,] 0.40  0.12  0.14 -0.14  0.05 -0.07 -0.72 -0.40 -0.08  0.31
 [8,] 0.37  0.25 -0.20  0.05  0.72 -0.06  0.09 -0.06 -0.03 -0.47
 [9,] 0.33  0.20 -0.41 -0.39 -0.10 -0.42  0.33  0.19  0.06  0.44
[10,] 0.30  0.41  0.10  0.10 -0.41  0.04 -0.27  0.59  0.06 -0.35

$prcompscore2
      [,1]  [,2]  [,3]  [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]  1.10  1.77 -1.11  0.46    0    0    0    0    0     0
[2,]  3.63 -0.29  0.75 -0.39    0    0    0    0    0     0
[3,] -2.45  1.47  1.04  0.14    0    0    0    0    0     0
[4,] -1.62 -0.80 -0.62 -1.15    0    0    0    0    0     0
[5,] -0.66 -2.16 -0.07  0.95    0    0    0    0    0     0

トップに戻る