R練習スクリプト - yoshida1983/protc2026 GitHub Wiki

Palmer Penguins: PCA & Clustering 実習

パッケージの準備

# パッケージのインストール（初回のみ、先頭の「#」を消して実行）
# install.packages(c("palmerpenguins", "tidyverse", "factoextra", "pheatmap"))

# パッケージの読み込み
library(palmerpenguins)
library(tidyverse)
library(factoextra)
library(pheatmap)

1) データ準備

# 解析に使用する変数を抽出し、NA（欠損値）を除去
peng <- penguins |>
  select(species,
         bill_length_mm,
         bill_depth_mm,
         flipper_length_mm,
         body_mass_g) |>
  drop_na()

# データの確認
glimpse(peng)

# PCAとクラスタリングに使用する数値データを抽出
X <- peng |>
  select(bill_length_mm,
         bill_depth_mm,
         flipper_length_mm,
         body_mass_g)

2) PCA（主成分分析）

# 各変数を標準化してPCAを実行
pca <- prcomp(X, center = TRUE, scale. = TRUE)

# 各主成分の寄与率を確認
summary(pca)

# スクリープロット（各主成分の寄与率）
fviz_eig(pca)

# PCA散布図（speciesで色分け）
fviz_pca_ind(pca,
             geom.ind = "point",
             habillage = peng$species,
             pointsize = 2) +
  labs(color = "Species",
       shape = "Species")

# 各変数の主成分への寄与を確認
fviz_pca_var(pca, col.var = "contrib")

3) k-means クラスタリング

# PCAのスコア（PC1, PC2）を取り出す
scores <- as.data.frame(pca$x[, 1:2])

# 乱数の初期値を固定し、解析結果を再現できるようにする
set.seed(123)

# クラスタ数 k の当たりをつける
# エルボー法
fviz_nbclust(scores, kmeans,
             method = "wss",
             nstart = 50)

# シルエット法
fviz_nbclust(scores, kmeans,
             method = "silhouette",
             nstart = 50)

# 図を参考に、ここでは k = 3 で実行
set.seed(123)
km <- kmeans(scores, centers = 3, nstart = 50)

# クラスタリング結果を可視化
fviz_cluster(km, data = scores, geom = "point")

# speciesとの対応を確認
tab <- table(Cluster = km$cluster, Species = peng$species)
tab

# 各クラスタで最も多い種を調べる（1 = Adelie, 2 = Chinstrap, 3 = Gentoo）
major <- apply(tab, 1, which.max)
major

# 2つのクラスタで多数派の種が同じ場合は、番号の付け直しが意味を持たない
if (anyDuplicated(major) > 0) {
  warning("同じ種が多数派のクラスタがあります。k や nstart を見直してください。")
}

# 多数派の種の順（Adelie → Chinstrap → Gentoo）になるようにcluster番号を付け直す
cluster_ordered <- factor(km$cluster,
                          levels = order(major),
                          labels = seq_along(major))

# 番号を付け直したclusterとspeciesの対応を確認
tab <- table(Cluster = cluster_ordered,
             Species = peng$species)
tab

tab_df <- as.data.frame(tab)

# clusterとspeciesの対応をタイル図で可視化（Cluster 1が上に来るよう fct_rev で反転）
ggplot(tab_df,
       aes(x = Species,
           y = fct_rev(Cluster),
           fill = Freq)) +
  geom_tile() +
  geom_text(aes(label = Freq), size = 5) +
  scale_fill_gradient(low = "white", high = "steelblue") +
  labs(x = "Species",
       y = "Cluster",
       fill = "Count",
       title = "Cluster vs Species") +
  theme_minimal()

4) 階層的クラスタリングとヒートマップ

# 元の数値データを標準化
X_scaled <- scale(X)
rownames(X_scaled) <- seq_len(nrow(X_scaled)) 

# サンプル間のユークリッド距離を計算
d <- dist(X_scaled)

# Ward法で階層的クラスタリング
hc <- hclust(d, method = "ward.D2")

# デンドログラムを表示
plot(hc,
     labels = FALSE,
     hang = -1,
     main = "Hierarchical Clustering (Ward.D2)")

# デンドログラムを k = 3 のクラスタに分割
grp <- cutree(hc, k = 3)

# PCA平面上にクラスタリング結果を表示
scores$HC_group <- factor(grp)

ggplot(scores, aes(PC1, PC2, color = HC_group)) +
  geom_point(size = 2) +
  coord_equal() +
  theme_minimal() +
  labs(color = "HC group",
       title = "Hierarchical clustering on PCA space")

# speciesとの対応を確認
table(HC_Group = grp, Species = peng$species)

# speciesの情報をアノテーションとして準備
annotation <- data.frame(
  Species = peng$species
)

rownames(annotation) <- rownames(X_scaled)

# ヒートマップ
pheatmap(t(X_scaled),
         cluster_cols = hc,
         cluster_rows = FALSE,
         annotation_col = annotation,
         show_colnames = FALSE,
         main = "Heatmap of Penguin Measurements")

# 形質もクラスタリングさせて、結果を比較
pheatmap(t(X_scaled),
         cluster_cols = hc,
         cluster_rows = TRUE,
         annotation_col = annotation,
         show_colnames = FALSE,
         main = "Heatmap of Penguin Measurements")
⚠️ **GitHub.com Fallback** ⚠️