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}第7講 主成分分析
評価と視覚化
準備
以下で利用する共通パッケージを読み込む．
主成分分析に用いる主な関数
分析に用いる主な関数としては以下がある．
stats::prcomp(): 主成分分析の実行base::summary(): 主成分負荷・寄与率を表示broom::tidy(): 主成分得点・負荷・寄与率を tibble クラスで取得broom::augment(): tibble クラスで主成分得点を計算ggfortify::autoplot(): 主成分得点による散布図およびバイプロット
#' データフレーム '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) jpamenity.csv を用いた例
人口関連データを整理する．
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(地方区分)) # 地方区分を出現順に因子化散布図により変数間の関係を確認する．
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)) 
主成分分析を行い，バイプロット(主成分得点の散布図と主成分方向の表示)を描画する．
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つのデータセットを使用する
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
データの読み込み方の例．
js_data <-
read_csv("data/japan_social2.csv",
col_types = list(地方区分 = col_factor())) # 地方区分を因子化列の type を指定したい場所が判っていて，先頭の方であれば以下でも可．
#' 1列目は自動判定(?)，2列目は因子(f)，3列目以降は自動判定
js_data <-
read_csv("data/japan_social2.csv", col_types = "?f")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 クラスに変換して整理しておく．
glimpse(MASS::UScereal) # 各変数の属性を確認する．
uc_data <-
MASS::UScereal |>
rownames_to_column(var = "product") |> # 行名の製品名を product 列に加える
as_tibble() # tibble クラスに変換寄与率・累積寄与率
問題
データセット
japan_social.csvMASS::UScereal
について，標準化の有無の違いで寄与率・累積寄与率がどのように異なるか確認しなさい．
toy_data |> prcomp() # 標準化を行わない場合
toy_data |> prcomp(scale. = TRUE) # 標準化を行う場合正式なオプション名は scale. であるが，sc = TRUE など他のオプションと区別できれば短縮表記も可能．
解答例
japan_social2.csv の場合
データを読み込む．
js_data <-
read_csv("data/japan_social2.csv",
col_types = list(地方区分 = col_factor())) # 地方区分を因子化データの性質を把握するために視覚化を行う．以下は全体の散布図および特定の2項目の散布図を描いた例である．
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 = 地方区分))) # 地方区分で色分け
いくつかの変数は相関が強いことがわかる．農業算出額 と 森林面積割合 で散布図を確認する．
js_data |>
ggplot(aes(x = 農業算出額, y = 森林面積割合, label = 都道府県名)) +
geom_point(colour = "blue") +
geom_text_repel(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 <- 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 ) # 標準化あり標準化しない場合の分析結果は以下のようになる． まず，寄与率および累積寄与率を確認する．
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
森林面積割合 -0.014203293 0.04817573 -3.546439e-04 -0.997512798 -0.0494515606
農業算出額 0.972868206 0.12084688 -1.971292e-01 -0.007963943 0.0003640023
人口割合 0.002220545 -0.01161045 2.068058e-05 0.048919679 -0.9987327627
土地生産性 0.221650278 -0.24672608 9.432651e-01 -0.015537081 0.0026195423
商品販売額 0.064745228 -0.96024297 -2.671906e-01 -0.047647214 0.0089675713負荷量が偏る傾向があり，各主成分はほぼ1つの変数に対応している．
寄与率および累積寄与率を tibble クラスで取得すれば， さまざまな表やグラフを作成することができる． 寄与率に関する情報は，各主成分ごとに標準偏差，寄与率，累積寄与率の列を持つデータフレームなので， 例えば以下のような表示が可能である．
#' 寄与率に関する表を作成
js_pca_raw |>
tidy("eigenvalues") |> # "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 |
#' 各主成分ごとに寄与率(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") 
主成分負荷量に関する情報は， 主成分と特性量を指定して負荷量を得る型式の長い表なので， 目的に応じて以下のように整理する必要がる．
#' 主成分負荷量に関する表を作成
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) # 表示を制限| column | PC1 | PC2 | PC3 | PC4 | PC5 |
|---|---|---|---|---|---|
| 森林面積割合 | −0.014 | 0.048 | 0.000 | −0.998 | −0.049 |
| 農業算出額 | 0.973 | 0.121 | −0.197 | −0.008 | 0.000 |
| 人口割合 | 0.002 | −0.012 | 0.000 | 0.049 | −0.999 |
| 土地生産性 | 0.222 | −0.247 | 0.943 | −0.016 | 0.003 |
| 商品販売額 | 0.065 | −0.960 | −0.267 | −0.048 | 0.009 |
#' 各特性量が第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 を利用する方法は以下のようになる．
js_pca_raw |>
autoplot(asp =1, # 縦横の比率を1
scale = 0, # 主成分得点のままの散布図
data = js_data, # ラベルの情報を取得するデータ
colour = "地方区分", # 地方ごとに色付け
label = TRUE, # 以下はラベルの付け方を指定
label.label = "都道府県名", # ラベルに用いる列名
label.repel = TRUE) # ラベルの重なりを低減
同様に標準化した場合の分析結果を見てみる．
js_pca |> # 寄与率に関する情報
tidy("eigenvalues") |>
gt() |>
fmt_number(columns = !1, decimals = 3)| PC | std.dev | percent | cumulative |
|---|---|---|---|
| 1 | 1.590 | 0.506 | 0.506 |
| 2 | 1.070 | 0.229 | 0.735 |
| 3 | 0.820 | 0.134 | 0.869 |
| 4 | 0.708 | 0.100 | 0.969 |
| 5 | 0.392 | 0.031 | 1.000 |
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) | column | PC1 | PC2 | PC3 | PC4 | PC5 |
|---|---|---|---|---|---|
| 森林面積割合 | −0.487 | 0.105 | −0.457 | 0.686 | −0.268 |
| 農業算出額 | 0.134 | 0.812 | 0.479 | 0.305 | 0.035 |
| 人口割合 | 0.585 | −0.151 | 0.045 | 0.164 | −0.778 |
| 土地生産性 | 0.355 | 0.485 | −0.742 | −0.290 | 0.069 |
| 商品販売額 | 0.526 | −0.269 | −0.095 | 0.571 | 0.562 |
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)
寄与率を表示するためには関数 stats::screeplot() を利用することもできる．詳細は ?screeplot を参照．
screeplot(js_pca) 
MASS::UScereal の場合
データを整理する．また，分析を行う前に適当な方法で視覚化をすることを推奨する．
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)))
標準化なしの分析は以下のようになる．
uc_pca_raw <-
uc_data |>
select(where(is.double)) |>
prcomp()
uc_pca_raw |>
tidy("eigenvalues") |>
gt() |>
fmt_number(columns = !1, decimals = 3)| PC | std.dev | percent | cumulative |
|---|---|---|---|
| 1 | 203.185 | 0.770 | 0.770 |
| 2 | 98.693 | 0.182 | 0.952 |
| 3 | 50.362 | 0.047 | 0.999 |
| 4 | 6.824 | 0.001 | 1.000 |
| 5 | 2.372 | 0.000 | 1.000 |
| 6 | 1.490 | 0.000 | 1.000 |
| 7 | 0.947 | 0.000 | 1.000 |
| 8 | 0.550 | 0.000 | 1.000 |
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) | column | PC1 | PC2 | PC3 | PC4 | PC5 | PC6 | PC7 | PC8 |
|---|---|---|---|---|---|---|---|---|
| calories | −0.179 | −0.146 | −0.965 | 0.045 | −0.093 | −0.012 | −0.004 | 0.070 |
| protein | −0.011 | 0.003 | −0.015 | −0.072 | −0.170 | −0.069 | 0.905 | −0.376 |
| fat | −0.003 | −0.001 | −0.015 | 0.077 | −0.396 | 0.371 | −0.354 | −0.758 |
| sodium | −0.493 | −0.842 | 0.220 | 0.008 | −0.002 | −0.001 | 0.000 | −0.001 |
| fibre | −0.027 | 0.020 | 0.010 | −0.076 | −0.187 | −0.919 | −0.218 | −0.257 |
| carbo | −0.015 | −0.027 | −0.110 | −0.693 | 0.621 | 0.042 | −0.077 | −0.336 |
| sugars | −0.009 | −0.001 | −0.046 | 0.707 | 0.620 | −0.104 | 0.032 | −0.318 |
| potassium | −0.851 | 0.519 | 0.078 | −0.006 | 0.013 | 0.033 | −0.002 | 0.011 |
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")
標準化ありの分析は以下のようになる．
uc_pca <-
uc_data |>
select(where(is.double)) |>
prcomp(scale. = TRUE)
uc_pca |>
tidy("eigenvalues") |>
gt() |>
fmt_number(columns = !1, decimals = 3)| PC | std.dev | percent | cumulative |
|---|---|---|---|
| 1 | 2.059 | 0.530 | 0.530 |
| 2 | 1.159 | 0.168 | 0.698 |
| 3 | 1.085 | 0.147 | 0.845 |
| 4 | 0.779 | 0.076 | 0.921 |
| 5 | 0.700 | 0.061 | 0.983 |
| 6 | 0.322 | 0.013 | 0.995 |
| 7 | 0.171 | 0.004 | 0.999 |
| 8 | 0.084 | 0.001 | 1.000 |
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) | column | PC1 | PC2 | PC3 | PC4 | PC5 | PC6 | PC7 | PC8 |
|---|---|---|---|---|---|---|---|---|
| calories | −0.410 | −0.401 | 0.196 | 0.007 | −0.203 | 0.053 | −0.350 | 0.683 |
| protein | −0.448 | 0.176 | 0.072 | −0.195 | −0.128 | −0.830 | 0.032 | −0.142 |
| fat | −0.268 | −0.444 | −0.270 | −0.602 | 0.497 | 0.123 | 0.051 | −0.173 |
| sodium | −0.347 | 0.068 | 0.120 | 0.617 | 0.691 | −0.037 | −0.035 | −0.019 |
| fibre | −0.383 | 0.467 | −0.198 | −0.097 | −0.137 | 0.356 | −0.604 | −0.276 |
| carbo | −0.288 | −0.205 | 0.686 | −0.003 | −0.226 | 0.295 | 0.235 | −0.459 |
| sugars | −0.194 | −0.463 | −0.552 | 0.456 | −0.370 | −0.037 | 0.074 | −0.304 |
| potassium | −0.415 | 0.365 | −0.233 | −0.045 | −0.105 | 0.278 | 0.668 | 0.322 |
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() にオプションを指定することで実行できる．
バイプロット表示
#' データフレームを分析する
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)base R には関数 biplot() が用意されており， 同様の描画を行うことができる．
#' 第1と第2主成分を利用した散布図
biplot(toy_pca)
#' 第2と第3主成分を利用した散布図
biplot(toy_pca, choices = c(2,3))
#' パラメタ s を変更 (既定値は1)
biplot(toy_pca, scale=0)問題
データセット
japan_social.csvMASS::UScereal
の主成分分析の結果を利用してバイプロットによる可視化を行いなさい．
- 標準化したデータでの主成分分析を行いなさい
- 第1主成分と第2主成分でのバイプロットを描きなさい
- 第2主成分と第3主成分でのバイプロットを描きなさい
解答欄
解答例
japan_social2.csv の場合
簡素なバイプロット表示の指定は以下のとおりである．
js_pca |>
autoplot(data = js_data, # 既定値では第1 vs 第2主成分
asp = 1, # 縦横比を1
label = TRUE, # ラベルの表示
label.label = "都道府県名", # 都道府県名をラベルとする
loadings = TRUE, # 主成分負荷の表示
loadings.label = TRUE) # 変数名の表示
指定可能なオプションの例は以下のようになる．
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主成分方向の正の向きには１人当たり農業産出額の上位県が集中している．
気になる項目をいくつか見てみる．
農業産出額を降順に並べる．
js_data |>
arrange(desc(農業算出額)) |> # 指定の列の降順に並べ替える
slice_head(n = 10) |> # 上位10件を選択
gt()| 都道府県名 | 地方区分 | 森林面積割合 | 農業算出額 | 人口割合 | 土地生産性 | 商品販売額 |
|---|---|---|---|---|---|---|
| 北海道 | 北海道 | 67.9 | 1150.6 | 4.23 | 96.8 | 283.3 |
| 宮崎県 | 九州 | 75.8 | 739.1 | 0.87 | 487.7 | 170.6 |
| 鹿児島県 | 九州 | 63.4 | 736.5 | 1.30 | 351.2 | 169.4 |
| 千葉県 | 関東 | 30.4 | 565.5 | 4.90 | 326.1 | 219.7 |
| 群馬県 | 関東 | 63.8 | 530.6 | 1.55 | 321.6 | 270.0 |
| 茨城県 | 関東 | 31.0 | 479.0 | 2.30 | 249.1 | 204.9 |
| 愛知県 | 中部 | 42.2 | 472.3 | 5.89 | 388.9 | 446.9 |
| 佐賀県 | 九州 | 45.2 | 468.7 | 0.66 | 230.3 | 137.9 |
| 熊本県 | 九州 | 60.4 | 456.6 | 1.41 | 285.5 | 172.5 |
| 沖縄県 | 九州 | 46.1 | 452.4 | 1.13 | 232.8 | 145.4 |
第2,3主成分を確認すると，第3主成分方向の負の向きには土地生産性の上位県が集中している．
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)
土地生産性を昇順に並べると，北海道の土地生産性は低いことがわかる．
js_data |>
arrange(土地生産性) |> # arrangeの既定値は昇順
slice_head(n = 10) |> # 昇順なので下位10件を選択
gt()| 都道府県名 | 地方区分 | 森林面積割合 | 農業算出額 | 人口割合 | 土地生産性 | 商品販売額 |
|---|---|---|---|---|---|---|
| 北海道 | 北海道 | 67.9 | 1150.6 | 4.23 | 96.8 | 283.3 |
| 秋田県 | 東北 | 70.5 | 268.7 | 0.81 | 98.5 | 153.3 |
| 富山県 | 中部 | 56.6 | 276.1 | 0.84 | 98.5 | 192.4 |
| 福井県 | 中部 | 73.9 | 216.1 | 0.62 | 98.5 | 167.3 |
| 滋賀県 | 関西 | 50.5 | 222.8 | 1.11 | 104.9 | 170.7 |
| 石川県 | 中部 | 66.0 | 271.3 | 0.91 | 112.0 | 222.9 |
| 宮城県 | 東北 | 55.9 | 299.9 | 1.84 | 125.3 | 365.9 |
| 山口県 | 中国 | 71.6 | 216.9 | 1.11 | 125.8 | 158.9 |
| 福島県 | 東北 | 67.9 | 236.4 | 1.51 | 127.1 | 184.5 |
| 島根県 | 中国 | 77.5 | 214.1 | 0.55 | 140.8 | 141.1 |
MASS::UScereal の場合
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主成分得点で散布図を描き，上と比較してみる．
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 とするとデータの座標は主成分得点となる．
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") # 凡例を表示しない