---
title: "第7講 主成分分析"
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(GGally)  
library(broom)      # 解析結果を tibble 形式に集約
library(gt)         # 表の作成
library(gtsummary)  # 分析結果の表の作成
library(ggfortify)  # バイプロット表示のため
library(ggrepel)
#' macOSのための日本語表示の設定
if(Sys.info()["sysname"] == "Darwin") { # macOSか調べる
    jp_font <- "HiraMaruProN-W4"
    theme_update(text = element_text(family = jp_font))
    update_geom_defaults("text", list(family = theme_get()$text$family))
    update_geom_defaults("label", list(family = theme_get()$text$family))
    update_geom_defaults("text_repel", list(family = theme_get()$text$family))
    update_geom_defaults("label_repel", list(family = theme_get()$text$family))
} else {jp_font <- NULL}
```

## 主成分分析に用いる主な関数

分析に用いる主な関数としては以下がある．

-   `stats::prcomp()`: 主成分分析の実行
-   `base::summary()`: 主成分負荷・寄与率を表示
-   `broom::tidy()`: 主成分得点・負荷・寄与率を tibble クラスで取得
-   `broom::augment()`: tibble クラスで主成分得点を計算
-   `ggfortify::autoplot()`: 主成分得点による散布図およびバイプロット

``` r
#' データフレーム 'toy_data' を分析
toy_pca <-
    toy_data |>
    select(計算対象の列) |> # 必要な列を選択するか不要な列を削除
    prcomp(必要なら標準化)  # 'scale.=TRUE' で標準化 
#' 分析結果の確認
summary(toy_pca)   # 主成分負荷量と寄与率を確認
toy_pca |>
    tidy("eigenvalues") # "scores","loadings","eigenvalues" で指定(別表記もあり)
#' 第1,2主成分得点で散布図を描く
toy_pca |>
    autoplot(scale = 0, x = 1, y = 2) 
```

::: callout-tip
### `jpamenity.csv` を用いた例

人口関連データを整理する．

```{r}
ja_data <-
    bind_cols(
        ## 都道府県名と地方名の取得
        read_csv(file = "data/prefecture.csv",
                 col_select = c(2,4)) |>
        set_names("都道府県名", "地方区分"), # 列の名称を"都道府県名"と"地方区分"に変更
        ## データの取得
        read_csv(file = "data/jpamenity.csv",
                 col_select = !1:2) |> # 不要な列を読み込まない
        slice(-1) |>                   # 不要な行を削除
        set_names(names(read_csv(file = "data/jpamenityitem.csv"))) # 簡略化した項目名に変更
    ) |> 
    select(1:10) |> # 人口関連のみ利用
    mutate(地方区分 = as_factor(地方区分)) # 地方区分を出現順に因子化
```

散布図により変数間の関係を確認する．

```{r}
ja_data |>
    ggpairs(columns = 3:10,  # 都道府県名・地方区分を除く
            legend = c(2,1), # 2行1列のグラフから凡例を作成
            upper = list(continuous = wrap("cor", size = 3.5)),
            lower = list(continuous = wrap("points", alpha = .5),
                         mapping = aes(colour = 地方区分))) +
    theme(text = element_text(size = 8)) 
```

主成分分析を行い，バイプロット(主成分得点の散布図と主成分方向の表示)を描画する．

```{r}
ja_pca <- ja_data |> 
    select(!1:2) |>                              # 人口関連データのみ
    prcomp(scale. = TRUE)                        # 主成分分析の実行
ja_pca |>
    autoplot(asp = 1,                            # 縦横比を設定
             data = ja_data,
             colour = "地方区分",                 # 地方ごとに色付け
             label = TRUE,                       # ラベルの表示
             label.label = "都道府県名", 
             label.repel = TRUE,                 # ラベルの表示を自動調整
             label.size = 3,
             loadings = TRUE,
             loadings.colour = "orchid",         # 負荷の表示
             loadings.label = TRUE,
             loadings.label.colour = "darkgray", # 負荷ラベルの表示
             loadings.label.size = 4) 
