統計データ解析

第7講 サンプルコード

Published

November 15, 2024

準備

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

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

寄与率・累積寄与率

問題

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

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

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

  • japan_social.csv の読み込み方の例

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

  • MASS::UScereal の整理の例

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

base R の data.frame 型なので tibble 型に変換しておく．

解答例

japan_social.csv の場合

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

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

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

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

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

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)

Tip

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

js_data |>
    ggplot(aes(x = Agri, y = Forest)) +
    geom_point(colour = "blue") +
    ggrepel::geom_text_repel(aes(label = Pref), colour = "brown")

主成分分析を実行する．

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

必要な列のみ選択すれば，分析は行うことができる．

js_pca_raw <- prcomp(js_data[-c(1,7)])            # 標準化なし
js_pca <- prcomp(js_data[-c(1,7)], scale. = TRUE) # 標準化あり

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

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)                # 標準化あり

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

summary(js_pca_raw) 
Importance of components:
                           PC1     PC2     PC3      PC4     PC5
Standard deviation     173.275 148.037 81.5231 12.97198 1.05151
Proportion of Variance   0.511   0.373  0.1131  0.00286 0.00002
Cumulative Proportion    0.511   0.884  0.9971  0.99998 1.00000

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

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

js_pca_raw$rotation 
                PC1         PC2           PC3          PC4           PC5
Forest -0.014203293  0.04817573 -3.546439e-04 -0.997512798 -0.0494515606
Agri    0.972868206  0.12084688 -1.971292e-01 -0.007963943  0.0003640023
Ratio   0.002220545 -0.01161045  2.068058e-05  0.048919679 -0.9987327627
Land    0.221650278 -0.24672608  9.432651e-01 -0.015537081  0.0026195423
Goods   0.064745228 -0.96024297 -2.671906e-01 -0.047647214  0.0089675713

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

寄与率および累積寄与率を tibble 形式で取得するには以下のようにすれば良い．

js_pca_raw |>    # 表 を作成
    tidy("d") |> # "d" または "eigenvalues" または "pcs"
    gt() |>      
    fmt_number(columns = !1, decimals = 3) # 1列目以外3桁
PC std.dev percent cumulative
1 173.275 0.511 0.511
2 148.037 0.373 0.884
3 81.523 0.113 0.997
4 12.972 0.003 1.000
5 1.052 0.000 1.000 
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 を利用する方法は以下のようになる．

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

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

summary(js_pca)
Importance of components:
                          PC1    PC2    PC3    PC4     PC5
Standard deviation     1.5904 1.0699 0.8196 0.7076 0.39190
Proportion of Variance 0.5059 0.2289 0.1343 0.1001 0.03072
Cumulative Proportion  0.5059 0.7348 0.8691 0.9693 1.00000
js_pca$rotation
              PC1        PC2         PC3        PC4         PC5
Forest -0.4871498  0.1045813 -0.45748795  0.6859649 -0.26815060
Agri    0.1339190  0.8115056  0.47912767  0.3045447  0.03483694
Ratio   0.5851294 -0.1511042  0.04467249  0.1640953 -0.77837539
Land    0.3547649  0.4851374 -0.74167904 -0.2897485  0.06885892
Goods   0.5258481 -0.2689436 -0.09517368  0.5708093  0.56238052
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")

Tip

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

screeplot(js_pca)

MASS::UScereal の場合

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

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

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

uc_pca_raw <- uc_data |>
    select(where(is.double)) |>
    prcomp() 
summary(uc_pca_raw)
Importance of components:
                          PC1     PC2      PC3     PC4    PC5     PC6     PC7
Standard deviation     203.19 98.6927 50.36206 6.82381 2.3718 1.48975 0.94740
Proportion of Variance   0.77  0.1817  0.04731 0.00087 0.0001 0.00004 0.00002
Cumulative Proportion    0.77  0.9517  0.99896 0.99983 0.9999 0.99998 0.99999
                           PC8
