---
title: "統計データ解析"
subtitle: "第6講 サンプルコード"
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)  # 分析結果の描画
```

## 人工データによる主成分分析

### 問題

数値実験により主成分分析の考え方を確認しなさい．

-   以下のモデルに従う人工データを生成する

    ``` r
    #' 観測データ (2次元) の作成 (aのスカラー倍に正規乱数を重畳)
    a <- c(1, 2)/sqrt(5) # 主成分負荷量 (単位ベクトル)
    n <- 100 # データ数
    toy_data <- # aのスカラー倍に正規乱数を重畳(行列として作成)
        runif(n,-1,1) %o% a + rnorm(2*n, sd=0.3) 
    colnames(toy_data) <- paste0("x", 1:2) # 行列に列名(x1,x2)を付与
    toy_data <- as_tibble(toy_data) 
    ```

-   観測データの散布図を作成

-   観測データから第1主成分負荷量を推定

    ``` r
    prcomp(toy_data) # 全ての主成分を計算する
    a_hat <- prcomp(toy_data)$rotation[,1] # 負荷量(rotation)の1列目が第1主成分
    ```

-   散布図上に主成分負荷量を描画

    ``` r
    geom_abline(slope = 傾き, intercept = 切片) # 指定の直線を追加できる  
    ```

### 解答例

2次元の人工データを生成する．

```{r}
set.seed(123123)
n <- 100 # データ数
a <- c(1, 2)/sqrt(5) # 主成分負荷量(単位ベクトル)の設定 (適宜変更せよ)
toy_data <- # aのスカラー倍に正規乱数を重畳(行列として作成)
    runif(n,-1,1) %o% a + rnorm(2*n, sd=0.3) 
colnames(toy_data) <- paste0("x", 1:2) # 行列に列名(x1,x2)を付与
toy_data <- as_tibble(toy_data) 
```

::: callout-tip
以下のようにして列名を書き直しても良い．

``` r
toy_data <-
    as_tibble(runif(n, -1, 1) %o% a + rnorm(2*n, sd = 0.3),
              .name_repair = "minimal") |>
    set_names(paste0("x",1:2))
```
:::

```{r}
p <-
    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(p)
```

$a$ 方向 $(1,2)$ に本質的な情報が集約されていることがわかる．

主成分負荷量の推定を行う．

```{r}
toy_pca <- prcomp(toy_data)
a_hat <- toy_pca$rotation[,1] # 第1主成分
```

第1主成分負荷量が $a$ に非常に近いことが確認できる． 乱数によっては符号が反対になることもある． 前の散布図上に推定された方向を描画することで，両者が近いことが視覚的にも確認できる．

```{r}
p +
    geom_abline(slope = a_hat[2]/a_hat[1], intercept = 0,
                colour = "orange", linetype = "dashed")
```

第1主成分得点を取得し，各データの第1主成分ベクトルを計算する．

```{r}
pc1 <- predict(toy_pca)[,1] 
toy_pc1 <- pc1 %o% a_hat
colnames(toy_pc1) <- paste0("x", 1:2)
toy_pc1 <- as_tibble(toy_pc1)
```

第1主成分ベクトルを元の散布図上で図示する．

```{r}
p +
    geom_point(data = toy_pc1, 
               aes(x = x1, y = x2),
               colour = "purple", shape = 18)
```

## 第1主成分の求め方

### 問題

第1主成分とGram行列の固有ベクトルの関係を調べなさい．

-   人工データを生成する

-   主成分分析を実行する

-   Gram 行列を計算し固有値・固有ベクトルを求める

    ``` r
    #' 中心化を行う
    X <- scale(toy_data, scale = FALSE)
    #' 詳細は '?base::scale' を参照
    #' Gram 行列を計算する
    G <- crossprod(X)
    #' 固有値・固有ベクトルを求める
    eigen(G) # 返り値 'values, vectors' を確認
    #' 詳細は '?base::eigen' を参照
    ```

### 解答例

3次元の人工データを作成する．

```{r}
set.seed(242424)
n <- 50 # データ数
d <- 3
a <- rnorm(d)
(a <- a/sqrt(sum(a^2))) # 主成分負荷量(単位ベクトル)を生成
toy_data <- runif(n, -1, 1) %o% a + rnorm(d*n, sd = 0.1)
colnames(toy_data) <- paste0("x", 1:d) # 行列に列名(x1,...,xd)を付与
toy_data <- as_tibble(toy_data)
```

散布図行列を用いて作成したデータの視覚化を行う．

```{r}
GGally::ggpairs(toy_data) 
```

推定された第1主成分負荷量を取得する．

```{r}
toy_pca <- prcomp(toy_data)
pc1 <- predict(toy_pca)[,1] 
```

::: callout-tip
以下の図は3次元のときのみ実行可能である．

```{r}
#| column: margin
s3d <- scatterplot3d::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)) # 固有値と主成分の標準偏差の関係
```

## 社会生活統計指標の分析

### 問題

都道府県別の社会生活統計指標のデータを用いて主成分分析を行いなさい．

-   データを読み込む

    ``` r
    js_data <- read_csv("data/japan_social.csv") |>
      mutate(Area = as_factor(Area)) # 地方区分を因子化
    ```

-   データの散布図行列を描く

-   各データの箱ひげ図を描き，変数の大きさを確認する

-   主成分負荷量を計算する

    ``` r
    js_pca <- prcomp(js_data[-c(1,7)], scale. = TRUE)
    #' '-c(1,7)' は都道府県名・地方区分を除く．関数 select() を利用することもできる
    #' 'scale.=TRUE' とすると変数を正規化してから解析する
    ```

### 解答例

データの読み込みを行う．

```{r}
js_data <- read_csv("data/japan_social.csv") |>
    mutate(Area = as_factor(Area)) # 地方区分を因子化
