主成分分析

伊藤暢浩研究室 統計ゼミ 2020 主成分分析 せつめいする人 安藤 圭祐

はじめに

3 講義の目標 主成分分析の基本的なしくみを理解する Excelを使った主成分分析に関する行列計算を身に着ける Rを使った主成分分析の方法とその結果の見方を学習する 主成分分析を自分の研究で利用可能なツールの1つとして獲得する

4 講義で使用するファイル 202009 5_4 主成分分析.pdf この講義資料 202009 5_4 主成分分析(演習ブック).xlsx 中編Excel演習で使用 202009 5_4 主成分分析(データ) tsutaya_rental.csv companies.csv 後編Rでの演習で使用 aichi_offices.csv SSDSE-2020A.xlsx 202009 5_4 主成分分析(作成用) 参考データ 主成分分析.pptx HackGen_v2.1.1.zip 資料作成時のデータ

主成分分析でできること

6 分析するデータの例 事業所数 農林 建設 製造業 電気・ガス 情報通信 … 名古屋市 58 8,654 10,896 116 2413 … 豊橋市 104 1,515 1,694 22 128 … 豊田市 74 1,434 1,564 21 101 … … 指標(次元)の数が多く扱いにくい

7 分析した結果の例 豊橋市 田原市 豊田市 稲沢市 岡崎市 名古屋市 都市圏で事業所が多い地域 一宮市 農 業 が 盛 ん な 地 域

8 主成分分析 principal component analysis; PCA 次元の高いデータ 次元の低いデータ 可視化 (データの要約) データの特徴をおさえた扱いやすい形に変換することができる

主成分分析の考え方

10 愛知県の市の男女別総人口 総人口 男性 女性 25 20 15 豊橋市 187,801 186,964 10 瀬戸市 63,189 65,857 5 豊田市 222,169 200,373 0 安城市 94,073 90,067 豊橋市 瀬戸市 男性 豊田市 女性 安城市

11 2つの変量の合成 総人口 豊橋市 次元削減 男性 女性 男女の平均人口 187,801 186,964 187,383 𝑥1 𝑥2 0.5𝑥1 + 0.5𝑥2 Dimentionality Reduction 多次元データの次元数を減らすこと 有用な情報が残るように次元削減をおこなうことでデータの特徴を捉えることができる

12 グラフで見る次元削減 女性(𝑥2 )／万人 25 豊田市 豊橋市 20 豊田市 瀬戸市 安城市 豊橋市 15 10 安城市 瀬戸市 6.5 5 𝑂 9.2 18.7 21.1 男女の平均人口(0.5𝑥1 + 0.5𝑥2 ) ／万人 5 10 15 20 25 男性(𝑥1 ) ／万人

13 ベクトルで考える情報の損失 男女の平均人口(軸) (次元削減した)豊橋市 情報の損失 (軸との距離) ＝ 豊橋市の男女の平均人口 豊橋市 = + = 豊橋市 豊橋市の女性の 総人口 豊橋市の男性の総人口 豊橋市の女性の総人口 豊橋市の男女の平均人口 + 情報の損失 男女の総人口の間の 差(偏り)の情報が失われた 豊橋市の男性の総人口

14 合成変量を再定義して情報の損失を軽減 0.500𝑥1 + 0.500𝑥2 0.500𝑥1 + 0.499𝑥2 元の合成変量(男女の平均人口) 𝑓: 0.500𝑥1 + 0.500𝑥2 情報の損失 再定義した合成変量 𝑓: 0.500𝑥1 + 0.499𝑥2 豊橋市 18.7,18.6 一般化した合成変量 𝑓: 𝑤1 𝑥1 + 𝑤2 𝑥2 = 𝑥1 𝑥2 係数ベクトル 𝑤1 𝑤2 を調整することで，情報の損失を軽減できる 𝑤1 𝑤2