Standard deviation     0.55042
Proportion of Variance 0.00001
Cumulative Proportion  1.00000
uc_pca_raw$rotation
                   PC1           PC2         PC3          PC4          PC5
calories  -0.179207473 -0.1462324241 -0.96480715  0.045234798 -0.092696691
protein   -0.011145612  0.0025160642 -0.01538123 -0.072061722 -0.169798607
fat       -0.002912531 -0.0006278695 -0.01548280  0.077006299 -0.396417270
sodium    -0.493173560 -0.8416992943  0.21967258  0.007757990 -0.001777641
fibre     -0.027380174  0.0202681877  0.01023689 -0.075965904 -0.187070010
carbo     -0.015056930 -0.0272462011 -0.10970338 -0.693413074  0.620733329
sugars    -0.008609442 -0.0013826562 -0.04640333  0.707209145  0.620447068
potassium -0.850577043  0.5186488067  0.07824275 -0.005786195  0.012896751
                   PC6           PC7           PC8
calories  -0.012034965 -0.0040178657  0.0696281364
protein   -0.069385729  0.9054169859 -0.3755184499
fat        0.370507919 -0.3543044317 -0.7575401217
sodium    -0.001439163 -0.0001605839 -0.0008059101
fibre     -0.918767097 -0.2184131271 -0.2571603627
carbo      0.042042740 -0.0770477377 -0.3363935184
sugars    -0.103827276  0.0320770971 -0.3175882779
potassium  0.032892502 -0.0016414154  0.0107594252
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")

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

uc_pca <- uc_data |>
    select(where(is.double)) |>
    prcomp(scale. = TRUE)
summary(uc_pca)
Importance of components:
                          PC1    PC2    PC3     PC4     PC5     PC6     PC7
Standard deviation     2.0595 1.1595 1.0847 0.77916 0.70038 0.32184 0.17080
Proportion of Variance 0.5302 0.1681 0.1471 0.07589 0.06132 0.01295 0.00365
Cumulative Proportion  0.5302 0.6982 0.8453 0.92120 0.98252 0.99547 0.99911
                           PC8
Standard deviation     0.08427
Proportion of Variance 0.00089
Cumulative Proportion  1.00000
uc_pca$rotation
                 PC1         PC2        PC3          PC4        PC5         PC6
calories  -0.4095924 -0.40065656  0.1958834  0.007486492 -0.2032066  0.05255351
protein   -0.4476405  0.17567026  0.0722067 -0.194640140 -0.1277506 -0.82959021
fat       -0.2683607 -0.44417203 -0.2702646 -0.602181956  0.4973648  0.12339688
sodium    -0.3474018  0.06764946  0.1195216  0.616762389  0.6907158 -0.03693360
fibre     -0.3829156  0.46677617 -0.1976082 -0.097098382 -0.1367666  0.35604740
carbo     -0.2882507 -0.20478232  0.6861866 -0.003106577 -0.2260198  0.29493081
sugars    -0.1938910 -0.46258366 -0.5521416  0.455559106 -0.3701522 -0.03653264
potassium -0.4145359  0.36462256 -0.2330796 -0.045479335 -0.1054703  0.27809848
                  PC7         PC8
calories  -0.35011143  0.68311726
protein    0.03187206 -0.14178644
fat        0.05108036 -0.17268441
sodium    -0.03505495 -0.01931262
fibre     -0.60446872 -0.27597329
carbo      0.23510504 -0.45909745
sugars     0.07382281 -0.30369583
potassium  0.66817773  0.32232225
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主成分でのバイプロットを描きなさい．
autoplot(prcompの結果, x = x軸成分, y = y軸成分) # 主成分の指定

解答例

japan_social.csv の場合

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

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

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

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主成分方向の正の向きには１人当たり農業産出額の上位県が集中している．

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

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