```
:::

## データセット

以下では2つのデータセットを使用する

::: callout-note
### `japan_social2.csv` の概要

総務省統計局より取得した都道府県別の社会生活統計指標の一部

`japan_social.csv` を日本語化したもの

-   都道府県名 :
-   地方区分 :
-   森林面積割合: 森林面積割合 (%) 2014年
-   農業産出額 : 就業者１人当たり農業産出額(販売農家）(万円) 2014年
-   人口割合 : 全国総人口に占める人口割合 (%) 2015年
-   土地生産性 : 土地生産性（耕地面積１ヘクタール当たり）(万円) 2014年
-   商品販売額 : 商業年間商品販売額［卸売業＋小売業］（事業所当たり）(百万円) 2013年
-   参考 : <https://www.e-stat.go.jp/SG1/estat/List.do?bid=000001083999&cycode=0>

データの読み込み方の例．

```{r}
#| eval: false
js_data <-
    read_csv("data/japan_social2.csv",
             col_types = list(地方区分 = col_factor())) # 地方区分を因子化
```
:::

::: column-margin
列の `type` を指定したい場所が判っていて，先頭の方であれば以下でも可．

```{r}
#| eval: false
#' 1列目は自動判定(?)，2列目は因子(f)，3列目以降は自動判定
js_data <-
    read_csv("data/japan_social2.csv", col_types = "?f")
```
:::

::: callout-note
### `MASS::UScereal` の概要

**Nutritional and Marketing Information on US Cereals**

The UScereal data frame has 65 rows and 11 columns. The data come from the 1993 ASA Statistical Graphics Exposition, and are taken from the mandatory F&DA food label. The data have been normalized here to a portion of one American cup.

詳細は `?MASS::UScereal` を参照

データの整理の仕方の例． 元のデータは data.frame クラスなので，tibble クラスに変換して整理しておく．

```{r}
#| eval: false
glimpse(MASS::UScereal)                    # 各変数の属性を確認する．
uc_data <-
    MASS::UScereal |>
    rownames_to_column(var = "product") |> # 行名の製品名を product 列に加える
    as_tibble()                            # tibble クラスに変換
```
:::

## 寄与率・累積寄与率

### 問題

データセット

-   `japan_social.csv`
-   `MASS::UScereal`

について，標準化の有無の違いで寄与率・累積寄与率がどのように異なるか確認しなさい．

``` r
toy_data |> prcomp()              # 標準化を行わない場合
toy_data |> prcomp(scale. = TRUE) # 標準化を行う場合
```

::: column-margin
正式なオプション名は `scale.` であるが，`sc = TRUE` など他のオプションと区別できれば短縮表記も可能．
:::

### 解答例

#### `japan_social2.csv` の場合

データを読み込む．

```{r}
js_data <-
    read_csv("data/japan_social2.csv",
             col_types = list(地方区分 = col_factor())) # 地方区分を因子化
```

データの性質を把握するために視覚化を行う．以下は全体の散布図および特定の2項目の散布図を描いた例である．

```{r}
js_data |> 
    ggpairs(columns = which(sapply(js_data,is.double)),    # 数値列を指定．列番号(3:7)でも指定可
            legend = c(2,1),                               # 2行1列のグラフから凡例を作成
            lower = list(continuous = wrap("points", alpha = .5),
                         mapping = aes(colour = 地方区分))) # 地方区分で色分け
```

いくつかの変数は相関が強いことがわかる．`農業算出額` と `森林面積割合` で散布図を確認する．

```{r}
js_data |>
    ggplot(aes(x = 農業算出額, y = 森林面積割合, label = 都道府県名)) +
    geom_point(colour = "blue") +
    geom_text_repel(colour = "brown")
```

主成分分析を実行する．

```{r}
js_pca_raw <-
    js_data |>
    select(where(is.double)) |> # 数値列を指定
    prcomp()                    # 標準化なし