15 合成変量を増やして情報の損失を軽減 第2合成変量 第1合成変量(男女の平均人口) 豊橋市の男女の平均人口 第1合成変量(男女の平均人口) 情報の損失 豊橋市 𝑓: 0.500𝑥1 + 0.500𝑥2 第2合成変量 𝑔: −0.900𝑥1 + 0.900𝑥2 異なる合成変量を追加することで情報の損失を軽減できる

16 情報の損失をゼロにする合成変量 女性人口に対する 男性人口の多さ 男女の平均人口 豊橋市の女性人口に対する 男性人口の多さ 豊橋市の男女の平均人口 情報の損失= 0 第1合成変量(男女の平均人口) 𝑓 = 0.500𝑥1 + 0.500𝑥2 豊橋市 第2合成変量(男女の人口の差) 𝑔 = 0.500𝑥1 − 0.500𝑥2 情報の損失がゼロとなるような合成変量は他の合成変量に直交する

17 主成分分析の考え方 有用な情報を残して次元を削減 目的 複数の変量を1つの合成変量に統合 情報の損失が発生 情報の損失の軽減 係数の調整 合成変量の追加 次元数増化 なるべく少ない数の「情報の損失が最小となるような係数を持つ合成変量」でデータを表現

主成分分析のしくみ

19 主成分分析の手順 1 情報の損失が最小となるような 係数を持つ合成変量を計算し定義 第2主成分得点 2 1で吸収できていない情報に対して その損失が最も小さくなる 直交する軸を計算し定義 第2主成分 第1主成分 3 1・2を繰り返して𝑛次元のデータに 対して𝑛本の軸を見つけたら終了 第1主成分得点

20 情報の損失と主成分得点の分散の関係 情報の損失が大きい 主成分得点の分散が小さい 情報の損失の最小化 情報の損失が小さい 主成分得点の分散が大きい 主成分得点の分散の最大化

21 主成分得点の分散の求め方 データ 𝑖 = 1, 2, … , 𝑛 各主成分のスケールを合わせた係数べクトル 指標(次元) 𝑗 = 1, 2, … , 𝑝 厳密には𝑥𝑖𝑗 − 𝜇𝑗 主成分に平行な単位ベクトル 𝕨 = 𝑤1 𝑤2 ⋯ 𝑤𝑝 𝑇 𝕩1 各データへのベクトル 𝕩𝑛 = 𝑥𝑛1 ⋯ 𝑥𝑛𝑝 𝑥𝑛2 𝕩2 主成分得点 𝕩𝑛 ∙ 𝕨 主成分得点の分散 𝕩3 𝑛 1 𝕩𝑖 ∙ 𝕨 𝑛 − 1 (主成分得点の平均は0) 𝑖=1 𝕨 2 主成分 𝕩1 ∙ 𝕨 (= 𝑥11 𝑤1 + 𝑥12 𝑤2 + ⋯ + 𝑥1𝑝 𝑤𝑝 )

22 「主成分得点の分散」の式の整理 1 𝑛−1 𝑛 𝑖=1 𝕩1 ∙ 𝕨 𝕩2 ∙ 𝕨 1 2 𝕩𝑛 ∙ 𝕨 = ⋮ 𝑛−1 𝕩𝑛 ∙ 𝕨 𝑇 𝕩1 𝕩2 ⋮ 𝕩𝑛 𝑇 1 = 𝕨𝑇 𝑛−1 1 = 𝕨𝑇 𝕏𝑇 𝕏𝕨 𝑛−1 𝕩1 ∙ 𝕨 𝕩2 ∙ 𝕨 ⋮ 𝕩𝑛 ∙ 𝕨 例．データ数𝑛 = 2，次元数𝑝 = 2の場合 𝕩1 ∙ 𝕨 𝕩2 ∙ 𝕨 𝕩1 𝕩2 ⋮ 𝕨 𝕩𝑛 𝕩1 𝕩2 ⋮ 𝕩𝑛 𝕏≡ 1 𝑥 𝑤 +𝑥 𝑤 = 𝑥11 𝑤1 + 𝑥12 𝑤2 21 1 22 2 𝑥11 𝑥12 𝑤1 = 𝑥 𝑤2 21 𝑥22 𝕩1 = 𝕩 𝕨 2 = 𝑥11 𝑥21 ⋯ 𝑥1𝑝 ⋯ 𝑥2𝑝 𝑥12 𝑥22 ⋮ 𝑥𝑛1 𝑥𝑛2 ⋯ 𝑥𝑛𝑝 データ表の行列𝕏の転置行列と𝕏の内積を含む𝑛−1 𝕏𝑇 𝕏という式が現れる

