高橋信『マンガ統計学因子分析編』の主成分分析（PCA）の素材を再分析する
factorのdataをtibbleで作る
MCAで得られた座標値で数量化する
- 総合順位を第一主成分負荷から見る
- 相関係数と固有値
  - 最適化したもの
  - 最適化前
eigen value 比較図示
- 表示させるデータの構成

実行日： 2018-03-03 15:52:52

本レポートは、Tokyo.R＃68（2018/03/03）でのLTの補足資料です。

library(FactoMineR)
library(factoextra)

## Loading required package: ggplot2

## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

library(explor)
library(tidyverse)

## Loading tidyverse: tibble
## Loading tidyverse: tidyr
## Loading tidyverse: readr
## Loading tidyverse: purrr
## Loading tidyverse: dplyr

## Conflicts with tidy packages ----------------------------------------------

## filter(): dplyr, stats
## lag():    dplyr, stats

高橋信『マンガ統計学因子分析編』の主成分分析（PCA）の素材を再分析する

ラーメン人気度分析 p105

データをmatrixで作る。

.d <- matrix(c(2,4,5,
               1,5,1,
               5,3,4,
               2,2,3,
               3,5,5,
               4,3,2,
               4,4,3,
               1,2,1,
               3,3,2,
               5,5,3),byrow=TRUE,10,3,
             dimnames = list(店名=c("二楽","夢田屋","地回","なのはな",
                                    "花ぶし","昇辰軒","丸蔵らあめん","海楽亭",
                                    "なるみ屋","奏月"),
                               要素=c("麺","具","スープ")))
.d

##               要素
## 店名         麺 具 スープ
##   二楽          2  4      5
##   夢田屋        1  5      1
##   地回          5  3      4
##   なのはな      2  2      3
##   花ぶし        3  5      5
##   昇辰軒        4  3      2
##   丸蔵らあめん  4  4      3
##   海楽亭        1  2      1
##   なるみ屋      3  3      2
##   奏月          5  5      3

res.PCA <- PCA(.d,graph=FALSE)
plot.PCA(res.PCA,axes = c(1,2),choix = "var")

plot.PCA(res.PCA,axes = c(1,3),choix = "var")

plot.PCA(res.PCA,axes = c(3,2),choix = "var")

plot.PCA(res.PCA,col.ind = 1:10)

第一主成分による順位

summary(res.PCA)

## 
## Call:
## PCA(X = .d, graph = FALSE) 
## 
## 
## Eigenvalues
##                        Dim.1   Dim.2   Dim.3
## Variance               1.573   0.814   0.613
## % of var.             52.428  27.134  20.438
## Cumulative % of var.  52.428  79.562 100.000
## 
## Individuals
##                  Dist    Dim.1    ctr   cos2    Dim.2    ctr   cos2  
## 二楽         |  1.721 |  0.750  3.581  0.190 |  0.550  3.714  0.102 |
## 夢田屋       |  2.343 | -1.027  6.702  0.192 |  1.993 48.817  0.724 |
## 地回         |  1.712 |  1.033  6.790  0.364 | -1.365 22.881  0.636 |
## なのはな     |  1.603 | -1.108  7.809  0.478 | -0.715  6.277  0.199 |
## 花ぶし       |  1.978 |  1.623 16.757  0.673 |  0.832  8.494  0.177 |
## 昇辰軒       |  1.104 | -0.292  0.541  0.070 | -0.784  7.547  0.504 |
## 丸蔵らあめん |  0.796 |  0.638  2.586  0.641 | -0.151  0.282  0.036 |
## 海楽亭       |  2.444 | -2.433 37.647  0.991 | -0.134  0.220  0.003 |
## なるみ屋     |  0.848 | -0.696  3.078  0.673 | -0.356  1.560  0.177 |
## 奏月         |  1.894 |  1.511 14.510  0.636 |  0.130  0.208  0.005 |
##               Dim.3    ctr   cos2  
## 二楽         -1.448 34.198  0.708 |
## 夢田屋        0.680  7.548  0.084 |
## 地回         -0.002  0.000  0.000 |
## なのはな     -0.911 13.547  0.323 |
## 花ぶし       -0.766  9.573  0.150 |
## 昇辰軒        0.721  8.473  0.426 |
## 丸蔵らあめん  0.452  3.339  0.323 |
## 海楽亭       -0.188  0.577  0.006 |
## なるみ屋      0.328  1.758  0.150 |
## 奏月          1.134 20.986  0.359 |
## 
## Variables
##           Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3    ctr
## 麺     |  0.717 32.662  0.514 | -0.545 36.539  0.297 |  0.435 30.799
## 具     |  0.655 27.261  0.429 |  0.712 62.348  0.508 |  0.252 10.392
## スープ |  0.794 40.077  0.630 | -0.095  1.114  0.009 | -0.600 58.809
##          cos2  
## 麺      0.189 |
## 具      0.064 |
## スープ  0.361 |

