大量のサンプルからデータを取得したときに、そのデータの回答を基に、回答傾向が似通っているサンプルをグルーピング（クラスター（群）わけ）したいということはよくある。今回はそのための分析方法ークラスター分析について説明する。

なお，クラスター分析には大きく分けて2つの方法がある。

  • 階層的クラスタリング法
  • k-means法

以下でこれらを説明する。

1 階層的クラスタリング法

1.1 基本的な考え方

階層的クラスタリング法は、回答傾向が似ているかどうかを回答間の距離で計算し、距離が近いもの同士を結合させていくという結合させていく、という方法である。ここでいう距離というのは、例えば2次元で考えるならば、以下の図で示される長さのことである。このような長さを全てのサンプルの組み合わせに対して計算させて、距離が近い者同士を順に結び付けていく。

上の図だと2次元（つまり変数が2つ）の場合の距離を示しているが、1次元（変数が1つ）でもよいし，3次元（変数が3つ）でもよい．，実際には多次元，つまり任意の個数の変数を分析対象に含めた場合でも同様に「距離」を計算することができる。それぞれの場合の距離の計算方法は以下の通りである。

\[ \begin{align} 一次元の場合の距離：& D = \sqrt{(x_1 - x_2)^2} \\ 二次元の場合の距離: & D = \sqrt{(x_1 - x_2)^2 + (y_1 - y_2)^2} \\ 三次元の場合の距離: & D = \sqrt{(x_1 - x_2)^2 + (y_1 - y_2)^2 + (z_1 - z_2)^2} \\ \vdots\\ n次元の場合の距離: & D = \sqrt{\sum_{i=1}^{n} (x_{1i} - x_{2i})^2} \end{align} \]

もちろん，実際にこのような数式を自分で打って計算させるわけではなく，Rには距離を計算させる関数が用意されているので，それを使う．

1.2 基本的なクラスター分析の実施手順

Rでクラスター分析を実施する際には，

(1)まずデータに含まれる各サンプル同士の距離データを計算しする

(2)その距離データを基に樹形図（デンドログラム）を描く

(3)デンドログラムを見ながらクラスター数を決定する

(4)決定したクラスター数でクラスター分けを行う

という流れとなる．

具体的な手順は以下の通りである。 dist()関数とhclust()関数を使う。 dist()関数はデータ間の距離を計算する関数であり，hclust()関数が階層的クラスタリングを行う関数である。 基本的な引数としては，それぞれ以下の引数を取る．

hclust()関数

  • 第1引数: dist()関数で計算した距離行列.
  • 第2引数 method（省略可）: クラスターの結合方法を指定する。デフォルトはward法である。

では実際に以下の例題を使ってクラスター分析を行ってみよう．

このデータは、2022年11月に著者のゼミ生が行った大学生向けに行った余暇の過ごし方についてのアンケート結果である。全体で222人が回答を寄せてくれた。このアンケートでは性格特性として「外向性」「親和性（協調性ともいう）」「勤勉性」「神経症傾向」「開放性（新しいものの受入れに積極的か消極的か）」のデータと取得している。この5つの性格特性を基に、222人のデータをグループ分けせよ。

1.2.1 データの読み込みと内容確認

まずはデータの読み込んで，内容を確認してみよう．

# データの読み込み
data <- read.csv("practice/example_14_1.csv")