js_pca <-
    js_data |> 
    select(where(is.double)) |> 
    prcomp(scale. = TRUE)       # 標準化あり
```

::: column-margin
以下のようにすると分析結果の得点の各行に県名 `都道府県名` をつけることができる．

```{r}
#| eval: false
js_pca_raw <- js_data |>
    column_to_rownames(var = "都道府県名") |> # '都道府県名'を行名に変換
    select(where(is.double)) |> 
    prcomp()                                # 標準化なし
js_pca <- js_data |> 
    column_to_rownames(var = "都道府県名") |> 
    select(where(is.double)) |> 
    prcomp(scale. = TRUE   )                # 標準化あり
```
:::

標準化しない場合の分析結果は以下のようになる． まず，寄与率および累積寄与率を確認する．

```{r}
summary(js_pca_raw) 
```

第1,2主成分でほとんど説明できることが示唆される．

主成分負荷量を取り出す．

```{r}
js_pca_raw$rotation 
```

負荷量が偏る傾向があり，各主成分はほぼ1つの変数に対応している．

寄与率および累積寄与率を tibble クラスで取得すれば， さまざまな表やグラフを作成することができる． 寄与率に関する情報は，各主成分ごとに標準偏差，寄与率，累積寄与率の列を持つデータフレームなので， 例えば以下のような表示が可能である．

```{r}
#' 寄与率に関する表を作成
js_pca_raw |>                              
    tidy("eigenvalues") |>                 # "d" または "eigenvalues" または "pcs"
    gt() |>      
    fmt_number(columns = !1, decimals = 3) # 1列目以外小数点以下3桁
#' 各主成分ごとに寄与率(percent)を表示
js_pca_raw |>
    tidy("eigenvalues") |> 
    ggplot(aes(x = PC, y = percent)) +     
    geom_bar(stat = "identity")
#' 各主成分ごとに累積寄与率(cumlative)を表示
js_pca_raw |>
    tidy("eigenvalues") |> 
    ggplot(aes(x = PC, y = cumulative)) +  
    geom_bar(stat = "identity")
#' まとめて並べて表示
js_pca_raw |>                              
    tidy("eigenvalues") |>
    pivot_longer(!PC) |>
    mutate(name = as_factor(name)) |>
    ggplot(aes(x = PC, y = value, fill = name)) +
    geom_bar(stat="identity") +
    facet_wrap(~ name, scales = "free_y") +
    theme(legend.position = "none") 
```

主成分負荷量に関する情報は， 主成分と特性量を指定して負荷量を得る型式の長い表なので， 目的に応じて以下のように整理する必要がる．

```{r}
#' 主成分負荷量に関する表を作成
js_pca_raw |>
    tidy("loadings") |>
    mutate(column = as_factor(column)) |>
    ggplot(aes(x = PC, y = value, fill = column)) +
    geom_bar(stat = "identity", width = 0.7, position = position_dodge(width = 0.7))
#' 各主成分負荷量の構成を表示
js_pca_raw |>                              
    tidy("loadings") |> 
    pivot_wider(names_from = PC,        
                names_prefix = "PC",    
                values_from = value) |> 
    gt() |>
    fmt_number(decimals = 3)               # 表示を制限
#' 各特性量が第1・第2主成分負荷量に占める方向を表示
js_pca_raw |>                              
    tidy("loadings") |> 
    pivot_wider(names_from = PC,        
                names_prefix = "PC",    
                values_from = value) |> 
    ggplot(aes(x = PC1, y = PC2, label = column)) +
    geom_point(colour = "blue", size = 2) + # 大きめの点にする
    geom_label_repel(colour = "royalblue")