res.PCA$ind$coord[,1] %>% sort

##       海楽亭     なのはな       夢田屋     なるみ屋       昇辰軒 
##   -2.4333611   -1.1082693   -1.0267387   -0.6957621   -0.2916428 
## 丸蔵らあめん         二楽         地回         奏月       花ぶし 
##    0.6377176    0.7504515    1.0334490    1.5107110    1.6234449

単純加算による順位

.d %>% margin.table(1) %>% sort

## 店名
##       海楽亭       夢田屋     なのはな     なるみ屋       昇辰軒 
##            4            7            7            8            9 
##         二楽 丸蔵らあめん         地回       花ぶし         奏月 
##           11           11           12           13           13

「夢田屋」と「なのはな」、「二楽」と「丸蔵らあめん」、「花ぶし」と「奏月」の順位が入れ替わっている。

相関係数と固有値（寄与度、累積寄与度）

round(cor(.d),3)

##           麺    具 スープ
## 麺     1.000 0.191   0.36
## 具     0.191 1.000   0.30
## スープ 0.360 0.300   1.00

round(res.PCA$eig,2)

##        eigenvalue percentage of variance cumulative percentage of variance
## comp 1       1.57                  52.43                             52.43
## comp 2       0.81                  27.13                             79.56
## comp 3       0.61                  20.44                            100.00

res.PCA$eig %>% margin.table(2) %>% unlist %>% .[1:2]

##             eigenvalue percentage of variance 
##                      3                    100

factorのdataをtibbleで作る

MCAには、dfで投入しますが、変数はfactorでないとだめです。
data.frameにしておかないと、次のcolumn_to_rownames(‘店名’) で「それ、deprecated！」と叱られる。

tribble(
  ~店名,~麺,~具,~スープ,
  #---------------
  "二楽",2,4,5,
  "夢田屋",1,5,1,
  "地回",5,3,4,
  "なのはな",2,2,3,
  "花ぶし",3,5,5,
  "昇辰軒",4,3,2,
  "丸蔵らあめん",4,4,3,
  "海楽亭",1,2,1,
  "なるみ屋",3,3,2,
  "奏月",5,5,3
  ) %>%  mutate_all(.funs=factor) %>% as.data.frame  %>% column_to_rownames('店名') -> .d.f

多重対応分析（MCA）を行う

res.MCA  <- MCA(.d.f,graph=FALSE)
plot.MCA(res.MCA,axes = c(1,2),autoLab="yes",col.ind = rep(2,10),
         col.var = c(rep(3,5),rep(4,4),rep(5,5)),
         title="店-麺具スープ 1:2軸",cex=0.9)

固有値

もろもろ統計量を確認するには、summary(res.MCA) です。

res.MCA$eig

##       eigenvalue percentage of variance cumulative percentage of variance
## dim 1 0.84596940              23.071893                          23.07189
## dim 2 0.72209248              19.693431                          42.76532
## dim 3 0.59311220              16.175787                          58.94111
## dim 4 0.54922958              14.978988                          73.92010
## dim 5 0.39331799              10.726854                          84.64695
## dim 6 0.32511466               8.866763                          93.51372
## dim 7 0.11885650               3.241541                          96.75526
## dim 8 0.07845628               2.139717                          98.89498
## dim 9 0.04051758               1.105025                         100.00000

res.MCA$eig %>% margin.table(2) %>% unlist %>% .[1:2]

##             eigenvalue percentage of variance 
##               3.666667             100.000000

変数の「順序性」「等間隔性」を確認する

