因子分析からクラスター分析までの流れ

マーケティングの基本的な考え方のひとつであるSTP戦略。
STPとはsegmentation、targeting、positioningの頭文字をとったものです。
消費者をニーズの異なるいくつかのグループに分類し、商品開発の対象を絞り込み、消費者から見た競合他社の商品との相対的な位置付けを定める、という考え方です。
そのS(segmentation)のステップではクラスター分析が良く使われます。


ここではSD法によるアンケートの回答から因子分析を行い、その因子得点をもとにクラスター分析を行うところまでをR言語でやってみたいと思います。


解説が不十分ですが、追々追記したいと思います。


データはここにある「COFFEE.txt」を使わせていただきます。
このデータは「経済・経営のための統計学」という書籍のデータです。


経済・経営のための統計学 (有斐閣アルマ)

経済・経営のための統計学 (有斐閣アルマ)


アンケートの内容はこちらで見られます。

「COFFEE.txt」ダウンロードしたらRのワーキングディレクトリに置いてください。

分析に使うのは「Q7.「缶コーヒー」について以下の考えがあります。それぞれについて、あなたのお考えに近い番号に○をつけてください。」という質問です。
選択肢はそれぞれ「5.非常にそうである」「4.そうである」「3.どちらともいえない」「2.そうではない」「1.まったくそうではない」となっています。


まずはデータの読み込み。

# タブ区切りのテキストファイル。欠損値は“.”で入力されている。
COFFEE <- read.delim("COFFEE.txt", na.strings=".")
Q7 <- na.omit(COFFEE[, c(1, 34:41)])  # 対象者番号とQ7の回答を取り出す。欠損値のある対象者は除く。
colnames(Q7) <- c("対象者番号", "味や品質にこだわる",
                  "時と場合により飲みわける", "気に入ったブランドはない",
                  "人よりよく知っている", "ひとより多く飲む", "特に気にしない",
                  "新製品が出ると試してみる",
                  "好きなブランドを人に勧めることがある")
FACE <- COFFEE[Q7[, 1], c(8, 75, 76)]  # 属性。ただしQ7で欠損値のある人を除く。
                                       # 対象者番号は連番になっていること。
colnames(FACE) <- c("飲用頻度", "年齢", "性別")

これでデータの準備ができました。


さっそく因子分析に入ります。
Rで因子分析を行う場合、factanal関数を使うのが一般的だと思いますが、ここではpsychパッケージのfa関数を使ってみます。
psychパッケージのfa関数に関しては以下を参照してください。


Rで因子分析:psychパッケージ - 人の嫌がることを進んでする
Rによる因子分析法の比較 - 人の嫌がることを進んでする

# 因子分析
library(psych)
fa.parallel(Q7[, -1], fm="ml")  # 平行分析で因子数を決める
                                # 実際にはもっと試行錯誤しながら決めます。
# 最尤法、プロマックス回転で因子分析
Q7.fa <- fa(Q7[, -1], nfactors=2, rotate = "promax", scores=TRUE, fm="ml")
print(Q7.fa, sort=TRUE)
# Factor Analysis using method =  ml
# Call: fa(r = Q7[, -1], nfactors = 2, rotate = "promax", scores = TRUE, 
#     fm = "ml")
#                                      item   ML2   ML1   h2   u2
# 好きなブランドを人に勧めることがある    8  0.78       1.00 0.00
# 新製品が出ると試してみる                7  0.67       0.16 0.84
# 人よりよく知っている                    4  0.63       0.26 0.74
# ひとより多く飲む                        5  0.59       0.55 0.45
# 特に気にしない                          6 -0.51       0.27 0.73
# 気に入ったブランドはない                3 -0.50       0.26 0.74
# 味や品質にこだわる                      1        0.93 0.47 0.53
# 時と場合により飲みわける                2        0.39 0.66 0.34
# 
#                 ML2  ML1
# SS loadings    2.42 1.19
# Proportion Var 0.30 0.15
# Cumulative Var 0.30 0.45
# 
#  With factor correlations of 
#      ML2  ML1
# ML2 1.00 0.49
# ML1 0.49 1.00
# 
# Test of the hypothesis that 2 factors are sufficient.
# 
# The degrees of freedom for the null model are  28  and the objective function was  2.52 with Chi Square of  190.28
# The degrees of freedom for the model are 13  and the objective function was  0.35 
# The number of observations was  80  with Chi Square =  25.64  with prob <  0.019 
# 
# Tucker Lewis Index of factoring reliability =  0.83
# Fit based upon off diagonal values = 0.97
# Measures of factor score adequacy             
#                                                ML2  ML1
# Correlation of scores with factors            0.91 1.00
# Multiple R square of scores with factors      0.84 0.99
# Minimum correlation of factor score estimates 0.67 0.98