# データの確認
head(data)
##   性別 学年   現在の住居形態
## 1   男  1年 一人暮らし・下宿
## 2   女  1年 一人暮らし・下宿
## 3   女  1年             実家
## 4   男  1年             実家
## 5   女  1年 一人暮らし・下宿
## 6   男  1年             実家
##   自身の金銭感覚..どちらか近いと思う方に答えてください. 外向性 親和性 勤勉性
## 1                                  浪費するタイプである    3.0    5.0    3.5
## 2                                  貯蓄するタイプである    2.5    4.5    3.5
## 3                                  浪費するタイプである    5.5    5.0    2.5
## 4                                  浪費するタイプである    6.0    7.0    3.0
## 5                                  貯蓄するタイプである    5.5    4.5    2.0
## 6                                  貯蓄するタイプである    3.5    4.5    3.5
##   神経症傾向 開放性 必要以外の睡眠をする.昼寝など.
## 1        4.0    2.5                              6
## 2        7.0    4.5                              5
## 3        6.5    4.5                              5
## 4        6.5    5.5                              6
## 5        3.5    4.5                              7
## 6        4.5    4.5                              7
##   家でYouTubeやNetflixなど動画を見て過ごす
## 1                                        7
## 2                                        7
## 3                                        7
## 4                                        5
## 5                                        7
## 6                                        6
##   パチンコや競艇などギャンブルに出かける ショッピングに行く
## 1                                      1                  5
## 2                                      1                  6
## 3                                      1                  6
## 4                                      1                  6
## 5                                      1                  6
## 6                                      1                  4
##   部活.サークルや地域クラブなどで運動する ボランティア活動に参加する
## 1                                       2                          1
## 2                                       1                          1
## 3                                       1                          1
## 4                                       7                          6
## 5                                       1                          2
## 6                                       7                          7
##   電話をする..相手からかかってきて話すものも含む. 旅行に行く
## 1                                               6          4
## 2                                               7          7
## 3                                               7          6
## 4                                               5          7
## 5                                               5          5
## 6                                               5          3
##   勉強や自己研鑽.自分磨き.をする ウォーキングなどの軽い運動を行う
## 1                              2                                3
## 2                              5                                2
## 3                              7                                1
## 4                              7                                7
## 5                              3                                2
## 6                              4                                3
##   読書.漫画を含む. 友人と会う 映画を見に行く 自身の余暇の内容には満足しているか
## 1                6          6              3                                  4
## 2                6          5              4                                  4
## 3                5          6              4                                  6
## 4                5          7              5                                  7
## 5                5          6              3                                  6
## 6                2          7              4                                  3
##   余暇を十分にとれていると感じているか 自身の余暇を改善したいと思いますか
## 1                                    4                                  7
## 2                                    7                                  2
## 3                                    7                                  7
## 4                                    7                                  1
## 5                                    7                                  5
## 6                                    6                                  4

質問項目文がそのまま変数名として入力されてしまっているデータであるため（本来はQ1やQ2という項目番号で整理すべき！！），非常に長ったらしい出力になってしまっているが，この結果から，5つの性格特性（外向性・親和性・勤勉性・神経症傾向・開放性）がそれぞれそのまま名前となった変数が含まれていることが確認できた．

性格特性の部分だけを取り出してみよう．

head(data[c("外向性", "親和性", "勤勉性", "神経症傾向", "開放性")])
##   外向性 親和性 勤勉性 神経症傾向 開放性
## 1    3.0    5.0    3.5        4.0    2.5
## 2    2.5    4.5    3.5        7.0    4.5
## 3    5.5    5.0    2.5        6.5    4.5
## 4    6.0    7.0    3.0        6.5    5.5
## 5    5.5    4.5    2.0        3.5    4.5
## 6    3.5    4.5    3.5        4.5    4.5

1.2.2 標準化と距離の計算

次に，距離を計算する．そのための関数としてdist()関数がある．引数は以下の通りである．

  • 引数: 分析対象の変数で構成されたデータフレーム
# データの標準化
data_std <- scale(data[c("外向性", "親和性", "勤勉性", "神経症傾向", "開放性")])

# 距離の計算
dist_data <- dist(data_std)

上記では，データを標準化してから距離を計算している．標準化を行うことで，各変数の尺度が異なる場合でも，距離の計算が適切に行われる．

ちなみに距離の計算は，すべての点に対して，自分以外のすべての点との距離が計算されており，それは以下のような行列で表現されている（縦横ともに5つまで表記しているが，実際には，222点のサンプルがあることから，222$$222の行列である．

##       1     2     3     4     5
## 1 0.000 2.790 2.990 4.075 2.521
## 2 2.790 0.000 2.171 3.556 3.495
## 3 2.990 2.171 0.000 2.254 2.391
## 4 4.075 3.556 2.254 0.000 3.638
## 5 2.521 3.495 2.391 3.638 0.000

1.2.3 階層的クラスタリング

次に，hclust()関数を使ってデンドログラムを描くためのデータを取得する．引数は以下の通り．

  • 第1引数: dist()関数で計算した距離行列
  • 第2引数 method（省略可）: クラスターの結合方法を指定する。省略した場合にはデフォルトはward.d2法となる。通常は省略でよい．

結果はhcというオブジェクトに収めた上で，plot()関数を使って樹形図（デンドログラムと呼ぶ）の形で可視化する．plot()関数の基本引数は以下の通り．

  • 第1引数: hclust()関数で得られたオブジェクト
  • 第2引数 hang（省略可）: クラスターの枝の先端をどれだけ下げるか．デフォルトは0.1．負の値にした場合，枝の先端が根元と同じ高さになる．