plot.MCA(res.MCA,axes = c(1,1),invisible = c("ind"),
         col.var = c(rep(3,5),rep(4,4),rep(5,5)))

plot.MCA(res.MCA,axes = c(2,2),invisible = c("ind"),
         col.var = c(rep(3,5),rep(4,4),rep(5,5)))

plot.MCA(res.MCA,axes = c(3,3),invisible = c("ind"),
         col.var = c(rep(3,5),rep(4,4),rep(5,5)))

クラスタの形成

res.MCA <- MCA(.d.f,graph=FALSE)
res.HCPC <- HCPC(res.MCA,nb.clust = 4,graph=FALSE)# nb.clust = xｘ を指定しないとエラー
plot.HCPC(res.HCPC)

plot.HCPC(res.HCPC,choice = "tree")

plot.HCPC(res.HCPC,choice = "map")

#plot.HCPC(res.HCPC,choice = "3D.map")

MCAで得られた座標値で数量化する

summary(res.MCA)

## 
## Call:
## MCA(X = .d.f, graph = FALSE) 
## 
## 
## Eigenvalues
##                        Dim.1   Dim.2   Dim.3   Dim.4   Dim.5   Dim.6
## Variance               0.846   0.722   0.593   0.549   0.393   0.325
## % of var.             23.072  19.693  16.176  14.979  10.727   8.867
## Cumulative % of var.  23.072  42.765  58.941  73.920  84.647  93.514
##                        Dim.7   Dim.8   Dim.9
## Variance               0.119   0.078   0.041
## % of var.              3.242   2.140   1.105
## Cumulative % of var.  96.755  98.895 100.000
## 
## Individuals
##                 Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3    ctr
## 二楽         | -0.190  0.426  0.009 | -1.566 33.983  0.613 | -0.162  0.443
## 夢田屋       | -1.278 19.299  0.474 |  0.971 13.068  0.274 | -0.321  1.736
## 地回         |  1.094 14.158  0.234 |  0.947 12.428  0.176 |  1.461 35.975
## なのはな     | -0.661  5.162  0.127 | -0.824  9.396  0.197 |  0.549  5.076
## 花ぶし       |  0.038  0.017  0.000 | -0.194  0.520  0.011 | -0.898 13.594
## 昇辰軒       |  1.185 16.612  0.408 |  0.265  0.975  0.020 | -0.618  6.450
## 丸蔵らあめん |  0.228  0.616  0.015 | -1.070 15.854  0.332 |  0.259  1.129
## 海楽亭       | -1.549 28.378  0.600 |  0.840  9.767  0.176 | -0.167  0.468
## なるみ屋     |  1.139 15.331  0.377 |  0.528  3.865  0.081 | -1.071 19.323
## 奏月         | -0.007  0.001  0.000 |  0.102  0.143  0.004 |  0.968 15.807
##                cos2  
## 二楽          0.007 |
## 夢田屋        0.030 |
## 地回          0.417 |
## なのはな      0.087 |
## 花ぶし        0.234 |
## 昇辰軒        0.111 |
## 丸蔵らあめん  0.019 |
## 海楽亭        0.007 |
## なるみ屋      0.333 |
## 奏月          0.325 |
## 
## Categories (the 10 first)
##             Dim.1    ctr   cos2 v.test    Dim.2    ctr   cos2 v.test  
## 麺_1     | -1.537 18.614  0.591 -2.305 |  1.066 10.486  0.284  1.599 |
## 麺_2     | -0.462  1.686  0.053 -0.694 | -1.406 18.261  0.494 -2.110 |
## 麺_3     |  0.640  3.227  0.102  0.960 |  0.197  0.358  0.010  0.295 |
## 麺_4     |  0.769  4.655  0.148  1.153 | -0.473  2.069  0.056 -0.710 |
## 麺_5     |  0.591  2.752  0.087  0.886 |  0.617  3.518  0.095  0.926 |
## 具_2     | -1.202 11.377  0.361 -1.802 |  0.009  0.001  0.000  0.014 |
## 具_3     |  1.239 18.146  0.658  2.433 |  0.683  6.459  0.200  1.341 |
## 具_4     |  0.021  0.003  0.000  0.031 | -1.551 22.218  0.602 -2.327 |
## 具_5     | -0.452  2.414  0.088 -0.887 |  0.345  1.648  0.051  0.677 |
## スープ_1 | -1.537 18.614  0.591 -2.305 |  1.066 10.486  0.284  1.599 |
##           Dim.3    ctr   cos2 v.test  
## 麺_1     -0.316  1.126  0.025 -0.475 |
## 麺_2      0.251  0.708  0.016  0.377 |
## 麺_3     -1.278 18.359  0.408 -1.917 |
## 麺_4     -0.234  0.613  0.014 -0.350 |
## 麺_5      1.577 27.953  0.622  2.365 |
## 具_2      0.248  0.692  0.015  0.372 |
## 具_3     -0.099  0.165  0.004 -0.194 |
## 具_4      0.063  0.044  0.001  0.094 |
## 具_5     -0.108  0.198  0.005 -0.213 |
## スープ_1 -0.316  1.126  0.025 -0.475 |
## 
## Categorical variables (eta2)
##          Dim.1 Dim.2 Dim.3  
## 麺     | 0.785 0.752 0.868 |
## 具     | 0.811 0.657 0.020 |
## スープ | 0.942 0.758 0.892 |