X <- Q7.fa$scores  # 因子得点をXに付値
head(X)  # データの確認
#          ML2        ML1
# 1  1.2668908 -0.5189983
# 2 -2.1310404  2.5071716
# 3  2.6297442 -0.2569733
# 4 -0.1980689  0.3260619
# 5 -0.7145165  0.6269295
# 6 -0.8986381 -0.3268858

とりあえずML1を「機能的関与因子」、ML2を「情緒的関与因子」と名づけてみます。


クラスター分析に入る前に、因子得点の分布の様子を確認します。
ヒストグラムはhist関数でも良いのですが、ここではMASSパッケージのtruehist関数を使ってみます。

X <- Q7.fa$scores  # 因子得点をXに付値
head(X)  # データの確認

# 因子得点をヒストグラムで確認
library(MASS)
par(mfrow=c(2,1))
truehist(X[,1], main="情緒的関与因子(ML2)", col="lightgreen")
truehist(X[,2], main="機能的関与因子(ML1)", col="lightblue")
par(mfrow=c(1,1))


次に散布図で確認します。

# 因子得点を散布図で確認
plot(X, pch=3, main="因子得点の散布図(plot)")


ちょっと脱線になりますが、Rではいろいろな散布図の表現が使えるので試してみます

smoothScatter(X, pch=3, main="因子得点の散布図(smoothScatter)")


library(hexbin)
plot(hexbin(X), colramp=function(n){plinrain(n, beg=90, end=15)},
     main="因子得点の散布図(hexbin)")



さて、それではクラスター分析に入りましょう。

# 階層型クラスター分析(ウォード法)
Q7.wa <- hclust(dist(X)^2, "ward")
plot(Q7.wa)
rect.hclust(Q7.wa, k=3, border=1:3) # クラスター数を3として枠を描く
Q7.wa.clust <- cutree(Q7.wa, 3)  # クラスター数を3としてクラスター番号を保存



これでクラスター分析自体は完了です。
実際にはどの手法が適切か、クラスター数を幾つにするかなど、ここでも試行錯誤が発生します。


あとは結果の確認です。

# 構成の確認。
table(Q7.wa.clust)
# Q7.wa.clust
#  1  2  3 
# 27 30 23 
pie(table(Q7.wa.clust), clockwise=TRUE, main="クラスターの構成")



散布図での確認します。

plot(X, pch=21, bg=Q7.wa.clust, main="散布図上で結果の確認")
# クラスターごとに楕円を重ね描き
library(cluster)
E1 <- ellipsoidhull(X[Q7.wa.clust==1, ])
E2 <- ellipsoidhull(X[Q7.wa.clust==2, ])
E3 <- ellipsoidhull(X[Q7.wa.clust==3, ])
lines(predict(E1), col=1)
lines(predict(E2), col=2)
lines(predict(E3), col=3)



最後にクラスターの解釈です。

# クラスターごとの平均値
(ans <- aggregate(X, by=list(Q7.wa.clust), FUN=mean))
#   Group.1        ML2        ML1
# 1       1  1.0042568 -0.0564751
# 2       2 -0.6966410  0.8983492
# 3       3 -0.2702479 -1.1054630

# 箱髭図
par(mfrow=c(1, 2))
boxplot(X[, 1] ~ Q7.wa.clust, main="情緒的関与因子(ML2)", xlab="クラスタ",
        col="lightgreen")
boxplot(X[, 2] ~ Q7.wa.clust, main="機能的関与因子(ML1)", xlab="クラスタ",
        col="lightblue")
par(mfrow=c(1, 1))


# 属性の確認
# 飲用頻度
FREQ <- table(Q7.wa.clust, FACE[, 1])
rownames(FREQ) <- c("clust1", "clust2", "clust3")
colnames(FREQ) <- c("毎日", "週に数回", "月に数回", "年に数回")
mosaicplot(FREQ, main="クラスタ毎の飲用頻度", ylab="飲用頻度", xlab="クラスタ",
           col=TRUE)

# 年代
AGE <- table(Q7.wa.clust, FACE[, 2])
rownames(AGE) <- c("clust1", "clust2", "clust3")
colnames(AGE) <- c("〜20", "21〜25", "26〜30", "31〜35", "36〜40", "41〜60",
                    "61〜")
mosaicplot(AGE, main="クラスタ毎の年齢", ylab="年齢", xlab="クラスタ",
           col=TRUE, las=1)

# 性別
SEX <- table(Q7.wa.clust, FACE[, 3])
rownames(SEX) <- c("clust1", "clust2", "clust3")
colnames(SEX) <- c("男性", "女性")
mosaicplot(SEX, main="クラスタ毎の性別", ylab="性別", xlab="クラスタ",
           col=TRUE)





とりあえずここまで。
ウォード法以外の手法も試してみたいと思います。