23 分散共分散行列 variance-covariance matrix Σ= 𝜎 2 (𝑥1 ) 𝐶(𝑥1 , 𝑥2 ) ⋯ 𝐶(𝑥1 , 𝑥𝑝 ) 𝐶(𝑥1 , 𝑥2 ) 𝜎 2 (𝑥2 ) ⋯ 𝐶(𝑥2 , 𝑥𝑝 ) ⋮ 𝐶(𝑥1 , 𝑥𝑝 ) 𝐶(𝑥2 , 𝑥𝑝 ) ⋯ 𝜎 2 (𝑥𝑝 ) 分散𝜎 𝑥𝑖 共分散𝐶 𝑥𝑖 , 𝑥𝑗 分散 𝑖 = 𝑗のとき𝜎 𝑥𝑖 = 𝐶 𝑥𝑖 , 𝑥𝑗 であるため 共分散行列とも呼ばれる 共分散 variance データの散らばり具合を表す指標 1 2 𝜎 𝑥 = 𝑛−1 𝐶 𝑥𝑖 , 𝑥𝑗 = 𝐶 𝑥𝑗 , 𝑥𝑖 より，対称行列 𝑛 𝑥𝑖 − 𝜇 2 𝑖=1 covariance 2変量𝑥，𝑦間の関係を表す指標の1つ 1 𝐶 𝑥, 𝑦 = 𝑛−1 𝑛 𝑥𝑖 − 𝜇 𝑦𝑖 − 𝜇 𝑖=1 ※ここでの𝑖， 𝑗は一般的な表現に基づきそれぞれ「何番目の行」，「何番目の列」という意味で使っています． 記号が足りないので許してください．

24 主成分得点の分散と分散共分散行列の関係 1 𝑇 𝕏 𝕏 𝑛 対角成分が𝜎(𝑥𝑖 ) 非対角成分が𝐶(𝑥𝑖 , 𝑥𝑗 ) 主成分得点の分散は分散共分散行列Σを使って表すことができる

25 「主成分得点の分散」を最大化するベクトル 1 𝑛−1 𝑛 𝕩𝑛 ∙ 𝕨 𝑖=1 2 1 = 𝕨𝑇 𝕏𝑇 𝕏𝕨 = 𝕨𝑇 Σ𝕨 𝑛−1 最大化する𝕨を求めたい 「𝕨が単位ベクトルである」という制約下で ラグランジュの未定乗数法(割愛) 固有方程式 Σ𝕨 = 𝜆𝕨 を満たす固有値と固有ベクトルの組(𝜆, 𝕨)を求める 行列式を解く必要がある べき乗法をはじめとする数値解析手法によって解く 解が複数となるが固有値が大きいものから第𝑗主成分の 分散𝜆と平行なベクトル𝕨となる すべての主成分は分散共分散行列Σから解析的に求めることができる

Excelをつかって 簡単な主成分分析

27 演習内容 愛知県の4つ市における男女別の人口のデータを使って 次の2つを求めてみよう． 1 分散共分散行列 𝜎 2 (𝑥1 ) 𝐶(𝑥1 , 𝑥2 ) 𝐶(𝑥1 , 𝑥2 ) 𝜎 2 (𝑥2 ) 2 第1主成分 20 10 𝑂 10 20