麺t <- res.MCA$var$coord[1:5,1]
具t <- res.MCA$var$coord[6:9,1]
スープt <- res.MCA$var$coord[10:14,1]

.d2 <- data.frame(麺t[.d.f$麺],具t[.d.f$具],スープt[.d.f$スープ])
rownames(.d2) <- c("二楽","夢田屋","地回","なのはな",
　　　　　　　　"花ぶし","昇辰軒", "丸蔵らあめん",
　　　　　　　　"海楽亭","なるみ屋","奏月")
.d2

##              麺t..d.f.麺. 具t..d.f.具. スープt..d.f.スープ.
## 二楽           -0.4624864   0.02087463          -0.08244833
## 夢田屋         -1.5368811  -0.45187617          -1.53688111
## 地回            0.5909102   1.23897692           1.18985460
## なのはな       -0.4624864  -1.20152575          -0.15942396
## 花ぶし          0.6398921  -0.45187617          -0.08244833
## 昇辰軒          0.7685653   1.23897692           1.26353808
## 丸蔵らあめん    0.7685653   0.02087463          -0.15942396
## 海楽亭         -1.5368811  -1.20152575          -1.53688111
## なるみ屋        0.6398921   1.23897692           1.26353808
## 奏月            0.5909102  -0.45187617          -0.15942396

res.PCA2 <- PCA(.d2,graph=FALSE)
plot.PCA(res.PCA2,axes = c(1,2),choix = "var")

plot.PCA(res.PCA2,axes = c(1,3),choix = "var")

plot.PCA(res.PCA2,axes = c(3,2),choix = "var")

plot.PCA(res.PCA2,choix = "ind",col.ind = 1:10,title="最適化尺度化後のPCA")

plot.PCA(res.PCA,choix = "ind",col.ind = 1:10,title="リッカート尺度でのPCA")

総合順位を第一主成分負荷から見る

res.PCA2$ind$coord[,1] %>% sort

##       海楽亭       夢田屋     なのはな         二楽         奏月 
##  -2.68365954  -2.21309375  -1.14459701  -0.32896007  -0.01279903 
##       花ぶし 丸蔵らあめん         地回     なるみ屋       昇辰軒 
##   0.06626627   0.39547004   1.89553518   1.97253395   2.05330397

相関係数と固有値

最適化したものでは、第一軸が84.6%、第2軸まであわせると97.0% の表現。リッカートスケールのものは、第一軸が、52.4%、第2軸までで、79.6％であることがわかる。

最適化したもの

round(cor(.d2),3)

##                      麺t..d.f.麺. 具t..d.f.具. スープt..d.f.スープ.
## 麺t..d.f.麺.                 1.00        0.630                0.820
## 具t..d.f.具.                 0.63        1.000                0.852
## スープt..d.f.スープ.         0.82        0.852                1.000

res.PCA2$eig

##        eigenvalue percentage of variance cumulative percentage of variance
## comp 1 2.53790821              84.596940                          84.59694
## comp 2 0.37124514              12.374838                          96.97178
## comp 3 0.09084665               3.028222                         100.00000