```

簡単な散布図の表示方法として `package::ggfortify` を利用する方法は以下のようになる．

```{r}
js_pca_raw |>
    autoplot(asp =1,                   # 縦横の比率を1
             scale = 0,                # 主成分得点のままの散布図
             data = js_data,           # ラベルの情報を取得するデータ
             colour = "地方区分",       # 地方ごとに色付け
             label = TRUE,             # 以下はラベルの付け方を指定
             label.label = "都道府県名", # ラベルに用いる列名
             label.repel = TRUE)       # ラベルの重なりを低減
```

同様に標準化した場合の分析結果を見てみる．

```{r}
js_pca |>                                  # 寄与率に関する情報
    tidy("eigenvalues") |>                
    gt() |>      
    fmt_number(columns = !1, decimals = 3)
js_pca |>
    tidy("eigenvalues") |>
    pivot_longer(!PC) |>
    mutate(name = as_factor(name)) |>
    ggplot(aes(x = PC, y = value, fill = name)) +
    geom_bar(stat="identity") +
    facet_wrap(~ name, scales = "free_y") +
    theme(legend.position = "none") 
js_pca |>                                  # 主成分負荷量に関する情報
    tidy("loadings") |> 
    pivot_wider(names_from = PC,        
                names_prefix = "PC",    
                values_from = value) |> 
    gt() |>
    fmt_number(decimals = 3)               
js_pca |>
    tidy("loadings") |>
    mutate(column = as_factor(column)) |>
    ggplot(aes(x = PC, y = value, fill = column)) +
    geom_bar(stat = "identity", width = 0.7, position = position_dodge(width = 0.7))
js_pca |>
    autoplot(asp =1,
             scale = 0,
             data = js_data,
             colour = "地方区分",
             label = TRUE,
             label.label = "都道府県名",
             label.repel = TRUE)
```

::: column-margin
寄与率を表示するためには関数 `stats::screeplot()` を利用することもできる．詳細は `?screeplot` を参照．

```{r}
screeplot(js_pca) 
```
:::

#### `MASS::UScereal` の場合

データを整理する．また，分析を行う前に適当な方法で視覚化をすることを推奨する．

```{r}
uc_data <-
    MASS::UScereal |> 
    rownames_to_column(var = "product") |>
    as_tibble() 
uc_data |>
    ggpairs(columns = which(sapply(uc_data,is.double)), # 数値列
            legend = c(2,1), 
            upper = list(continuous = wrap("cor", size = 3.5)),
            lower = list(continuous = wrap("points", alpha = .5),
                         mapping = aes(colour = mfr)))
```

標準化なしの分析は以下のようになる．

```{r}
uc_pca_raw <-
    uc_data |>
    select(where(is.double)) |>
    prcomp() 
uc_pca_raw |> 
    tidy("eigenvalues") |>                
    gt() |>      
    fmt_number(columns = !1, decimals = 3)
uc_pca_raw |>
    tidy("eigenvalues") |>
    pivot_longer(!PC) |>
    mutate(name = as_factor(name)) |>
    ggplot(aes(x = PC, y = value, fill = name)) +
    geom_bar(stat="identity") +
    facet_wrap(~ name, scales = "free_y") +
    theme(legend.position = "none") 
uc_pca_raw |>
    tidy("loadings") |> 
    pivot_wider(names_from = PC,        
                names_prefix = "PC",    
                values_from = value) |> 
    gt() |>
    fmt_number(decimals = 3)               
uc_pca_raw |>
    tidy("loadings") |>
    mutate(column = as_factor(column)) |>
    ggplot(aes(x = PC, y = value, fill = column)) +
    geom_bar(stat = "identity", width = 0.7, position = position_dodge(width = 0.7))
uc_pca_raw |>
    autoplot(asp =1,
             scale = 0,
             data = uc_data,
             colour = "mfr",          # メーカー毎に色付け
             label = TRUE,
             label.label = "product", # 製品名のラベル
             label.repel = TRUE) +    # 重なりを低減
    labs(colour = "company")
```

標準化ありの分析は以下のようになる．

```{r}
uc_pca <-
    uc_data |>
    select(where(is.double)) |>
    prcomp(scale. = TRUE)