28 使用するデータ 愛知県の男女別の人口(万人) 男性 女性 豊橋市 18.8 18.7 瀬戸市 6.3 6.6 豊田市 22.2 20.0 安城市 9.4 9.0 ExcelのC4:D7にすでに記載

29 演習の手順 1 各項目の平均値を計算 2 データと平均値を使って分散共分散行列を計算 3 べき乗法で固有ベクトルを計算

30 べき乗法 分散共分散行列 𝜎 2 (𝑥1 ) 𝐶(𝑥1 , 𝑥2 ) 8.5 6.3 Σ= = 6.3 11.7 𝐶(𝑥1 , 𝑥2 ) 𝜎 2 (𝑥2 ) Σ∙𝕧 1 −1 Σ∙𝕧 2.2 −5.4 Σ∙𝕧 −15.3 −49.3 第1主成分 Σ∙𝕧 −440.9 −673.5 −7991.4 −10658.5 𝕧を更新していくと第1主成分と 平行なベクトルに限りなく近づく 計算するときは𝕧の大きさが 大きくなっていくため 𝕧 で割る 初期ベクトル

Rをつかって 主成分分析

32 主成分の計算方法 4つの市の男女別人口の主成分の計算 R (Console) > data <- rbind( # データの入力 + c(18.8, 18.7), # 豊橋市 + c( 6.3, 6.6), # 瀬戸市 + c(22.2, 20.0), # 豊田市 + c( 9.4, # 安城市 9.0) + ) # 主成分の計算 > prcomp(data) 主成分分析の計算はこの1行だけ Rをつかうととても簡単に主成分分析ができる

33 主成分分析の結果の見方 4つの市の男女別人口の主成分の計算結果 R (Console) Standard deviations (1, .., p=2): 主成分得点の標準偏差 [1] 10.0813386 第1主成分 0.5476735 第2主成分 Rotation (n x k) = (2 x 2) PC1 PC2 [1,] -0.7451377 0.6669106 [2,] -0.6669106 -0.7451377 第1主成分 主成分が持つ情報量 第2主成分 固有ベクトル(変換行列) 主成分(Principal Component)の向き(傾き) 第1主成分の固有ベクトルがExcelで 計算した結果と同じであることを確認しよう

34 分析するデータ TSUTAYA DISCAS週間総合ランキング上位20位 在庫枚数 … ランキング レンタル開始日 STRAY SHEEP / 米津玄師 1 2020/08/22 358 5807 … Traveler(通常盤) / Official髭男dism 2 2019/10/26 364 218 … eyes(通常盤) / milet 3 2020/06/20 155 158 … クロマティカ / レディー・ガガ 4 2019/08/29 45 25 … 1位登録者数 …

35 データの作成 1 順序に意味がある(数値)データを選択 順序に意味がある ランキング 在庫数 全国評価 1位登録者数 レンタル開始日 2位登録者数 順序に意味がない(数値でない) ジャンル 2 日付を1900年元日からの日数に変換 2020年08月22日 44065

36 データの種類 全国評価 質的データ ランキング 分類や区別のためのデータ ジャンル レンタル開始日 在庫数量 1位登録者数 2位登録者数 比較のために 工夫して数値化 量的データ 数値として意味のあるデータ 足し引きができるデータ

37 作業ディレクトリの確認と変更 作業ディレクトリ working directory Rのデータの読み込み先となるディレクトリ(フォルダ) 作業ディレクトリの確認 > getwd() R (Console) # 作業ディレクトリの確認コマンド [1] "E:/Workspace" 作業ディレクトリの変更 > setwd("E:/Workspace") R (Console) # 作業ディレクトリの変更コマンド

38

データの読み込み
作業ディレクトリに演習ファイル「tsutaya_rental.csv」を配置

CSVデータの読み込み

R (Console)

> data <- readcsv(

# データの読み込み

+

"tsutaya_rental.csv",

#

ファイル名

+

header=T,

#

ヘッダ行の有無(T: 有, F: 無)

+

row.names=1)

#

行名の指定(列番号)

39 データの確認 変数の確認(R Studio) 変数の確認 R (Console) > data # CSVデータの確認 CD名...アーティスト名 1 STRAY SHEEP / 米津玄師 2 Traveler(通常盤) / Official髭男dism ...

40 スケールが異なる指標をもつデータ 単位が異なる 0～6000 指標間のばらつきが 大きく異なる 主成分に指標が反映されない 0～400

41 変数の標準化 女性(𝑥2 )／万人 1位登録者数(𝑥2 )／人 25 25 20 20 15 15 10 10 5 5 𝑂 5 10 15 20 25 男性(𝑥1 )／万人 スケールが同じ 𝑂 1位登録者数(𝑥2 ) 0 5 10 15 20 25 在庫数(𝑥1 )／枚 スケールが異なる ス ケ ー リ ン グ 0 在庫数(𝑥1 ) スケールを合わせる

42 相関行列 correlation matrix 1 𝑟(𝑥1 , 𝑥2 ) ⋯ 𝑟(𝑥1 , 𝑥𝑝 ) 𝑟(𝑥1 , 𝑥2 ) 1 ⋯ 𝑟(𝑥2 , 𝑥𝑝 ) ⋮ 𝑟(𝑥1 , 𝑥𝑝 ) 𝑟(𝑥2 , 𝑥𝑝 ) ⋯ 1 ℝ= 相関係数𝑟(𝑥𝑖 , 𝑥𝑗 ) 相関係数 correlation coefficient 𝑟 𝑥𝑖 , 𝑥𝑗 = 𝑟 𝑥𝑗 , 𝑥𝑖 より，対称行列 𝑖 = 𝑗のとき𝑟 𝑥𝑖 , 𝑥𝑗 = 1 相関係数の定義から，ℝ = 分散共分散係数 -1以上1以下の値を取り 2変量𝑥，𝑦間の線形な関係の強弱を示す指標 𝐶(𝑥, 𝑦) 𝑟(𝑥, 𝑦) = 𝜎 𝑥 𝜎(𝑦) Σ= Σ 𝜎 𝑥𝑖 𝜎 𝑥𝑗 variance-covariance matrix 𝜎 2 (𝑥1 ) 𝐶(𝑥1 , 𝑥2 ) ⋯ 𝐶(𝑥1 , 𝑥𝑝 ) 𝐶(𝑥1 , 𝑥2 ) 𝜎 2 (𝑥2 ) ⋯ 𝐶(𝑥2 , 𝑥𝑝 ) ⋮ 𝐶(𝑥1 , 𝑥𝑝 ) 𝐶(𝑥2 , 𝑥𝑝 ) ⋯ 𝜎 2 (𝑥𝑝 )

43 固有ベクトルから解釈する主成分 スケーリングと分析結果の表示 R (Console) > result <-prcomp(data, scale=T) # スケーリング(scale)して主成分分析 > result # 変数の内容を表示 ... 値の大きい成分から解釈 PC1 PC2 ... 0.41451571 -0.1070368 ... レンタル開始日 -0.04882473 0.6441358 ... 在庫枚数 -0.49154089 -0.2935740 ... 1位登録者 -0.51055576 0.2164189 ... 2位登録者 -0.50333112 0.2418657 ... 全国評価 -0.26489815 -0.6181508 ... ランキング 第1主成分は登録者の少なさ 第2主成分は レンタル開始日の新しさと 全国評価の低さを併せた指標

44 主成分負荷量 principal component loading 変量と主成分の相関係数 主成分がどの変量と強く関係しているかがわかる 変量𝑥𝑖 に対する第𝑚主成分の主成分負荷量𝑟(𝑥𝑖 , 𝑤𝑚 ) (𝑖) スケーリングをしていない場合 スケーリングをしていない場合 𝑟(𝑥𝑖 , 𝑤𝑚 ) = 𝑟(𝑥𝑖 , 𝑤𝑘 ) = 𝜆𝑚 𝑤𝑚 𝜎 𝑥𝑖 𝜆𝑚 (𝑖) 𝑤𝑚 𝜆𝑘 𝑤𝑚 (𝑖) 第𝑚主成分の固有値(分散) 第𝑚主成分の固有ベクトルの 第𝑖成分 ※ここまで定義してきた記号とは一部違います． ※なぜ相関係数となるのか？については説明しません．ぜひ自分で確認してみてください．

45 主成分負荷量から解釈する主成分 主成分負荷量の表示 R (Console) # 主成分負荷量を表示 > sweep(result$rotation, 2, result$sdev, "*") PC1 PC2 ... 0.69881805 -0.1446553 ... レンタル開始日 -0.08231196 0.8705199 ... 在庫枚数 -0.82867220 -0.3967518 ... 第1主成分は 1位登録者 -0.86072872 0.2924802 ... 登録者の少なさに加えて 2位登録者 -0.84854894 0.3268704 ... 在庫枚数の少なさも強く関係 全国評価 -0.44658284 -0.8354024 ... ランキング 固有ベクトルの成分の値，主成分負荷量ともに大きな値

46 分析結果の2次元プロット 2次元プロット R (Console) > par(family="HiraKakuProN-W3") # 描画用フォントの設定(Macユーザ向け) > biplot(result) # 2次元平面にプロット > rownames(result$x) <- NULL # 行名の削除(要再プロット)

47 2次元プロットの見方 固有ベクトル 第1主成分が比較的低い 第2主成分が比較的高い 第1主成分 第2主成分 においては平均的 第2主成分が比較的低い 主成分得点 変量(回転前の軸)

48 2次元プロットの解釈 評レ 価ン がタ 低ル い開 始 さ れ て か ら 日 が 浅 く 「STRAY SHEEP / 米津玄師」は レンタル開始直後で評価は安定していないが ランキング上位登録者が多く在庫も多数ある 「Traveler / Official髭男dism」などは レンタル開始から日が経っているが 比較的評価が高い位置に留まっている 第2主成分が比較的低い 「紅蓮華 / LiSA」などは 指標全体から平均的なグループであり 評価も平均的な値になっている ランキング登録者が少なく 在庫枚数も多く確保されていない

49 スクリープロットによる次元削減 2.5 固 有 値 2.0 1.5 1.0 0.5 0 1 2 3 4 5 6 第𝑚主成分 固有値が落ち込んでいる第4主成分までを利用

50 寄与率による次元削減 寄与率を使う 説明できている情報量の割合で次元を削減 寄与率，累積寄与率の表示 R (Console) # 寄与率，累積寄与率の表示 > summary(result) Importance of components: PC1 PC2 PC3 PC4 ... 1.6859 1.3515 0.9503 0.52170 ... Proportion of Variance 0.4737 0.3044 0.1505 0.04536 ... Cumlative Proportion 0.4737 0.7781 0.9286 0.97396 ... Standard deviation 寄与率 累積寄与率 80%以上情報が表現されている 省略

51 主成分分析を実践してみよう 次のデータの好きな方を使って主成分分析をおこない 主成分に対する解釈を考えてみよう 1 愛知県の市区町村における事業所数 2 伊藤研卒業生と関わりのあるIT企業10社 興味のある人は教育用標準データセット2020年版(SSDSE-2020A.xlsx)などの オープンデータから分析をしてみよう

まとめ

53 主成分分析の理解しておくべきポイント 1 データを効率的に説明できる扱いやすい形に 変換したい場合に主成分分析を利用しよう※ 2 主成分は分散共分散行列・相関行列から得られる 3 主成分の解釈・言語化は容易ではなく 望ましい結果が毎回必ず得られるわけではない ※共通する要因を見つけたいときは因子分析おこなおう