# 階層的クラスタリング
hc <- hclust(dist_data)

# デンドログラムの描画（hangをデフォルトにした場合）
plot(hc)

# デンドログラムの描画 その2（hangを-1にした場合）
plot(hc, hang=-1)

デンドログラムは，データ間の距離に基づいてクラスタリングされた結果を木構造で表現したものである．デンドログラムを見ることで，データがどのようにクラスタリングされたかを視覚的に把握することができる．すなわち，Height=0の箇所から，距離が近いもの同士が結びついていき，最終的に全てのデータが1つのクラスターにまとまるまで結合が進む様子がわかる．

1.2.4 クラスターの数の決定

デンドログラムはあくまで結びつきの過程を示しているだけであり，デンドログラムを見ながらどのHeigtのラインで樹形図をカットするかを分析者が決定する必要がある．

カットの基準としては以下の2つの基準がある．

  • デンドログラムの上下の分離ラインがハッキリしている箇所で切る
  • でき上るクラスター数がその後の分析での利用の観点から適切な数になる箇所で切る

例えば，今回のデンドログラムであれば，比較的上下での分離ラインがハッキリしていて，出来上がるクラスター数として手ごろな数になる箇所として，以下の図のような赤のラインで切るというのは，一つの見方として十分ありえるだろう．赤いラインで切ると大きく5つのクラスタができる（クラスタの境目を青いラインで描いている）．

このように，最終的には分析者がデータの特性や目的に応じて適切なクラスター数を決定する必要がある．

1.2.5 クラスター分け

クラスター数を決定したら，cutree()関数を使って実際にデータのクラスター分けを行う．cutree()関数の引数は以下の通りである．

  • 第1引数: hclust()関数で得られたオブジェクト
  • 第2引数 k: クラスター数

結果として，各サンプルがどのクラスターに配属されるかを示すベクトルが返される．

# クラスター分け
cluster_res <- cutree(hc, k=5)
print(cluster_res)
##   [1] 1 2 2 2 3 1 4 2 1 1 1 3 1 2 1 3 3 1 2 3 4 2 2 3 4 5 5 5 2 2 1 5 1 1 2 5 1
##  [38] 5 2 5 3 1 1 1 1 1 1 2 3 1 3 1 2 5 2 2 1 2 5 1 5 1 5 2 1 3 1 1 1 5 1 2 2 2
##  [75] 2 2 1 1 1 1 4 4 1 4 3 1 5 3 3 1 1 5 5 3 1 2 2 2 5 1 2 1 3 5 1 2 1 2 2 5 5
## [112] 1 1 5 1 3 5 3 2 2 5 1 4 1 1 5 2 5 2 1 5 5 1 1 2 5 4 5 2 2 2 5 1 5 2 3 1 2
## [149] 3 2 2 1 3 1 5 1 3 5 2 1 5 1 4 2 5 5 3 1 2 2 1 3 5 1 2 3 1 1 3 1 1 1 3 2 2
## [186] 1 4 1 3 3 2 2 1 3 4 1 3 2 2 1 5 3 1 3 1 1 1 1 3 5 5 1 1 5 5 1 5 2 5 1 4 1

以下では，得られたクラスター分けのデータを元のデータに追加して，クラスターごとのサンプル数をtable()関数を使って確認している．

# クラスターの割り当て
data$cluster_res <- cluster_res

# クラスターごとのサンプル数
table(data$cluster_res)
## 
##  1  2  3  4  5 
## 80 54 33 12 43

1.2.6 クラスターの解釈

最後に，各クラスターの特徴を把握するために，各クラスターの平均値を計算してみる．なお，比較しやすくするために，各性格特性値は全サンプルの平均で中心化（要するに各サンプルの値-全体平均）している

# 中心化
data$m外向性 <- data$外向性 - mean(data$外向性)
data$m親和性 <- data$親和性 - mean(data$親和性)
data$m勤勉性 <- data$勤勉性 - mean(data$勤勉性)
data$m神経症傾向 <- data$神経症傾向 - mean(data$神経症傾向)
data$m開放性 <- data$開放性 - mean(data$開放性)

#クラスタ番号を要因型に
data$cluster_res <- as.factor(data$cluster_res)

# ロング型に変換
df <- pivot_longer(data, 
        cols = c("m外向性", "m親和性", "m勤勉性", "m神経症傾向", "m開放性"), 
        names_to = "性格特性", 
        values_to = "値")

