--- title: "第6講 主成分分析" subtitle: "基本的な考え方" date: "`r Sys.Date()`" format: html: toc: true html-math-method: katex self-contained: true grid: margin-width: 350px execute: echo: true warning: false reference-location: margin citation-location: margin tbl-cap-location: margin fig-cap-location: margin editor: visual --- ## 準備 以下で利用する共通パッケージを読み込む． ```{r} library(conflicted) # 関数名の衝突を警告 conflicts_prefer( # 優先的に使う関数を指定 dplyr::filter(), dplyr::select(), dplyr::lag(), ) library(tidyverse) library(broom) # 解析結果を tibble クラスに集約 library(gt) # 表の作成 library(gtsummary) # 分析結果の表の作成 library(ggfortify) # 分析結果の描画 library(GGally) library(ggrepel) ``` ## 主成分分析に用いる主な関数 主成分分析に必要な関数を説明する． Rの標準的な関数として - `stats::prcomp()` - `stats::princomp()` があり，計算法に若干の違いがある． - 数値計算の観点からみると関数 `prcomp()` が優位 - 関数 `princomp()` はS言語(R言語の元となる商用言語)との互換性を重視した実装 本講義では関数 `stats::prcomp()` を利用する． ::: callout-note ### 主成分分析の実行 関数 `stats::prcomp()` には以下の2通りの使い方がある． データフレームの全ての列を用いる場合 ``` r prcomp(x, retx = TRUE, center = TRUE, scale. = FALSE, tol = NULL, rank. = NULL, ...) #' x: 必要な変数のみからなるデータフレーム #' center: 中心化(平均0)を行って処理するか否か #' scale.: 規格化(分散1)を行って処理するか否か ``` 式(formula)を用いて列を指定する場合 ``` r prcomp(formula, data = NULL, subset, na.action, ...) #' formula: ~ 変数名 (解析の対象を + で並べる) 左辺はないので注意 #' data: 必要な変数を含むデータフレーム #' 詳細は '?stats::prcomp' を参照 ``` ::: 分析結果を参照する方法は `base R` と `tidyverse` で さまざまな形で用意されている． ::: callout-note ### 結果を得る基本関数 分析結果の概要を確認するには， 関数 `stats::prcomp()` の返値を関数 `base::print()` で単に表示するか， 関数 `base::summary()` を用いれば良い． ``` r print(x, ...) #' x: prcomp の返値 #' command line では単に print は書かなくてもよい summary(object, ...) #' object: prcomp の返値 #' 詳細は '?base::summary' を参照 ``` 関数 `stats::prcomp()` の返値にどのような情報が含まれているかを見るには 関数 `pillar::glimpse()` を利用する． ``` r glimpse(x, width = NULL, ...) #' x: prcomp の返値 #' width: 表示の長さを指定．既定値(NULL)は表示環境が可能な範囲 #' 詳細は '?pillar::glimpse' を参照 #' 関数 utils::str() も同様な働きをする ``` ::: 分析結果を `tibble` クラスで整理して取得するには `package::broom` に用意されている以下の2つの関数を利用する． ::: callout-note ### 主成分得点・負荷・寄与率の取得 ``` r tidy(x, matrix = "u", ...) #' x: prcomp が出力したオブジェクト #' matrix: 結果として取り出す行列 #' 主成分得点: "u", "samples", "scores" または "x" (既定値) #' 主成分負荷: "v", "rotation", "loadings" または "variables" #' 寄与率: "d", "eigenvalues" または "pcs" #' 詳細は '?broom::tidy.prcomp' を参照 ``` ::: ::: callout-note ### 主成分得点の計算 ``` r augment(x, data = NULL, newdata, ...) #' x: prcomp が出力したオブジェクト #' data: 元のデータ (通常不要) #' newdata: 新たに主成分得点を計算するデータフレーム #' 詳細は '?broom::augment.prcomp' を参照 ``` ::: ::: column-margin 主成分得点を計算する関数は `base R` にも用意されている． ``` r predict(object, newdata, ...) #' object: prcomp が出力したオブジェクト #' newdata: 主成分得点を計算するデータフレーム #' 詳細は '?stats::prcomp' または '?stats::predict.prcomp' を参照 ``` ::: 視覚化は関数 `broom::augment()` などで取得したデータフレームを用いて行うことができるが， 主成分分析において標準的なバイプロットは関数 `ggfortify::autoplot()` を用いれば良い． ::: callout-note ### バイプロット ``` r autoplot( object, data = NULL, scale = 1, x = 1, y = 2, variance_percentage = TRUE, ... ) #' object: prcomp が出力したオブジェクト #' data: 分析に用いたデータフレーム(分析に用いた列以外を利用する場合) #' scale: 得点と負荷のスケール #' x,y: xy軸に用いる主成分．規定値は第1と第2主成分 #' variance_percentage: 寄与率の表示 #' 詳細は '?ggfortify::autoplot.prcomp' を参照 ``` ::: ## データセット 今回配布したデータセットの内容は以下のとおりである． ::: callout-tip ### `japan_social.csv` の概要 総務省統計局より取得した都道府県別の社会生活統計指標の一部 - Pref : 都道府県名 - Forest : 森林面積割合 (%) 2014年 - Agri : 就業者１人当たり農業産出額(販売農家）(万円) 2014年 - Ratio : 全国総人口に占める人口割合 (%) 2015年 - Land : 土地生産性（耕地面積１ヘクタール当たり）(万円) 2014年 - Goods : 商業年間商品販売額［卸売業＋小売業］（事業所当たり）(百万円) 2013年 - Area : 地方区分 - 参考 : データを読み込む際に， 視覚化に利用するために 地方区分を因子化しておくと良い． ```{r} #| eval: false js_data <- read_csv("data/japan_social.csv") |> mutate(Area = as_factor(Area)) # 地方区分を因子化 ``` ::: ## 人工データによる主成分分析 ### 問題 数値実験により主成分分析の考え方を確認しなさい． - 以下のモデルに従う人工データを生成する ```{r} #| eval: false #' 観測データ (2次元) の作成 (aのスカラー倍に正規乱数を重畳) a <- (c(1,2)/sqrt(5)) |> # 主成分負荷量 (単位ベクトル) set_names(c("x1","x2")) # 後の行列計算のために成分名を付ける n <- 100 # データ数 #' 行列として作成して tibble クラスに変換 toy_data <- # aのスカラー倍に正規乱数を重畳 as_tibble(runif(n,-1,1) %o% a + rnorm(2*n, sd=0.3)) ``` - 観測データの散布図を作成する - 観測データから第1主成分負荷量を推定する ```{r} #| eval: false toy_pca <- toy_data |> prcomp() # 全ての列を用いて主成分を計算する a_hat <- toy_pca |> tidy("loadings") |> # 負荷量の取得 filter(PC == 1) |> # 第1主成分の選択 pull(value) # 値のみベクトルとして取得 ``` - 散布図上に主成分負荷量を描画する 参考 : 主成分方向を図示するには，傾きを計算して直線を引けば良い ``` r geom_abline(slope = 傾き, intercept = 切片) # 指定の直線を追加できる ``` ### 解答欄 ```{r} ``` ### 解答例 2次元の人工データを生成する． ```{r} set.seed(123123) n <- 100 # データ数 a <- (c(1,2)/sqrt(5)) |> # 主成分負荷量(単位ベクトル)の設定 (適宜変更せよ) set_names(paste0("x",1:2)) # 規則的な文字列を生成する例 toy_data <- # aのスカラー倍に正規乱数を重畳(行列として作成) as_tibble(runif(n,-1,1) %o% a + rnorm(2*n, sd=0.3)) ``` 生成したデータの散布図を描く． ```{r} gg <- # 図を上描きして用いるためにオブジェクトに代入する toy_data |> ggplot(aes(x = x1, y = x2)) + geom_point(colour = "blue", shape = 4) + geom_abline(slope = a[2]/a[1], # 真の傾き=真の負荷量のy成分/x成分 intercept = 0, colour = "red") + # 主成分負荷量の図示 xlim(c(-2,2)) + ylim(c(-2,2)) # xy軸を揃えて直交関係を見易くする #' 関数 coord_fixed() を用いて縦横比を1としてもよいが，領域はデータに依存するので注意する print(gg) ``` $a$ 方向 $(1,2)$ に本質的な情報が集約されていることがわかる． 主成分負荷量の推定を行う． ```{r} toy_pca <- toy_data |> prcomp() a_hat <- toy_pca |> tidy("loadings") |> # 負荷量の取得 filter(PC == 1) |> # 第1主成分の選択 pull(value) |> # 値のみベクトルとして取得 set_names(paste0("x",1:2)) # 成分に名前を付けておく ``` 第1主成分負荷量が $a$ に非常に近いことが確認できる． 乱数によっては符号が反対になることもある． 前の散布図上に推定された方向を描画することで，両者が近いことが視覚的にも確認できる． ```{r} gg + geom_abline(slope = a_hat[2]/a_hat[1], intercept = 0, colour = "orange", linetype = "dashed") ``` 第1主成分得点を取得し，各データの第1主成分ベクトルを計算する． ```{r} pc1 <- augment(toy_pca) |> # 主成分得点のデータフレーム pull(.fittedPC1) # ベクトルとして抽出 toy_pc1 <- as_tibble(pc1 %o% a_hat) # 元の2次元に戻す ``` 各データの第1主成分を元の散布図上で図示する． ```{r} gg + geom_point(data = toy_pc1, aes(x = x1, y = x2), colour = "purple", shape = 18) ``` ## 第1主成分の求め方 ### 問題 第1主成分とGram行列の固有ベクトルの関係を調べなさい． - 人工データを生成する - 主成分分析を実行する - Gram 行列を計算し固有値・固有ベクトルを求める ```{r} #| eval: false #' 中心化を行う X <- scale(toy_data, scale = FALSE) #' 詳細は '?base::scale' を参照 #' Gram 行列を計算する G <- crossprod(X) #' 固有値・固有ベクトルを求める eigen(G) # 返り値 'values, vectors' を確認 #' 詳細は '?base::eigen' を参照 ``` ### 解答欄 ```{r} ``` ### 解答例 3次元の人工データを作成する． ```{r} set.seed(242424) n <- 50 # データ数 d <- 3 a <- rnorm(d) |> set_names(paste0("x",1:d)) # 成分名(x1,...,xd)を付与 (a <- a/sqrt(sum(a^2))) # 主成分負荷量(単位ベクトル)を生成 toy_data <- as_tibble(runif(n, -1, 1) %o% a + rnorm(d*n, sd = 0.1)) ``` 散布図行列を用いて作成したデータの視覚化を行う． ```{r} ggpairs(toy_data) ``` 推定された第1主成分負荷量を取得する． ```{r} toy_pca <- prcomp(toy_data) pc1 <- augment(toy_pca) |> pull(.fittedPC1) ``` ::: column-margin 以下の図は3次元のときのみ実行可能である． ```{r} library(scatterplot3d) s3d <- scatterplot3d(toy_data, type = "h", ## asp=1, # 軸の比率を揃える場合 highlight.3d = TRUE) s3d$points3d(pc1 %o% a, col = "blue") ``` ::: 主成分負荷量の推定を固有値分解と比較する． ```{r} toy_eigen <- # 固有値分解 eigen(crossprod(scale(toy_data, scale = FALSE))) toy_pca$rotation # 主成分負荷量 toy_eigen$vectors # 固有ベクトル (符号を除いて主成分負荷量と一致) toy_pca$sdev # 主成分の標準偏差 sqrt(toy_eigen$values/(n-1)) # 固有値と主成分の標準偏差の関係 ``` ::: column-margin `package::broom` を利用して主成分負荷量や標準偏差を 行列やベクトルとして取得するには以下のようにすればよい． ```{r} #' 主成分負荷量 toy_pca |> tidy("loadings") |> # 負荷量を取得 pivot_wider(names_from = PC, # 主成分ごとの列に変換 names_prefix = "PC") #' 標準偏差 toy_pca |> tidy("eigenvalues") |> # 標準偏差を取得 pull(std.dev) # 値のみ取り出す ``` `tibble`クラスを`matrix`クラスに変換するには以下のようにすればよい． ```{r} toy_pca |> tidy("loadings") |> pivot_wider(names_from = PC, names_prefix = "PC") |> select(starts_with("PC")) |> as.matrix() ``` ::: ## 実データによる分析 ### 問題 データセット `japan_social.csv` を用いて主成分分析を行いなさい． - データを読み込む - データの散布図行列を描く - 各データの箱ひげ図を描き，変数の大きさを確認する - 主成分負荷量を計算する ```{r} #| eval: false js_pca <- js_data |> select(!c(1,7)) |> # '!c(1,7)' は都道府県名・地方区分を除く prcomp(scale. = TRUE) # 'scale.=TRUE' とすると変数を正規化してから解析する ``` ### 解答欄 ```{r} ``` ### 解答例 データの読み込みを行う． ```{r} js_data <- read_csv("data/japan_social.csv") |> mutate(Area = as_factor(Area)) # 地方区分を因子化 ``` 散布図を用いてデータを視覚化する． ```{r} js_data |> select(!c(Pref,Area)) |> # 都道府県名・地方区分は削除 ggpairs() # 散布図．いくつかの変数は相関が高いことがわかる ``` ::: column-margin 地方ごとに色を変える場合は以下のようにすればよい． ```{r} js_data |> ggpairs(columns = 2:6, # 都道府県名・地方区分は除く legend = c(2,1), # 2行1列のグラフから凡例を作成 lower = list(continuous = wrap("points", alpha = .5), mapping = aes(colour = Area))) ``` ::: 箱ひげ図を用いて変数のばらつきを調べると，変数ごとのばらつきに大きな違いがあることがわかる． ```{r} js_data |> pivot_longer(where(is.double)) |> # 実数値(都道府県名・地方区分以外)をまとめる mutate(name = as_factor(name)) |> # name列に現れる順番どおりboxplotを並べるため ggplot(aes(x = name, y = value)) + # 既定値の name と value をxy軸に用いる geom_boxplot(aes(fill = name), # 変数ごとに色を変える show.legend = FALSE) # 凡例は表示しない ``` ::: column-margin 箱ひげ図で判り難い場合には violin plot を利用して密度を表示しても良い． ```{r} js_data |> pivot_longer(where(is.double)) |> mutate(name = as_factor(name)) |> ggplot(aes(x = name, y = value)) + geom_violin(aes(fill = name), show.legend = FALSE) ``` ::: 密度表示を行うと以下のようになる． ```{r} js_data |> pivot_longer(where(is.double)) |> # 都道府県名・地方区分以外をまとめる mutate(name = as_factor(name)) |> ggplot(aes(x = value)) + # value を x 軸に用いる geom_density(aes(fill = name), # 変数ごとに色を変える alpha = 0.4) # 重なっても良いように半透明にする ``` 変数ごとに標準化(平均0分散1)したデータフレームを作成するには以下のようにすれば良い． ```{r} js_data |> mutate(across(where(is.double), \(x)c(scale(x)))) |> ## 実数値(都道府県名・地方区分以外)の列を標準化する pivot_longer(where(is.double)) |> mutate(name = as_factor(name)) |> ggplot(aes(x = name, y = value)) + geom_boxplot(aes(fill = name), show.legend = FALSE) ``` ::: column-margin データフレームの特定の列のみ変換を行うには関数 `dplyr::across()` を利用する． 関数を組み合わせることで様々な条件が表現できる． 具体的な例は `?dplyr::across` や `?tidyselect::where` を参照のこと． 関数 `base::scale()` は常に行列を返すため，そのまま用いると期待どおりに動かない． この例ではベクトルを作成する関数 `base::c()` (関数 `base::as.vector()` でも可)と無名関数を用いてこれを解消している． 無名関数の代わりにラムダ式 `~c(scale(.))` などを用いることもできる． ::: ::: column-margin 標準化したデータを有効数字3桁で表示するには，例えば以下のようにすればよい． ```{r} js_data |> # mutate(across(where(is.double), \(x)signif(c(scale(x)), digits = 3))) ``` ::: 標準化したデータの密度推定は以下のようになる． ```{r} js_data |> mutate(across(where(is.double), \(x)c(scale(x)))) |> pivot_longer(where(is.double)) |> mutate(name = as_factor(name)) |> ggplot(aes(x = value)) + geom_density(aes(fill = name), alpha = 0.4) ``` 主成分分析を行う． 変数の標準化は関数 `stats::prcomp()` にオプションを指定することで実行できるので， 元のデータをそのまま渡せば良い． ```{r} js_pca <- js_data |> select(where(is.double)) |> # 実数値の列(都道府県名・地方区分)を抽出 prcomp(scale. = TRUE) # 変数のばらつきを規格化して分析 js_pca # 結果の表示．関数 print() を用いてもよい ``` 主成分負荷量を表にするには，例えば以下のようにすればよい． ```{r} js_pca |> tidy("loadings") |> # "v", "rotation", "loadings" または "variables" pivot_wider(names_from = PC, # 横長の形式に変換する names_prefix = "PC", # PC列の番号で新しい列を作る values_from = value) |> # valueは既定値なので指定しなくても良い gt() |> fmt_number(decimals = 3) ``` 主成分方向から読み取れることは 1. 第1: 人の多さに関する成分(正の向きほど人が多い) 2. 第2: 農業生産力に関する成分(正の向きほど高い) 第1，第2主成分得点を利用して2次元の地図を作成することができる． 関数 `ggrepel::geom_text_repel()` を用いると，ある程度自動的に文字が重ならないように表示することができる． ```{r} js_pca |> augment(data = js_data) |> # データと主成分得点のデータフレーム ggplot(aes(x = .fittedPC1, y = .fittedPC2)) + geom_point(aes(colour = Area), # 地方区分ごとに色を変える shape = 19, size = 2) + geom_text_repel(aes(label = Pref), colour = "darkgray") ``` ::: column-margin 関数 `ggplot2::geom_text()` を用いる場合は， 文字と記号の重なりを防ぐためのオプションを利用する． ```{r} js_pca |> augment(data = js_data) |> # データと主成分得点のデータフレーム ggplot(aes(x = .fittedPC1, y = .fittedPC2)) + geom_point(aes(colour = Area), # 地方区分ごとに色を変える shape = 19, size = 2) + geom_text(aes(label = Pref), # 県名を重ね書きする colour = "darkgray", hjust = 0, nudge_x = 0.1, check_overlap = TRUE) ``` ::: 主成分得点による地図は関数 `ggfortify::autoplot()`を利用することもできる． ```{r} js_pca |> autoplot(scale = 0, # バイプロットの設定 (主成分得点での表示) data = js_data, # 必要な情報を含むデータフレーム colour = "Area", # 地方区分ごとに色付 label = TRUE, # ラベルを付加 label.label = "Pref", # 都道府県名を追加 label.repel = TRUE) # 横にずらして表示 ``` 寄与率を表示する．詳しくは次週説明する． ```{r} summary(js_pca) ``` ::: column-margin 2変数での解析例 (AgriとLandを取り上げる，その他の組み合わせでも試みよ) は以下のようになる. 多変数での視覚化は次週詳しく説明する． ```{r} js_pca2 <- js_data |> select(c(Agri,Land)) |> # 必要な変数だけのデータフレーム prcomp(scale. = TRUE) a_hat <- tidy(js_pca2, "loadings") |> filter(PC == 1) |> pull(value, name = column) # column列から名前を取ったベクトルを作成 js_pc1 <- as_tibble((augment(js_pca2) |> pull(.fittedPC1)) %o% a_hat) js_data |> mutate(across(where(is.double), \(x)c(scale(x)))) |> ggplot(aes(x = Agri, y = Land)) + geom_point(colour = "blue", shape = 4) + geom_abline(slope = a_hat[2]/a_hat[1], # 主成分負荷量(方向)の図示 intercept = 0, colour = "orange") + geom_point(data = js_pc1, # 第1主成分を元の散布図上で図示 colour = "purple", shape = 18) + coord_fixed() # 縦横比を1に指定 ``` ::: ::: column-margin 配布したファイル `prefecture.csv` を利用すると都道府県名などを日本語にして分析することができる． ```{r} js_data <- bind_cols( # 2つのデータフレームを結合 read_csv(file = "data/prefecture.csv", col_select = c(2,4)) |> set_names("県名","地方区分"), # 列の名称を"県名"と"地方区分"に変更 read_csv(file = "data/japan_social.csv", col_select = 2:6) |> set_names("森林面積割合","農業算出額","人口割合","土地生産性","商品販売額") )|> # 簡略化した項目名に変更 mutate(地方区分 = as_factor(地方区分)) # 地方区分を出現順に因子化 #' macOSのための日本語表示の設定(ここから) if(Sys.info()["sysname"] == "Darwin") { # macOSか調べる jp_font <- "HiraMaruProN-W4" theme_update(text = element_text(family = jp_font)) } else {jp_font <- NULL} #' macOSのための日本語表示の設定(ここまで) js_data |> select(where(is.double)) |> prcomp(scale. = TRUE) |> autoplot(asp = 1, # 縦横比を設定 data = js_data, colour = "地方区分", # 地方ごとに色付け label = TRUE, # ラベルの表示 label.label = "県名", label.repel = TRUE, # ラベルの表示を自動調整 label.family = jp_font, label.size = 3, loadings = TRUE, # 負荷の表示 loadings.colour = "orchid", loadings.label = TRUE, loadings.label.colour = "darkgray", # 負荷ラベルの表示 loadings.label.family = jp_font, loadings.label.size = 4) ``` :::