第11回（1月8日） Task Check and Weekly Assignment
To Do □ クラスター分析について理解する □ 非階層的クラスター分析を行う □ クラスターの特徴をデータに基づき記述する。
Assignment □ 検査データを因子分析した結果から，因子得点を従属変数としてクラスター分析をする。 □ 各クラスターの特徴を検証。
※提出はソースコードでお願いします！
前回同様，5教科 いつものファイルの読み込みと下処理。済んでいる場合はやらなくてよいです。 ※Workspaceにsampleが入っている場合＝済んでいる場合です。
sample <- read.csv("sample(mac).csv", na.strings = "*")
sample$sex <- factor(sample$sex, labels = c("male", "female"))
因子分析に使うところだけサブセットを作っておきましょうか。
fac_data <- subset(sample, select = c("kokugo", "sansuu", "rika", "eigo", "syakai"))
因子分析は次のようにするのでした。
library(psych)
result.fa <- fa(fac_data, nfactors = 2, fm = "ml", rotate = "varimax", scores = T)
## Loading required package: GPArotation
print(result.fa, sort = T, digit = 3)
## Factor Analysis using method = ml
## Call: fa(r = fac_data, nfactors = 2, rotate = "varimax", scores = T,
## fm = "ml")
## Standardized loadings (pattern matrix) based upon correlation matrix
## item ML2 ML1 h2 u2
## kokugo 1 0.987 -0.144 0.995 0.005
## eigo 4 0.921 -0.105 0.859 0.141
## syakai 5 0.677 0.095 0.467 0.533
## sansuu 2 0.053 0.996 0.995 0.005
## rika 3 -0.085 0.696 0.491 0.509
##
## ML2 ML1
## SS loadings 2.290 1.517
## Proportion Var 0.458 0.303
## Cumulative Var 0.458 0.761
## Proportion Explained 0.602 0.398
## Cumulative Proportion 0.602 1.000
##
## Test of the hypothesis that 2 factors are sufficient.
##
## The degrees of freedom for the null model are 10 and the objective function was 3.309 with Chi Square of 319.4
## The degrees of freedom for the model are 1 and the objective function was 0.077
##
## The root mean square of the residuals (RMSR) is 0.026
## The df corrected root mean square of the residuals is 0.118
## The number of observations was 100 with Chi Square = 7.328 with prob < 0.00679
##
## Tucker Lewis Index of factoring reliability = 0.7925
## RMSEA index = 0.2587 and the 90 % confidence intervals are 0.1061 0.4352
## BIC = 2.723
## Fit based upon off diagonal values = 0.994
## Measures of factor score adequacy
## ML2 ML1
## Correlation of scores with factors 0.998 0.998
## Multiple R square of scores with factors 0.995 0.995
## Minimum correlation of possible factor scores 0.990 0.990
さて，早速本題と行きたいのだけど，今回の因子得点は欠損値があるので，それを平均値で補完しておく。 欠損がないデータだったら，こんなことしなくていいからね。
result.fa$scores <- ifelse(is.na(result.fa$scores), 0, result.fa$scores)
ここかほんとに本題。 非階層的クラスター分析をやります。
非階層的クラスター分析は，ひとまずは，データのサイズが大きくなったときにやる，と理解しておいてください。使われているジャンルは，計量社会学とか，マーケティング分野などです。
手法はk-means法が有名。関数もそのままkmeans。先にクラスター数を決めておかないといけない。理論的仮説があればそれに沿えばいいし，解釈しやすい答えを求めてトライアル＆エラーでもいいかとおもいます。
result.km <- kmeans(result.fa$scores, 4)
result.km$cluster
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
## 3 2 3 3 1 4 4 2 3 3 2 4 3 1 1 1 4 4
## 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
## 2 1 1 2 2 3 3 2 4 3 2 1 2 2 2 4 1 2
## 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54
## 1 3 3 1 1 2 4 3 2 3 4 4 4 4 4 3 4 1
## 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72
## 2 1 3 3 2 4 3 3 3 3 3 3 3 4 1 4 1 4
## 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90
## 2 4 4 4 2 1 2 3 1 4 1 1 1 4 4 4 4 2
## 91 92 93 94 95 96 97 98 99 100
## 3 3 3 3 1 4 3 1 3 4
result.km$size
## [1] 22 20 30 28
あとは，クラスターの中身を探って解釈する，というところは同じ。
sample$cl <- result.km$cluster
# 記述統計で確認したり
describeBy(sample[6:10], sample$cl)
## group: 1
## var n mean sd median trimmed mad min max range skew kurtosis
## kokugo 1 22 72.95 8.52 71.5 72.39 7.41 61 91 30 0.61 -0.74
## sansuu 2 22 77.18 3.00 77.5 77.00 2.22 73 86 13 0.75 1.13
## rika 3 22 53.41 4.27 53.0 53.28 4.45 47 62 15 0.27 -0.93
## syakai 4 22 56.95 14.35 50.0 56.06 14.83 39 86 47 0.49 -1.22
## eigo 5 22 69.23 8.82 68.5 68.67 7.41 52 92 40 0.68 0.33
## se
## kokugo 1.82
## sansuu 0.64
## rika 0.91
## syakai 3.06
## eigo 1.88
## --------------------------------------------------------
## group: 2
## var n mean sd median trimmed mad min max range skew kurtosis
## kokugo 1 20 79.40 6.78 78.5 78.81 5.93 70 94 24 0.62 -0.74
## sansuu 2 20 67.90 2.90 68.5 67.88 3.71 63 73 10 -0.07 -1.27
## rika 3 20 48.00 4.21 46.5 47.69 3.71 43 56 13 0.54 -1.22
## syakai 4 20 55.65 9.63 56.0 55.38 8.15 39 78 39 0.28 -0.31
## eigo 5 20 76.80 8.69 74.5 76.50 6.67 61 94 33 0.33 -0.72
## se
## kokugo 1.52
## sansuu 0.65
## rika 0.94
## syakai 2.15
## eigo 1.94
## --------------------------------------------------------
## group: 3
## var n mean sd median trimmed mad min max range skew
## kokugo 1 30 51.30 7.16 52.5 51.75 6.67 34 62 28 -0.47
## sansuu 2 30 75.17 3.54 74.5 74.83 2.22 69 84 15 0.87
## rika 3 30 53.37 4.61 53.0 53.33 4.45 43 66 23 0.22
## syakai 4 30 42.50 11.37 44.5 43.12 12.60 20 60 40 -0.41
## eigo 5 30 45.40 9.70 45.5 45.71 9.64 25 65 40 -0.25
## kurtosis se
## kokugo -0.61 1.31
## sansuu 0.38 0.65
## rika 0.41 0.84
## syakai -1.10 2.08
## eigo -0.42 1.77
## --------------------------------------------------------
## group: 4
## var n mean sd median trimmed mad min max range skew
## kokugo 1 27 61.19 6.69 61.0 61.17 7.41 46 75 29 -0.07
## sansuu 2 27 65.59 4.19 66.0 65.65 4.45 58 75 17 -0.10
## rika 3 27 46.70 5.57 48.0 46.96 5.93 34 55 21 -0.33
## syakai 4 28 46.54 11.04 45.5 46.17 14.83 30 71 41 0.34
## eigo 5 28 55.96 8.51 55.5 55.58 8.15 41 80 39 0.56
## kurtosis se
## kokugo -0.35 1.29
## sansuu -0.68 0.81
## rika -1.02 1.07
## syakai -1.00 2.09
## eigo 0.41 1.61
boxplot(kokugo ~ cl, data = sample)
# クロス集計表で確認したり
chisq.test(xtabs(~sex + cl, data = sample))
##
## Pearson's Chi-squared test
##
## data: xtabs(~sex + cl, data = sample)
## X-squared = 4.305, df = 3, p-value = 0.2304
# クラスターごとに検定したりして。
source("anovakun_431.txt")
sampleK <- subset(sample, select = c("cl", "kokugo"))
anovakun(sampleK, "As", 4, peta = TRUE)
##
## [ As-Type Design ]
##
## This output was generated via anovakun 4.3.1 at R version 2.15.2.
## It was executed on Mon Jan 7 18:03:15 2013.
##
##
## << DESCRIPTIVE STATISTICS >>
##
## --------------------------------------------
## A N Mean S.D.
## --------------------------------------------
## 2 20 79.4000 6.7777
## 1 22 72.9545 8.5216
## 3 30 51.3000 7.1638
## 4 27 61.1852 6.6912
## --------------------------------------------
##
##
## << ANOVA TABLE >>
##
## == This data is UNBALANCED!! ==
## == Type III SS is applied. ==
## == The number of removed case is 1. ==
##
## ----------------------------------------------------------------------------------------------------
## Source SS df MS F-ratio p-value p.eta^2
## ----------------------------------------------------------------------------------------------------
## A 11536.5987 3 3845.5329 72.3399 0.0000 *** 0.6955
## Error 5050.1286 95 53.1592
## ----------------------------------------------------------------------------------------------------
## Total 16586.7273 98 +p < .10, *p < .05, **p < .01, ***p < .001
##
##
## << POST ANALYSES >>
##
## < MULTIPLE COMPARISON for FACTOR A >
##
## == Shaffer's Modified Sequentially Rejective Bonferroni Procedure ==
## == The factor < A > is analysed as independent means. ==
## == Alpha level is 0.05. ==
##
## --------------------------------------------
## A N Mean S.D.
## --------------------------------------------
## 2 20 79.4000 6.7777
## 1 22 72.9545 8.5216
## 3 30 51.3000 7.1638
## 4 27 61.1852 6.6912
## --------------------------------------------
##
## ---------------------------------------------------------------------------
## Pair Interval t-value df p adj.p
## ---------------------------------------------------------------------------
## 2-3 28.1000 13.3508 95 0.0000 0.0000 2 > 3 *
## 1-3 21.6545 10.5811 95 0.0000 0.0000 1 > 3 *
## 2-4 18.2148 8.4680 95 0.0000 0.0000 2 > 4 *
## 1-4 11.7694 5.6203 95 0.0000 0.0000 1 > 4 *
## 3-4 -9.8852 5.1109 95 0.0000 0.0000 3 < 4 *
## 2-1 6.4455 2.8613 95 0.0052 0.0052 2 > 1 *
## ---------------------------------------------------------------------------
##
##
## output is over --------------------///
階層的クラスタ分析とどう違うか，ということを考えてみてもいいかも。
課題は，http://www1.doshisha.ac.jp/~mjin/data/ にある職業適性能力検査データをつかって，因子分析をし（抽出方法の指定，回転の指定もする，因子得点も算出）クラスター分析してください。