# ggplotで可視化　グループごとに各性格特性の平均値を棒グラフで表示
ggplot(df,aes(x=性格特性, y=値,fill=性格特性)) +
  geom_bar( stat="summary", fun.data = mean_se) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) + # x軸のラベルを縦向きに
  facet_wrap(~cluster_res, nrow=1)

この結果から，各クラスターの特徴が把握できる．例えば，クラスター1は開放性や外向性，勤勉性も低いことから内向的であまり活動的でない性格であることがうかがえる．

逆にクラスター3などは極めて外向的で新しいものに対して積極的で,神経質なところもあまりないが，そこまで親和的というわけでもない，ということがわかる．

このようにクラスター分けを行ったら，最後にそのクラスターがどういうクラスターなのかを，各クラスターの平均値などの統計量を用いて解釈することが重要である．

エルボー法

1つ目の基準でのカットをサポートする方法として，エルボー法という方法がある．これは，クラスター数を変化させたときのクラスタ内誤差平方和（Within-Cluster Sum of Squares, WCSS）をプロットし，そのグラフの曲線が肘（エルボー）のようになっているところ，すなわちその点を境に急激に値が高くなっている点を適切なクラスター数とする方法である．

エルボーグラフを描くにはfactoextraパッケージに含まれるfviz_nbclust()関数を使う．この関数は以下の引数を取る．

  • 第1引数: 分析対象のデータ(標準化済みのものがよい)
  • 第2引数: クラスタリング手法（ここでは階層的クラスタリングを行ったので，hcutを指定）
  • 第3引数 method: クラスター数を決定するための基準を指定する．ここでは”wss”を指定する．
# エルボーグラフのプロット
library(factoextra)
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
fviz_nbclust(data_std, hcut, method = "wss")

この結果からは，あまりはっきりとはクラスター数を読み取れない．

このように，エルボー法の場合，肘が明確に見える場合もあれば，見えにくい場合もあるので，見えにくい場合には，他の基準も併用して判断することが望ましい．

シルエット係数

クラスタ数を決定する際の指標としてシルエット係数とそれを図示したシルエット図というものが利用されることがある．

シルエット係数（Silhouette Coefficient）は、以下の計算に基づいて各データポイントのクラスタ適合度を評価する指標である．

\[ \begin{align} s(i) &= \frac{b(i) - a(i)}{max(a(i), b(i))}\\ \text{ただし} \\ a(i) &: 各データポイントi について、同じクラスタ内の他の点との平均距離\\ b(i) &: 各データポイントi に対して、最も近い異なるクラスタの点との平均距離 \end{align} \]

要するに，自クラスターに所属する他のデータ点らとの距離の平均と，最も近くにある他のクラスターに所属するデータ点らとの距離の平均を比較しているものである．

この\(s(i)\)は-1から1の値をとり，1に近いほどそのデータポイントが適切なクラスタに所属していることを示す．逆に0に近いとクラスタの境界上に存在していることを示し，負の値の場合，不適切なクラスタに分類されている可能性があることを示す．

Rではclusterライブラリを用いて，シルエット係数の計算とシルエット図の描画することができる．用いるのはsilhouette()関数とplot()関数で，それぞれ引数は以下の通りである．

  • silhouette()関数
    • 第1引数: クラスタリング結果
    • 第2引数: 距離を求めるdist()関数の出力結果
  • plot()関数
    • 第1引数: silhouette()関数の結果
    • main: グラフのタイトル
    • col: クラスターの色の設定．通常は「1:クラスター数」を設定
    • border: クラスターの境界線の色の設定．通常はNA．
library(cluster)
# クラスター分けまですでに済んでいる前提

# シルエット係数の計算
sil <- silhouette(cluster_res, dist_data)

# シルエット図の描画
plot(sil, main = "シルエット図", col = 1:5, border = NA)

シルエット図からいずれのクラスターにもシルエット係数が負となってしまっているデータがあることがわかり，5つのクラスターに分けるのは適切とは言えないことが示唆されている．

さらに，全体のクラスタリングの質を評価するには、シルエット係数の平均値を確認する．上の図の下の方に記載されているが，以下のようにsummray()関数の出力から得ることもできる．

summary_res<-summary(sil)
summary_res$avg.width
## [1] 0.1187263

この値は以下の基準で評価するとされている[参考]。

  • 値が 0.71 以上: クラスタリングは良好。
  • 値が 0.51 ～ 0.7: クラスタリングは程々。
  • 値が 0.26 ～ 0.5: クラスタリングは許容範囲だが，他も試したほうが良い．
  • 値が 0.25 未満: クラスタリングは不適切。

