---
title: "統計データ解析"
subtitle: "第7講 サンプルコード"
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
editor_options: 
  chunk_output_type: console
---

## 準備

以下で利用する共通パッケージを読み込む．

```{r}
library(conflicted) # 関数名の衝突を警告
conflicts_prefer(   # 優先的に使う関数を指定
  dplyr::filter(),
  dplyr::select(),
  dplyr::lag(),
  )
library(tidyverse)  
library(broom)      # 解析結果を tibble 形式に集約
library(gt)         # 表の作成
library(gtsummary)  # 分析結果の表の作成
library(ggfortify)  # biplot表示のため
```

## 寄与率・累積寄与率

### 問題

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

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

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

-   `japan_social.csv` の読み込み方の例

    ``` r
    js_data <- read_csv("data/japan_social.csv") |>
        mutate(Area = as_factor(Area))
    ```

-   `MASS::UScereal` の整理の例

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

    ::: column-margin
    `base R` の `data.frame` 型なので tibble 型に変換しておく．
    :::

### 解答例

#### `japan_social.csv` の場合

総務省統計局の都道府県別の社会生活統計指標データの各列は以下のとおりである．

-   Pref: 都道府県名
-   Forest: 森林面積割合 (%) 2014年
-   Agri: 就業者１人当たり農業産出額(販売農家）(万円) 2014年
-   Ratio: 全国総人口に占める人口割合 (%) 2015年
-   Land: 土地生産性（耕地面積１ヘクタール当たり）(万円) 2014年
-   Goods: 商業年間商品販売額［卸売業＋小売業］（事業所当たり）(百万円) 2013年
-   Area: 地方区分

データを読み込む． このとき色分けのため，地方区分は因子化しておく．

```{r}
js_data <- read_csv("data/japan_social.csv") |>
    mutate(Area = as_factor(Area))
```

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

```{r}
js_data |> # 散布図．いくつかの変数は相関強いことがわかる
    GGally::ggpairs(columns = 2:6, # 都道府県名・地方区分は除く
                    legend = c(2,1), # 2行1列のグラフから凡例を作成
                    lower = list(continuous = GGally::wrap("points", alpha = .5),
                                 mapping = aes(colour = Area)))
js_data |>
    ggplot(aes(x = Agri, y = Forest)) +
    geom_point(colour = "blue") +
    geom_text(aes(label = Pref),
              colour = "brown",
              vjust = 1, nudge_y = 0.0, check_overlap = TRUE)
```

::: callout-tip
点と文字が被らないように座標をずらし，密集したところは消している． `package::ggrepel` を利用すれば自動調整してくれる．

```{r}
#| column: margin
js_data |>
    ggplot(aes(x = Agri, y = Forest)) +
    geom_point(colour = "blue") +
    ggrepel::geom_text_repel(aes(label = Pref), 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
js_pca_raw <- prcomp(js_data[-c(1,7)])            # 標準化なし
js_pca <- prcomp(js_data[-c(1,7)], scale. = TRUE) # 標準化あり
```

各行に県名 `Pref` をつけるには以下のようにすれば良い．

``` r
js_pca_raw <- js_data |>
    column_to_rownames(var = "Pref") |> # 'Pref'を行名に変換
    select(where(is.double)) |> 
    prcomp()                            # 標準化なし
js_pca <- js_data |> 
    column_to_rownames(var = "Pref") |> 
    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("d") |> # "d" または "eigenvalues" または "pcs"
    gt() |>      
    fmt_number(columns = !1, decimals = 3) # 1列目以外小数点以下3桁
js_pca_raw |>
    tidy("d") |> 
    ggplot(aes(x = PC, y = percent)) + # 各主成分(PC)ごとに寄与率(percent)を表示
    geom_bar(stat = "identity")
js_pca_raw |>
    tidy("d") |> 
    ggplot(aes(x = PC, y = cumulative)) + # 各主成分(PC)ごとに累積寄与率(cumlative)を表示
    geom_bar(stat = "identity")
```

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

```{r}
autoplot(js_pca_raw, scale = 0,
         data = js_data,                     # ラベルの情報を取得するデータ
         label = TRUE, label.label = "Pref") # ラベルの付け方を指定
```

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

```{r}
summary(js_pca)
js_pca$rotation
js_pca |>
    tidy("d") |> 
    ggplot(aes(x = PC, y = percent)) + # 寄与率(percent)を表示
    geom_bar(stat = "identity")
js_pca |>
    tidy("d") |> 
    ggplot(aes(x = PC, y = cumulative)) + # 累積寄与率(cumlative)を表示
    geom_bar(stat = "identity")
autoplot(js_pca, scale = 0,
         data = js_data, label = TRUE, label.label = "Pref")
```

::: callout-tip
寄与率を表示するためには関数stats::screeplot()を利用してもよい． 詳細は `?screeplot` を参照のこと．

```{r}
#| column: margin
screeplot(js_pca) 
```
:::

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

元のデータは data.frame 形式なので，tibble に変換して整理しておく． また，適当な方法で視覚化をすることを推奨する．

```{r}
uc_data <- MASS::UScereal |> 
    rownames_to_column(var = "product") |> as_tibble() 
```

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

```{r}
uc_pca_raw <- uc_data |>
    select(where(is.double)) |>
    prcomp() 
summary(uc_pca_raw)
uc_pca_raw$rotation
uc_pca_raw |>
    tidy("d") |> 
    ggplot(aes(x = PC, y = percent)) + geom_bar(stat = "identity")
uc_pca_raw |>
    tidy("d") |> 
    ggplot(aes(x = PC, y = cumulative)) + geom_bar(stat = "identity")
autoplot(uc_pca_raw, scale = 0,
         data = uc_data, label = TRUE, label.label = "product") 
```

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

```{r}
uc_pca <- uc_data |>
    select(where(is.double)) |>
    prcomp(scale. = TRUE)
summary(uc_pca)
uc_pca$rotation
uc_pca |>
    tidy("d") |> 
    ggplot(aes(x = PC, y = percent)) + geom_bar(stat = "identity")
uc_pca |>
    tidy("d") |> 
    ggplot(aes(x = PC, y = cumulative)) + geom_bar(stat = "identity")
autoplot(uc_pca_raw, scale = 0,
         data = uc_data, label = TRUE, label.label = "product") 
```

## 主成分分析の視覚化

### 問題

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

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

``` r
autoplot(prcompの結果, x = x軸成分, y = y軸成分) # 主成分の指定
```

### 解答例

#### `japan_social.csv` の場合

簡素な biplot の指定は以下のとおりである．

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

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

```{r}
autoplot(js_pca, 
         data = js_data,                  # 地方区分などの補助情報を渡す
         colour = "Area",                 # 地方区分ごとに色付けする
         shape = 19,                      # 以下データ点の修飾．点の形
         size = 1,                        # 点の大きさ
         label = TRUE,                    # 以下ラベルの設定
         label.label = "Pref",            # ラベルの指定
         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(Agri) |> gt()
```

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

```{r}
autoplot(js_pca, x = 2, y = 3,
         data = js_data, colour = "Area",
         label = TRUE, label.label = "Pref", label.repel = TRUE,
         loadings = TRUE, loadings.label = TRUE, loadings.label.repel = TRUE)
```

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

```{r}
js_data |> arrange(desc(Land)) |> gt()
```

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

```{r}
autoplot(uc_pca,
         data = uc_data, 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) +
    theme(legend.position = "none")
autoplot(uc_pca, x = 2, y = 3, # 第2 vs 第3
         data = uc_data, 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) +
    theme(legend.position = "none")
```

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

```{r}
augment(uc_pca, data = uc_data) |>
    ggplot(aes(x = .fittedPC1, y = .fittedPC2, label = product, colour = mfr)) +
    geom_point(shape = 19, size = 1) +
    ggrepel::geom_text_repel(size = 3, max.overlaps = 40) +
    theme(legend.position = "none")
```

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

```{r}
options(ggrepel.max.overlaps = 40)
autoplot(uc_pca, scale = 0, 
         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")
```