```

散布図を用いてデータを視覚化する．

```{r}
js_data |> # 散布図．いくつかの変数は相関強いことがわかる
    select(!c(Pref,Area)) |>  # 都道府県名・地方区分は削除
    GGally::ggpairs() 
```

::: callout-tip
地方ごとに色を変える場合は以下のようにすればよい．

```{r}
#| column: margin
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)))
```
:::

箱ひげ図を用いて変数のばらつきを調べると，変数ごとのばらつきに大きな違いがあることがわかる．

```{r}
js_data |> 
    pivot_longer(where(is.double)) |>  # 実数値(都道府県名・地方区分以外)をまとめる
    mutate(name = as_factor(name)) |>  # 列の順番どおりboxplotを並べる
    ggplot(aes(x = name, y = value)) + # 既定値の name と value を利用
    geom_boxplot(aes(fill = name), show.legend = FALSE) # 変数ごとに色を変える
```

変数ごとに標準化(平均0分散1)したデータフレームを作成するには以下のようにすれば良い．

```{r}
js_data |> 
    mutate(across(where(is.double), \(x)c(scale(x)))) |>
    pivot_longer(where(is.double)) |>  # 実数値をまとめる
    mutate(name = as_factor(name)) |>  # 列の順番どおりboxplotを並べる
    ggplot(aes(x = name, y = value)) + # 既定値の name と value を利用
    geom_boxplot(aes(fill = name), show.legend = FALSE) # 変数ごとに色を変える
```

::: column-margin
データフレームの特定の列のみ変換を行うには関数 `dplyr::across()` を利用する． 関数を組み合わせることで様々な条件が表現できる． 具体的な例は `?dplyr::across` や `?tidyselect::where` を参照のこと．
:::

::: column-margin
関数 `base::scale()` は常に行列を返すため，そのまま用いると期待どおりに動かない． この例ではベクトルを作成する関数 `base::c()` (関数 `base::as.vector()` でも可)と無名関数を用いてこれを解消している． 無名関数の代わりにラムダ式 `~c(scale(.))` などを用いることもできる．
:::

::: callout-tip
箱ひげ図で判り難い場合には密度を表示しても良い．

```{r}
#| column: margin
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_violin(aes(fill = name), show.legend = FALSE) 
```
:::

主成分分析を行う． 変数の標準化は関数 `stats::prcomp()` にオプションを指定することで実行できるので， 元のデータをそのまま渡せば良い．

```{r}
js_pca <- js_data |>
    select(where(is.double)) |> # 実数値の列(都道府県名・地方区分)を抽出
    prcomp(scale. = TRUE) # 変数のばらつきを規格化
js_pca
```

主成分負荷量を表にするには，例えば以下のようにすればよい．

```{r}
js_pca |>
    tidy("v") |> # "v", "rotation", "loadings" または "variables"
    pivot_wider(names_from = PC,        # 横長の形式に変換する
                names_prefix = "PC",    # PC列の番号で新しい列を作る
                values_from = value) |> 
    gt() |>
    fmt_number(decimals = 3)
```

主成分方向から読み取れることは

1.  第1: 人の多さに関する成分(正の向きほど人が多い)
2.  第2: 農業生産力に関する成分(正の向きほど高い)

第1，第2主成分得点を利用して2次元の地図を作成することができる．

```{r}
augment(js_pca, 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)
```

寄与率を表示する．詳しくは次週説明する．

```{r}
summary(js_pca) 
```

::: callout-tip
主成分得点による地図は関数 `ggfortify::autoplot()`を利用することもできる．

```{r}
#| column: margin
autoplot(js_pca,
         scale = 0,            # バイプロットの設定 (主成分得点での表示)
         data = js_data,       # 必要な情報を含むデータフレーム
         colour = "Area",      # 地方区分ごとに色付
         label = TRUE,         # ラベルを付加
         label.label = "Pref", # 都道府県名を追加
         label.repel = TRUE)   # 横にずらして表示
```
:::

2変数での解析例 (AgriとLandを取り上げる，その他の組み合わせでも試みよ) は以下のようになる. 多変数での視覚化は次週詳しく説明する．

```{r}
js_pca2 <- js_data |>
    select(c(Agri,Land)) |>
    prcomp(scale. = TRUE)
a_hat <- js_pca2$rotation[,1] # 主成分負荷量のベクトルを取得
js_pc1 <- predict(js_pca2)[,1] %o% a_hat # 第1主成分を行列・ベクトルで計算
colnames(js_pc1) <- paste0("x", 1:2) # 列名の付与
js_pc1 <- as_tibble(js_pc1) # データフレーム化
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主成分を元の散布図上で図示
               aes(x = x1, y = x2),
               colour = "purple", shape = 18) +
    coord_fixed() # 縦横比を1に指定
```