この結果から，やはり5つのクラスターに分けた今回のクラスタリングは不適切であることがわかる．

また，このシルエット係数の平均値を用いて，最適なクラスター数を探ることができる． 以下はクラスター数を2から10まで変化させたときのシルエット係数の平均値を出力させたものである．

for(i in 2:10){
  cluster_res <- cutree(hc, k=i)
  sil <- silhouette(cluster_res, dist_data)
  sil_avg <- summary(sil)$avg.width
  cat("クラスタ数:", i, " 平均シルエット幅:", sil_avg, "\n")
}
## クラスタ数: 2  平均シルエット幅: 0.1185557 
## クラスタ数: 3  平均シルエット幅: 0.09729979 
## クラスタ数: 4  平均シルエット幅: 0.1071872 
## クラスタ数: 5  平均シルエット幅: 0.1187263 
## クラスタ数: 6  平均シルエット幅: 0.1145162 
## クラスタ数: 7  平均シルエット幅: 0.1155542 
## クラスタ数: 8  平均シルエット幅: 0.111742 
## クラスタ数: 9  平均シルエット幅: 0.1131134 
## クラスタ数: 10  平均シルエット幅: 0.08367037

この結果から，クラスター数は5がまだ比較的マシである，という結果であった．

1.3 階層的クラスタの種類

上記は階層的クラスタリングの基本的な手法であるが，階層的クラスタリングにはインプットする距離データをどういうものにするのか，デンドログラムを作成する際のグルーピングアルゴリズムに何を使うのか，いくつかの種類がある．

1.3.1 距離データの算出方法

冒頭に示した高校数学で学習する方法で計算される距離を幾何学的距離，もしくはユークリッド距離とよぶ．dist()関数を使用して距離を計算する際には，デフォルトではユークリッド距離が用いられる．

ただ，dist()関数には第2引数でmethod引数があり，ここに以下の表のモノを指定することで，ユークリッド距離以外の距離を計算することができる．

method 距離尺度 計算方法 特徴・適用場面
"euclidean"（デフォルト） ユークリッド距離 2つのサンプルの各変数の差の2乗和の平方根 \[ d(x, y) = \sqrt{\sum_{i} (x_i - y_i)^2} \] 直線距離を測る標準的な距離尺度。連続変数に適用され、主にクラスタリングや主成分分析（PCA）に使用される。
"manhattan" マンハッタン距離 2つのサンプルの各変数の差の絶対値の和 \[ d(x, y) = \sum_{i} |x_i - y_i| \] 「L1距離」とも呼ばれ、碁盤の目のような距離を測る。外れ値に強く、金融データや経済データの解析に適用される。
"maximum" チェビシェフ距離 2つのサンプルの各変数の差の絶対値の最大値 \[ d(x, y) = \max_{i} |x_i - y_i| \] 最大の要素間差を基準にするため、「最悪ケース」に対応。品質管理や最大誤差を評価する場面で使用される。
"canberra" キャンベラ距離 各変数の差の絶対値を、それぞれの変数の絶対値の和で割ったものの総和 \[ d(x, y) = \sum_{i} \frac{|x_i - y_i|}{|x_i| + |y_i|} \] 小さな値に対して敏感な距離尺度。環境データや医療データなど、値の範囲が広いデータに適用される。
"binary" バイナリ距離 2つのサンプルの各変数が同じ場合は1、異なる場合は0として、その合計を変数数で割ったものを，さらに1から引く． \[ \text{similarity} = \frac{\sum (x_i = y_i)}{\text{total elements}} \] \[ d(x, y) = 1 - \text{similarity} \] 0/1 のバイナリデータ（例：アンケートの Yes/No）に適用。テキストマイニングやユーザー行動分析に有用。
"minkowski" （利用する場合には，さらに第3引数としてp=で任意の指数を与える．） ミンコフスキー距離 ユークリッド距離とマンハッタン距離の一般化形。パラメータ \(p\) を用いる \[ d(x, y) = \left( \sum_{i} |x_i - y_i|^p \right)^{1/p} \] \(p\) の値によって距離の性質を調整可能（例：\(p=1\) はマンハッタン、\(p=2\) はユークリッド）。機械学習やパターン認識に適用。

1.3.2 グルーピングアルゴリズム

一方，グルーピングアルゴリズムとしては，以下のものがある．省略した場合は"complete"が使われる．

