library(conflicted) # 関数名の衝突を警告
conflicts_prefer( # 優先的に使う関数を指定
dplyr::filter(),
dplyr::select(),
dplyr::lag(),
)
library(tidyverse)
library(broom) # 解析結果を tibble 形式に集約
library(gt) # 表の作成
library(gtsummary) # 分析結果の表の作成
library(ggfortify) # 分析結果の描画統計データ解析
第6講 サンプルコード
準備
以下で利用する共通パッケージを読み込む．
人工データによる主成分分析
問題
数値実験により主成分分析の考え方を確認しなさい．
以下のモデルに従う人工データを生成する
#' 観測データ (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主成分負荷量を推定
prcomp(toy_data) # 全ての主成分を計算する a_hat <- prcomp(toy_data)$rotation[,1] # 負荷量(rotation)の1列目が第1主成分散布図上に主成分負荷量を描画
geom_abline(slope = 傾き, intercept = 切片) # 指定の直線を追加できる
解答例
2次元の人工データを生成する．
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)
Tip
以下のようにして列名を書き直しても良い．
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))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) に本質的な情報が集約されていることがわかる．
主成分負荷量の推定を行う．
toy_pca <- prcomp(toy_data)
a_hat <- toy_pca$rotation[,1] # 第1主成分第1主成分負荷量が a に非常に近いことが確認できる． 乱数によっては符号が反対になることもある． 前の散布図上に推定された方向を描画することで，両者が近いことが視覚的にも確認できる．
p +
geom_abline(slope = a_hat[2]/a_hat[1], intercept = 0,
colour = "orange", linetype = "dashed")
第1主成分得点を取得し，各データの第1主成分ベクトルを計算する．
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主成分ベクトルを元の散布図上で図示する．
p +
geom_point(data = toy_pc1,
aes(x = x1, y = x2),
colour = "purple", shape = 18)
第1主成分の求め方
問題
第1主成分とGram行列の固有ベクトルの関係を調べなさい．
人工データを生成する
主成分分析を実行する
Gram 行列を計算し固有値・固有ベクトルを求める
#' 中心化を行う X <- scale(toy_data, scale = FALSE) #' 詳細は '?base::scale' を参照 #' Gram 行列を計算する G <- crossprod(X) #' 固有値・固有ベクトルを求める eigen(G) # 返り値 'values, vectors' を確認 #' 詳細は '?base::eigen' を参照
解答例
3次元の人工データを作成する．
set.seed(242424)
n <- 50 # データ数
d <- 3
a <- rnorm(d)
(a <- a/sqrt(sum(a^2))) # 主成分負荷量(単位ベクトル)を生成[1] 0.1370560 -0.4962935 0.8572680toy_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)散布図行列を用いて作成したデータの視覚化を行う．
GGally::ggpairs(toy_data) 
推定された第1主成分負荷量を取得する．
toy_pca <- prcomp(toy_data)
pc1 <- predict(toy_pca)[,1]
Tip
以下の図は3次元のときのみ実行可能である．
s3d <- scatterplot3d::scatterplot3d(toy_data, type="h",
## asp=1, # 軸の比率を揃える場合
highlight.3d=TRUE)
s3d$points3d(pc1 %o% a, col="blue")
主成分負荷量の推定を固有値分解と比較する．
toy_eigen <- # 固有値分解
eigen(crossprod(scale(toy_data, scale=FALSE)))
toy_pca$rotation # 主成分負荷量 PC1 PC2 PC3
x1 0.1352560 -0.6540848 0.7442304
x2 -0.5069376 0.5996896 0.6191823
x3 0.8513049 0.4610265 0.2504686toy_eigen$vectors # 固有ベクトル (符号を除いて主成分負荷量と一致) [,1] [,2] [,3]
[1,] 0.1352560 0.6540848 0.7442304
[2,] -0.5069376 -0.5996896 0.6191823
[3,] 0.8513049 -0.4610265 0.2504686toy_pca$sdev # 主成分の標準偏差[1] 0.5056561 0.1256128 0.1142378sqrt(toy_eigen$values/(n-1)) # 固有値と主成分の標準偏差の関係[1] 0.5056561 0.1256128 0.1142378社会生活統計指標の分析
問題
都道府県別の社会生活統計指標のデータを用いて主成分分析を行いなさい．
データを読み込む
js_data <- read_csv("data/japan_social.csv") |> mutate(Area = as_factor(Area)) # 地方区分を因子化データの散布図行列を描く
各データの箱ひげ図を描き，変数の大きさを確認する
主成分負荷量を計算する
js_pca <- prcomp(js_data[-c(1,7)], scale. = TRUE) #' '-c(1,7)' は都道府県名・地方区分を除く．関数 select() を利用することもできる #' 'scale.=TRUE' とすると変数を正規化してから解析する
解答例
データの読み込みを行う．
js_data <- read_csv("data/japan_social.csv") |>
mutate(Area = as_factor(Area)) # 地方区分を因子化散布図を用いてデータを視覚化する．
js_data |> # 散布図．いくつかの変数は相関強いことがわかる
select(!c(Pref,Area)) |> # 都道府県名・地方区分は削除
GGally::ggpairs() 
Tip
地方ごとに色を変える場合は以下のようにすればよい．
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 |>
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)したデータフレームを作成するには以下のようにすれば良い．
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) # 変数ごとに色を変える
データフレームの特定の列のみ変換を行うには関数 dplyr::across() を利用する． 関数を組み合わせることで様々な条件が表現できる． 具体的な例は ?dplyr::across や ?tidyselect::where を参照のこと．
関数 base::scale() は常に行列を返すため，そのまま用いると期待どおりに動かない． この例ではベクトルを作成する関数 base::c() (関数 base::as.vector() でも可)と無名関数を用いてこれを解消している． 無名関数の代わりにラムダ式 ~c(scale(.)) などを用いることもできる．
Tip
箱ひげ図で判り難い場合には密度を表示しても良い．
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() にオプションを指定することで実行できるので， 元のデータをそのまま渡せば良い．
js_pca <- js_data |>
select(where(is.double)) |> # 実数値の列(都道府県名・地方区分)を抽出
prcomp(scale. = TRUE) # 変数のばらつきを規格化
js_pcaStandard deviations (1, .., p=5):
[1] 1.5903931 1.0698965 0.8195653 0.7076020 0.3918975
Rotation (n x k) = (5 x 5):
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("v") |> # "v", "rotation", "loadings" または "variables"
pivot_wider(names_from = PC, # 横長の形式に変換する
names_prefix = "PC", # PC列の番号で新しい列を作る
values_from = value) |>
gt() |>
fmt_number(decimals = 3)| column | PC1 | PC2 | PC3 | PC4 | PC5 |
|---|---|---|---|---|---|
| Forest | −0.487 | 0.105 | −0.457 | 0.686 | −0.268 |
| Agri | 0.134 | 0.812 | 0.479 | 0.305 | 0.035 |
| Ratio | 0.585 | −0.151 | 0.045 | 0.164 | −0.778 |
| Land | 0.355 | 0.485 | −0.742 | −0.290 | 0.069 |
| Goods | 0.526 | −0.269 | −0.095 | 0.571 | 0.562 |
主成分方向から読み取れることは
- 第1: 人の多さに関する成分(正の向きほど人が多い)
- 第2: 農業生産力に関する成分(正の向きほど高い)
第1，第2主成分得点を利用して2次元の地図を作成することができる．
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)
寄与率を表示する．詳しくは次週説明する．
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
Tip
主成分得点による地図は関数 ggfortify::autoplot()を利用することもできる．
autoplot(js_pca,
scale = 0, # バイプロットの設定 (主成分得点での表示)
data = js_data, # 必要な情報を含むデータフレーム
colour = "Area", # 地方区分ごとに色付
label = TRUE, # ラベルを付加
label.label = "Pref", # 都道府県名を追加
label.repel = TRUE) # 横にずらして表示
2変数での解析例 (AgriとLandを取り上げる，その他の組み合わせでも試みよ) は以下のようになる. 多変数での視覚化は次週詳しく説明する．
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に指定