uc_pca |> 
    tidy("eigenvalues") |>                
    gt() |>      
    fmt_number(columns = !1, decimals = 3)
uc_pca |>
    tidy("eigenvalues") |>
    pivot_longer(!PC) |>
    mutate(name = as_factor(name)) |>
    ggplot(aes(x = PC, y = value, fill = name)) +
    geom_bar(stat="identity") +
    facet_wrap(~ name, scales = "free_y") +
    theme(legend.position = "none") 
uc_pca |>
    tidy("loadings") |> 
    pivot_wider(names_from = PC,        
                names_prefix = "PC",    
                values_from = value) |> 
    gt() |>
    fmt_number(decimals = 3)               
uc_pca |>
    tidy("loadings") |>
    mutate(column = as_factor(column)) |>
    ggplot(aes(x = PC, y = value, fill = column)) +
    geom_bar(stat = "identity", width = 0.7, position = position_dodge(width = 0.7))
uc_pca |>
    autoplot(asp = 1,
             scale = 0,
             data = uc_data,
             colour = "mfr",
             label = TRUE,
             label.label = "product",
             label.repel = TRUE) +
    labs(colour = "company")
```

## 主成分分析の視覚化

主成分分析の分析結果を可視化するバイプロット表示は， 関数 `ggfortify::autoplot()` にオプションを指定することで実行できる．

::: callout-note
### バイプロット表示

``` r
#' データフレームを分析する
toy_pca <-
    toy_data |>
    select(計算対象の列を指定) |>
    prcomp(必要なら標準化)
#' 第1と第2主成分を利用したバイプロット表示
toy_pca |>
    autoplot(label = TRUE,          # 各データのラベルを表示
             loadings = TRUE,       # 主成分負荷による成分方向の表示
             loadings.label = TRUE) # 成分名の表示
#' 第2と第3主成分を利用した散布図
toy_pca |>
    autoplot(loadings = TRUE,
             x = 2, y = 3)          # 主成分の指定
#' パラメタ s を変更 (既定値は1)
toy_pca |>
    autoplot(scale = 0,             # パラメタ s の指定
             loadings = TRUE)
```
:::

::: column-margin
`base R` には関数 `biplot()` が用意されており， 同様の描画を行うことができる．

``` r
#' 第1と第2主成分を利用した散布図
biplot(toy_pca)
#' 第2と第3主成分を利用した散布図
biplot(toy_pca, choices = c(2,3))
#' パラメタ s を変更 (既定値は1)
biplot(toy_pca, scale=0)
```
:::

### 問題

データセット

-   `japan_social.csv`
-   `MASS::UScereal`

の主成分分析の結果を利用してバイプロットによる可視化を行いなさい．

-   標準化したデータでの主成分分析を行いなさい
-   第1主成分と第2主成分でのバイプロットを描きなさい
-   第2主成分と第3主成分でのバイプロットを描きなさい

### 解答欄

```{r}

```

### 解答例

#### `japan_social2.csv` の場合

簡素なバイプロット表示の指定は以下のとおりである．

```{r}
js_pca |>
    autoplot(data = js_data,           # 既定値では第1 vs 第2主成分
             asp = 1,                  # 縦横比を1
             label = TRUE,             # ラベルの表示
             label.label = "都道府県名", # 都道府県名をラベルとする
             loadings = TRUE,          # 主成分負荷の表示
             loadings.label = TRUE)    # 変数名の表示
```

指定可能なオプションの例は以下のようになる．

```{r}
js_pca |>
    autoplot(data = js_data,                  # 地方区分などの補助情報を渡す
             asp = 1,
             colour = "地方区分",               # 地方区分ごとに色付けする
             shape = 19,                      # 以下データ点の修飾．点の形
             size = 1,                        # 点の大きさ
             label = TRUE,                    # 以下ラベルの設定
             label.label = "都道府県名",        # ラベルの指定
             label.repel = TRUE,              # 表示が重ならないように調整
             label.size = 3,                  # 文字の大きさ
             loadings = TRUE,                 # 以下主成分負荷の修飾
             loadings.colour = "orange",      # 色
             loadings.label = TRUE,           # 負荷量のラベル
             loadings.label.repel = TRUE,     # 表示が重ならないように調整
             loadings.label.colour = "brown", # 負荷量のラベルの色
             loadings.label.size = 4)         # 負荷量のラベルの大きさ