method クラスタリング手法 特徴・適用場面
"ward.D2" Ward法 クラスター内の分散の増加量を最小化するようにクラスターを結合する方法。クラスター内の分散が小さく、クラスター間の分散が大きいクラスタリングを行う。
"ward.D" Ward法 "ward.D2"と同じだが，距離データにユークリッド距離の2乗を与える時にはこちらを使う．
"single" シングルリンク法 クラスター内の最も近い2点間の距離をクラスター間の距離とする方法。クラスター内の最も近い2点を基準にクラスタリングを行う。
"complete"
（デフォルト値）		 コンプリートリンク法 クラスター内の最も遠い2点間の距離をクラスター間の距離とする方法。クラスター内の最も遠い2点を基準にクラスタリングを行う。
"average" ウォーキングリンク法 クラスター内の全ての点間の距離の平均をクラスター間の距離とする方法。クラスター内の全ての点間の距離の平均を基準にクラスタリングを行う。
"mcquitty" エンブリア法 クラスター内の全ての点間の距離の最大値をクラスター間の距離とする方法。クラスター内の全ての点間の距離の最大値を基準にクラスタリングを行う。
"median" メディアンリンク法 クラスター内の全ての点間の距離の中央値をクラスター間の距離とする方法。クラスター内の全ての点間の距離の中央値を基準にクラスタリングを行う。
"centroid" セントロイド法 クラスター内の全ての点間の距離の平均をクラスター間の距離とする方法。クラスター内の全ての点間の距離の平均を基準にクラスタリングを行う。

これらの違いによって，クラスタリングの結果が変わることがあるので，思ったような結果とならなかった場合には，これらの設定を変えてみるのも良いだろう．ただし，特に距離データの算出方法は，データの特性に合わせて適切なものを選ぶ必要があるので，注意が必要である．

2 k-means法

2.1 基本的な考え方

k-means法は，以下の手順でクラスタリングを行う．

  1. クラスタ数 k を決める。
  2. k個の初期クラスタ中心をランダムに選ぶ。
  3. 各データ点を最も近いクラスタに割り当てる。
  4. 新しいクラスタ中心を計算（各クラスタ内の平均値）。
  5. クラスタ中心の位置が収束するまで 3 ～ 4 を繰り返す。

k-means法は階層化クラスタリング法に比べて高速なアルゴリズムでり，扱う変数やサンプル数が多い場合にも適用しやすい．

ただ，この手法は，クラスタ中心をランダムに選ぶため，初期値によっては異なるクラスタリング結果が得られることがある．そのため，複数回の計算を行い，最も良い結果を選ぶという方法が取られることが多い．また，最初にクラスタ数を決めておく必要がある．

2.2 k-means法の実施手順

先ほどの例題を今度はk-means法でクラスタリングしてみよう．

2.2.1 データの読み込み

すでに読み込んでいるデータであるが，一応改めて読み直しておこう． また，先ほどはhead()を使ってデータに含まれる変数を確認したが，今度は変数名だけを出力させてみよう．

data <- read.csv("practice/example_14_1.csv")
colnames(data)
##  [1] "性別"                                                 
##  [2] "学年"                                                 
##  [3] "現在の住居形態"                                       
##  [4] "自身の金銭感覚..どちらか近いと思う方に答えてください."
##  [5] "外向性"                                               
##  [6] "親和性"                                               
##  [7] "勤勉性"                                               
##  [8] "神経症傾向"                                           
##  [9] "開放性"                                               
## [10] "必要以外の睡眠をする.昼寝など."                       
## [11] "家でYouTubeやNetflixなど動画を見て過ごす"             
## [12] "パチンコや競艇などギャンブルに出かける"               
## [13] "ショッピングに行く"                                   
## [14] "部活.サークルや地域クラブなどで運動する"              
## [15] "ボランティア活動に参加する"                           
## [16] "電話をする..相手からかかってきて話すものも含む."      
## [17] "旅行に行く"                                           
## [18] "勉強や自己研鑽.自分磨き.をする"                       
## [19] "ウォーキングなどの軽い運動を行う"                     
## [20] "読書.漫画を含む."                                     
## [21] "友人と会う"                                           
## [22] "映画を見に行く"                                       
## [23] "自身の余暇の内容には満足しているか"                   
## [24] "余暇を十分にとれていると感じているか"                 
## [25] "自身の余暇を改善したいと思いますか"

2.2.2 k-means法の実行

k-means法を実行する関数はkmeans()関数である．引数は以下の通りである．

  • 第1引数: 分析対象の変数で構成されたデータフレーム(標準化済みのものがよい)
  • 第2引数 centers: クラスター数
  • 第3引数 nstart（省略可）: 初期値を変えて何回計算するか．省略すると1に設定される．