js_data |> arrange(Agri) |> gt()
Pref Forest Agri Ratio Land Goods Area
Nara 76.8 207.0 1.07 182.7 147.0 Kansai
Shimane 77.5 214.1 0.55 140.8 141.1 Chugoku
Fukui 73.9 216.1 0.62 98.5 167.3 Chubu
Osaka 30.1 216.3 6.96 238.8 451.2 Kansai
Yamaguchi 71.6 216.9 1.11 125.8 158.9 Chugoku
Shiga 50.5 222.8 1.11 104.9 170.7 Kansai
Fukushima 67.9 236.4 1.51 127.1 184.5 Tohoku
Kagawa 46.4 249.5 0.77 242.9 232.9 Shikoku
Tottori 73.3 249.9 0.45 187.6 162.2 Chugoku
Wakayama 76.4 251.1 0.76 278.4 136.4 Kansai
Okayama 68.0 254.8 1.51 184.9 207.8 Chugoku
Hyogo 66.7 261.2 4.35 197.7 212.5 Kansai
Kyoto 74.2 267.8 2.05 212.5 196.7 Kansai
Tokyo 34.8 268.5 10.63 404.7 1062.6 Kanto
Akita 70.5 268.7 0.81 98.5 153.3 Tohoku
Ishikawa 66.0 271.3 0.91 112.0 222.9 Chubu
Toyama 56.6 276.1 0.84 98.5 192.4 Chubu
Nagano 75.5 280.0 1.65 211.3 194.4 Chubu
Gifu 79.0 283.7 1.60 192.1 167.9 Chubu
Hiroshima 71.8 286.2 2.24 192.2 304.6 Chugoku
Yamanashi 77.8 287.4 0.66 325.3 156.2 Chubu
Ehime 70.3 288.5 1.09 231.6 179.4 Shikoku
Miyagi 55.9 299.9 1.84 125.3 365.9 Tohoku
Niigata 63.5 308.6 1.81 141.9 205.5 Chubu
Mie 64.3 310.6 1.43 174.3 170.1 Kansai
Tokushima 75.2 315.4 0.59 313.5 134.5 Shikoku
Kanagawa 38.8 322.8 7.18 396.4 246.1 Kanto
Saitama 31.9 324.7 5.72 247.0 244.7 Kanto
Iwate 74.9 334.3 1.01 155.2 179.4 Tohoku
Kochi 83.3 354.2 0.57 339.9 137.9 Shikoku
Oita 70.7 360.1 0.92 222.8 148.3 Kyushu
Shizuoka 63.1 375.8 2.91 314.5 211.4 Chubu
Fukuoka 44.5 381.0 4.01 255.6 295.7 Kyushu
Yamagata 68.7 396.3 0.88 174.1 157.5 Tohoku
Tochigi 53.2 402.6 1.55 199.6 204.3 Kanto
Nagasaki 58.4 428.9 1.08 296.0 154.0 Kyushu
Aomori 63.8 444.7 1.03 186.0 183.0 Tohoku
Okinawa 46.1 452.4 1.13 232.8 145.4 Kyushu
Kumamoto 60.4 456.6 1.41 285.5 172.5 Kyushu
Saga 45.2 468.7 0.66 230.3 137.9 Kyushu
Aichi 42.2 472.3 5.89 388.9 446.9 Chubu
Ibaraki 31.0 479.0 2.30 249.1 204.9 Kanto
Gumma 63.8 530.6 1.55 321.6 270.0 Kanto
Chiba 30.4 565.5 4.90 326.1 219.7 Kanto
Kagoshima 63.4 736.5 1.30 351.2 169.4 Kyushu
Miyazaki 75.8 739.1 0.87 487.7 170.6 Kyushu
Hokkaido 67.9 1150.6 4.23 96.8 283.3 Hokkaido

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

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)

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