```

第1主成分方向の正の向きには大都市をもつ県が集中している． また，人口割合, 商品販売額および森林面積割合は，１人当たり農業産出額とほぼ直交しており, 両者に関連はあまりないといえそう． 第2主成分方向の正の向きには１人当たり農業産出額の上位県が集中している．

気になる項目をいくつか見てみる．

農業産出額を降順に並べる．

```{r}
js_data |>
    arrange(desc(農業算出額)) |> # 指定の列の降順に並べ替える
    slice_head(n = 10) |>      # 上位10件を選択
    gt()
```

第2,3主成分を確認すると，第3主成分方向の負の向きには土地生産性の上位県が集中している．

```{r}
js_pca |>
    autoplot(x = 2, y = 3,
             data = js_data,
             asp = 1,
             colour = "地方区分",
             label = TRUE,
             label.label = "都道府県名",
             label.repel = TRUE,
             loadings = TRUE,
             loadings.label = TRUE,
             loadings.label.repel = TRUE)
```

土地生産性を昇順に並べると，北海道の土地生産性は低いことがわかる．

```{r}
js_data |>
    arrange(土地生産性) |>  # arrangeの既定値は昇順
    slice_head(n = 10) |> # 昇順なので下位10件を選択
    gt()
```

#### `MASS::UScereal` の場合

```{r}
uc_pca |>
    autoplot(data = uc_data,
             asp = 1,
             colour = "mfr", # メーカー毎に色付け
             shape = 19,
             size = 1,
             label = TRUE,
             label.label = "product",
             label.repel = TRUE,
             loadings = TRUE,
             loadings.colour = "orange",
             loadings.label = TRUE,
             loadings.label.repel = TRUE, 
             loadings.label.colour = "brown",
             loadings.label.size = 4) +
    labs(colour = "company")
uc_pca |>
    autoplot(x = 2, y = 3, # 第2 vs 第3
             data = uc_data,
             asp = 1,
             colour = "mfr", 
             shape = 19,
             size = 1, 
             label = TRUE,
             label.label = "product",
             label.repel = TRUE,
             loadings = TRUE,
             loadings.colour = "orange",
             loadings.label = TRUE,
             loadings.label.repel = TRUE, 
             loadings.label.colour = "brown",
             loadings.label.size = 4) +
    labs(colour = "company")
```

第1,2主成分得点で散布図を描き，上と比較してみる．

```{r}
uc_pca |>
    augment(data = uc_data) |>
    ggplot(aes(x = .fittedPC1, y = .fittedPC2,
               label = product, colour = mfr)) +
    geom_point(shape = 19, size = 1) +
    geom_text_repel(size = 3, max.overlaps = 40) +
    theme(aspect.ratio=1) # 縦横比を1
```

配置は変わらないが，座標軸が異なることがわかる． 関数 `autoplot()` において `scale=0` とするとデータの座標は主成分得点となる．

```{r}
options(ggrepel.max.overlaps = 40)
uc_pca |>
    autoplot(scale = 0, 
             asp = 1,
             data = uc_data,
             colour = "mfr", 
             shape = 19,
             size = 1, 
             label = TRUE,
             label.label = "product",
             label.repel = TRUE,
             label.size = 3,
             loadings = TRUE,
             loadings.colour = "orange",
             loadings.label = TRUE,
             loadings.label.repel = TRUE, 
             loadings.label.colour = "brown",
             loadings.label.size = 4) +
    theme(legend.position = "none") # 凡例を表示しない
```