特に第3引数については，k-means法では初期値をどこに採るかによって結果が変わることや，通常では初期値はランダムに選ばれることになるため，計算する度に異なる結果が返される．そのためnstartに10や25などの値を設定して，複数回の計算を行い，その中で最も良い結果（クラスタ内誤差平方和が最小のもの）を選ぶという方法が取られることが多い．

また，k-means法を実際にRで行う際には，初期値の選択の再現性を持たせるために，set.seed()関数を使って乱数のシードを固定することが推奨される． これにより，k-means()を実行する度に結果が変わってしまうことを防ぐことができる．set.seed()関数の引数には任意の数値を指定する．この値が変わらない限り，同じ乱数が生成されるため，k-means()の結果も同じになる．

誤差平方和とは

各データ点とそのデータ点が属するクラスタの中心との距離の2乗を，すべてのデータ点について計算し、それらを足し合わせたもの．数式で表すと以下の通りとなる． \[ \sum_{i=1}^{n} \sum_{j=1}^{k} r_{ij} ||\boldsymbol{x}_i - \boldsymbol{\mu}_j||^2 \] ただし，\(r_{ij}\) はデータ点 \(\boldsymbol{x}_i\) がクラスタ \(j\) に属する場合は 1，そうでない場合は 0 となる指示関数であり，\(\boldsymbol{\mu}_j\) はクラスタ \(j\) の中心点である．また \(||\cdots||\) はユークリッド距離（ベクトルの大きさ）を表す．

この値が小さいということは、それだけ各クラスタが凝集していることを意味し、クラスタのまとまりの良さの指標となる． 
# 乱数のシードを固定
set.seed(13)

# データの標準化
data_std <- scale(data[c("外向性", "親和性", "勤勉性", "神経症傾向", "開放性")])

# k-means法の実行
kmeans_res <- kmeans(data_std, centers=5, nstart=25)

# クラスターの割り当て
data$cluster_res_kmeans <- kmeans_res$cluster

# クラスターごとのサンプル数
table(data$cluster_res_kmeans)
## 
##  1  2  3  4  5 
## 58 25 49 34 56

上の例では，k-means()の結果をkmeans_resとして受け取っている．この中のcluster要素に各サンプルがどのクラスタに属するかが格納されている．

あとは先と同様に，これを元のデータに追加して，クラスターごとのサンプル数を確認している．

先ほどの結果とクラスに含まれるサンプル数が変わっていることが確認できる．なお，クラスター番号自体は，k-means法の実行時にランダムに割り当てられるため，先ほどのクラスター番号とは対応はしない．

2.2.3 クラスターの解釈

先ほどと同様に，各クラスターの特徴を把握するために，各クラスターの平均値を計算してみる．

# 中心化
data$m外向性 <- data$外向性 - mean(data$外向性)
data$m親和性 <- data$親和性 - mean(data$親和性)
data$m勤勉性 <- data$勤勉性 - mean(data$勤勉性)
data$m神経症傾向 <- data$神経症傾向 - mean(data$神経症傾向)
data$m開放性 <- data$開放性 - mean(data$開放性)

#クラスタ番号を要因型に
data$cluster_res_kmeans <- as.factor(data$cluster_res_kmeans)

# ロング型に変換
df <- pivot_longer(data, 
        cols = c("m外向性", "m親和性", "m勤勉性", "m神経症傾向", "m開放性"), 
        names_to = "性格特性", 
        values_to = "値")

# ggplotで可視化　グループごとに各性格特性の平均値を棒グラフで表示
ggplot(df,aes(x=性格特性, y=値,fill=性格特性)) +
  geom_bar( stat="summary", fun.data = mean_se) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) + # x軸のラベルを縦向きに
  facet_wrap(~cluster_res_kmeans, nrow=1)

先ほどの結果と比較して，クラスタ1は先ほどのクラスタ1が，クラスタ5は先ほどのクラスタ5が，それぞれ比較的類似した特性を持ったクラスタとして形成されているが，クラスタ2と4では開放性・外向性はクラスタ3と似てはいるものの，その他では大きく異なっている．また，クラスタ3と先ほどのクラスタ2と神経症傾向の高さでは似ているが，外向性や親和性では逆向きといえる性格特性となっている．

このようにk-means法と階層的クラスタリング法では，クラスタリングの結果が異なる．