res.PCA2$eig %>% margin.table(2)

##                        eigenvalue            percentage of variance 
##                            3.0000                          100.0000 
## cumulative percentage of variance 
##                          281.5687

#fviz_eig(res.PCA2)

最適化前

round(cor(.d),3)

##           麺    具 スープ
## 麺     1.000 0.191   0.36
## 具     0.191 1.000   0.30
## スープ 0.360 0.300   1.00

res.PCA$eig

##        eigenvalue percentage of variance cumulative percentage of variance
## comp 1  1.5728539               52.42846                          52.42846
## comp 2  0.8140083               27.13361                          79.56207
## comp 3  0.6131378               20.43793                         100.00000

res.PCA$eig %>% margin.table(2)

##                        eigenvalue            percentage of variance 
##                            3.0000                          100.0000 
## cumulative percentage of variance 
##                          231.9905

#fviz_eig(res.PCA)

eigen value 比較図示

表示させるデータの構成

res.PCA$eig %>% .[,2:3] -> .dd1
res.PCA2$eig %>%  .[,2:3] -> .dd2
colnames(.dd1) <- c("固有値","累積％")
colnames(.dd2) <- c("固有値","累積％")

#.dd1 %>% as.data.frame %>%  rownames_to_column("軸") %>% mutate(区分="整数尺度") -> .dd1.a
#.dd2 %>% as.data.frame %>% rownames_to_column("軸")  %>% mutate(区分="最適尺度") -> .dd2.a
#rbind(.dd1.a,.dd2.a) %>% gather(key = 種類, -軸,-区分,value=値) -> .data

.dd1 %>% as.data.frame %>%  rownames_to_column("軸") %>% mutate(区分="整数尺度") %>%
  bind_rows(.dd2 %>% as.data.frame %>% rownames_to_column("軸")  %>% 
              mutate(区分="最適尺度")) %>% gather(key = 種類, -軸,-区分,value=値) -> .data

.data %>% filter(種類=="固有値") %>% 
  ggplot(aes(x=軸,y=値)) + geom_bar(stat="identity",position = "dodge",fill="lightblue") -> p
p <- p + geom_line(aes(x=軸,y=値,group=種類)) + geom_point()
p <- p + xlab("主軸") + ylab("固有値％")
p <- p + facet_grid(. ~ factor(区分,levels = c("整数尺度","最適尺度"))) + 
  theme_bw(base_family = "sans") + theme(strip.background = element_rect(colour="black", fill="#CCCCFF"))
p <- p + geom_text(aes(label=round(値,2),vjust=-0.3))
p

その数量化、大丈夫ですか？整数尺度と最適化尺度

kazuo.fujimoto2007@gmail.com

2018年2月28日

高橋信『マンガ統計学因子分析編』の主成分分析（PCA）の素材を再分析する

ラーメン人気度分析 p105

第一主成分による順位

単純加算による順位

相関係数と固有値（寄与度、累積寄与度）

factorのdataをtibbleで作る

多重対応分析（MCA）を行う

固有値

変数の「順序性」「等間隔性」を確認する

クラスタの形成

MCAで得られた座標値で数量化する

総合順位を第一主成分負荷から見る

相関係数と固有値

最適化したもの

最適化前

eigen value 比較図示

表示させるデータの構成

その数量化、大丈夫ですか？ 整数尺度と最適化尺度

kazuo.fujimoto2007@gmail.com

2018年2月28日

高橋信『マンガ統計学 因子分析編』の主成分分析（PCA）の素材を再分析する

ラーメン人気度分析 p105

第一主成分による順位

単純加算による順位

相関係数と固有値 （寄与度、累積寄与度）

factorのdataをtibbleで作る

多重対応分析（MCA）を行う

固有値

変数の「順序性」「等間隔性」を確認する

クラスタの形成

MCAで得られた座標値で数量化する

総合順位を第一主成分負荷から見る

相関係数と固有値

最適化したもの

最適化前

eigen value 比較図示

表示させるデータの構成

その数量化、大丈夫ですか？整数尺度と最適化尺度

高橋信『マンガ統計学因子分析編』の主成分分析（PCA）の素材を再分析する

相関係数と固有値（寄与度、累積寄与度）