js_data |> arrange(desc(Land)) |> gt()
Pref Forest Agri Ratio Land Goods Area
Miyazaki 75.8 739.1 0.87 487.7 170.6 Kyushu
Tokyo 34.8 268.5 10.63 404.7 1062.6 Kanto
Kanagawa 38.8 322.8 7.18 396.4 246.1 Kanto
Aichi 42.2 472.3 5.89 388.9 446.9 Chubu
Kagoshima 63.4 736.5 1.30 351.2 169.4 Kyushu
Kochi 83.3 354.2 0.57 339.9 137.9 Shikoku
Chiba 30.4 565.5 4.90 326.1 219.7 Kanto
Yamanashi 77.8 287.4 0.66 325.3 156.2 Chubu
Gumma 63.8 530.6 1.55 321.6 270.0 Kanto
Shizuoka 63.1 375.8 2.91 314.5 211.4 Chubu
Tokushima 75.2 315.4 0.59 313.5 134.5 Shikoku
Nagasaki 58.4 428.9 1.08 296.0 154.0 Kyushu
Kumamoto 60.4 456.6 1.41 285.5 172.5 Kyushu
Wakayama 76.4 251.1 0.76 278.4 136.4 Kansai
Fukuoka 44.5 381.0 4.01 255.6 295.7 Kyushu
Ibaraki 31.0 479.0 2.30 249.1 204.9 Kanto
Saitama 31.9 324.7 5.72 247.0 244.7 Kanto
Kagawa 46.4 249.5 0.77 242.9 232.9 Shikoku
Osaka 30.1 216.3 6.96 238.8 451.2 Kansai
Okinawa 46.1 452.4 1.13 232.8 145.4 Kyushu
Ehime 70.3 288.5 1.09 231.6 179.4 Shikoku
Saga 45.2 468.7 0.66 230.3 137.9 Kyushu
Oita 70.7 360.1 0.92 222.8 148.3 Kyushu
Kyoto 74.2 267.8 2.05 212.5 196.7 Kansai
Nagano 75.5 280.0 1.65 211.3 194.4 Chubu
Tochigi 53.2 402.6 1.55 199.6 204.3 Kanto
Hyogo 66.7 261.2 4.35 197.7 212.5 Kansai
Hiroshima 71.8 286.2 2.24 192.2 304.6 Chugoku
Gifu 79.0 283.7 1.60 192.1 167.9 Chubu
Tottori 73.3 249.9 0.45 187.6 162.2 Chugoku
Aomori 63.8 444.7 1.03 186.0 183.0 Tohoku
Okayama 68.0 254.8 1.51 184.9 207.8 Chugoku
Nara 76.8 207.0 1.07 182.7 147.0 Kansai
Mie 64.3 310.6 1.43 174.3 170.1 Kansai
Yamagata 68.7 396.3 0.88 174.1 157.5 Tohoku
Iwate 74.9 334.3 1.01 155.2 179.4 Tohoku
Niigata 63.5 308.6 1.81 141.9 205.5 Chubu
Shimane 77.5 214.1 0.55 140.8 141.1 Chugoku
Fukushima 67.9 236.4 1.51 127.1 184.5 Tohoku
Yamaguchi 71.6 216.9 1.11 125.8 158.9 Chugoku
Miyagi 55.9 299.9 1.84 125.3 365.9 Tohoku
Ishikawa 66.0 271.3 0.91 112.0 222.9 Chubu
Shiga 50.5 222.8 1.11 104.9 170.7 Kansai
Akita 70.5 268.7 0.81 98.5 153.3 Tohoku
Toyama 56.6 276.1 0.84 98.5 192.4 Chubu
Fukui 73.9 216.1 0.62 98.5 167.3 Chubu
Hokkaido 67.9 1150.6 4.23 96.8 283.3 Hokkaido

MASS::UScereal の場合

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主成分得点で散布図を描き，上と比較してみる．

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 とするとデータの座標は主成分得点となる．

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")