いずれが良いかは特に定めがあるわけではないので，出てきた結果の解釈のしやすさや，その後の分析での利用可能性を考慮して選択する．

クラスタ数の決定

先ほどの階層的クラスタリングと同様に，k-means法でもエルボー法やシルエット係数でのクラスタ数の決定ができる．

まずはエルボー法．

library(factoextra)
# エルボー法
fviz_nbclust(data_std, kmeans, method = "wss")

さっきとそれほど変わらない結果になっているが，4-5あたりで一旦なだらかになっている．

次いで，シルエット係数を算出して図示してみる．

library(cluster)
# クラスター分けまですでに済んでいる前提
# シルエット係数の計算
dist_data <- dist(data_std)
sil <- silhouette( kmeans_res$cluster, dist_data)
# シルエット図の描画
plot(sil, main = "シルエット図", col = 1:5, border = NA)

fviz_nbclust(data_std, kmeans, method="silhouette")

この結果から，今回のクラスタリングは先ほどの階層的クラスタリングよりも負の値となっているものが少なくなっているのがわかる．

for(i in 2:10){
  set.seed(13)
  kmeans_res <- kmeans(data_std, centers=i, nstart=25)

  sil <- silhouette(kmeans_res$cluster, dist_data)
  sil_avg <- summary(sil)$avg.width
  cat("クラスタ数:", i, " 平均シルエット幅:", sil_avg, "\n")
}
## クラスタ数: 2  平均シルエット幅: 0.1982031 
## クラスタ数: 3  平均シルエット幅: 0.1853639 
## クラスタ数: 4  平均シルエット幅: 0.1758217 
## クラスタ数: 5  平均シルエット幅: 0.1852062 
## クラスタ数: 6  平均シルエット幅: 0.1633759 
## クラスタ数: 7  平均シルエット幅: 0.1788516 
## クラスタ数: 8  平均シルエット幅: 0.1736349 
## クラスタ数: 9  平均シルエット幅: 0.1731275 
## クラスタ数: 10  平均シルエット幅: 0.174598

この結果から，今回のデータにおいては，k-means法ではクラスタ数が2の時が最も良いとなる． ただ，こちらで示した基準によると，2でも基準を満たしていないため，今回の結果は基準に照らせばあまり良いとは言えない．

かといって，シルエット平均幅で基準を満たさなかったからといって，クラスタ数を変えるというのもあまり良い方法ではない．シルエット幅は使う変数が大きくなると値が小さくなる傾向があったりもするし，シルエット幅が高い結果でも、ビジネスや研究において「意味のある分割」とは限らない．最終的には出てきた結果の利用可能性や解釈可能性で判断することが重要である．

次のデータは主成分分析 ・因子分析の時の課題で出したデータである．

  • Q1: 何かを選ぶ時は些細なことは気にしない
  • Q2: 重要なことを決めるとき，些細な問題にこだわり多くの時間を費やす
  • Q3: 決断を先送りにしようとする
  • Q4: 重要な問題は，直前まで決定を延ばす
  • Q5: 他の人が何を選ぶのか気になる
  • Q6: 自分が選んだものが他の人と違うと不安になる
  • Q7: 決断した後，その決断に後悔しない
  • Q8: 決めた後で，間違ってしまったと思うことがよくある

このデータは因子分析を実施したところ，因子数4，Varimax回転，Pa推定の場合に，次のような因子構造となる．

  • 因子1（熟慮）: Q1, Q2
  • 因子2（先延ばし）: Q3, Q4
  • 因子3（他者参照）: Q5, Q6
  • 因子4（不安）: Q7, Q8

実際にこの結果は文献(齊藤・緑川, 2016)の通りであり，もともとこれらの質問項目は文献に記載されている項目から因子負荷量の大きいものを2つずつ選んだものである．

さて，このデータについて，以下の問に答えよ．

  1. 上記の因子について，それぞれ紐づいた質問の回答値を平均したものとをその因子の尺度として，熟慮，先延ばし，他者参照，不安の4つの尺度データを作成せよ．なお，逆転項目が含まれていることに注意すること．
  2. 作成した4つの因子の尺度データを用いて，階層的クラスタリングを実施し，デンドログラムを表示させよ．
  3. デンドログラムを基に，各自なりの判断でクラスター分けを実施し，各クラスタに含まれるサンプル数を示せ．
  4. クラスター分けの結果を元に，各クラスタの平均値から各クラスターの特徴を述べよ．
  5. 以上のことをk-means法でも実施し，結果を階層的クラスタリングの場合と比較せよ．