第10回（12月18日） 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)
ここかほんとに本題。
階層的クラスター分析をやります。階層的クラスター分析は，距離をもとに分析するので，距離行列をつくります。
kyori <- dist(result.fa$scores, method = "euclidean")
この距離行列をもとに階層的クラスター分析。特にパッケージは必要としません。
result.hcl <- hclust(kyori, method = "ward")
plot(result.hcl)
自分でいくつのクラスターにするのかを決めて，カットする。
ans <- cutree(result.hcl, k = 4)
カットした結果を元のデータセットに追加しよう
sample$cl <- ans
head(sample)
## ID class sex height weight kokugo sansuu rika syakai eigo cl
## 1 1 A male 170.2 64.16 34 74 43 20 28 1
## 2 2 B male 165.7 75.63 82 63 44 54 72 2
## 3 3 C male 157.8 62.64 50 74 55 26 44 1
## 4 4 A male 161.6 69.57 57 75 55 46 44 1
## 5 5 B male 161.1 60.23 74 73 54 41 65 3
## 6 6 C male 156.2 54.99 NA 58 38 47 58 4
あとは，クラスターの中身を探って解釈。
# 記述統計で確認したり
describeBy(sample[6:10], sample$cl)
## group: 1
## var n mean sd median trimmed mad min max range skew
## kokugo 1 21 49.43 7.01 50 49.76 7.41 34 61 27 -0.35
## sansuu 2 21 73.81 2.02 74 73.94 1.48 69 77 8 -0.52
## rika 3 21 53.43 3.98 53 53.71 4.45 43 60 17 -0.64
## syakai 4 21 39.48 10.19 44 39.82 11.86 20 54 34 -0.30
## eigo 5 21 44.76 11.12 45 44.82 11.86 25 65 40 -0.10
## kurtosis se
## kokugo -0.78 1.53
## sansuu -0.03 0.44
## rika 0.12 0.87
## syakai -1.29 2.22
## eigo -0.94 2.43
## --------------------------------------------------------
## group: 2
## var n mean sd median trimmed mad min max range skew kurtosis
## kokugo 1 21 81.95 6.55 81 81.82 8.90 72 94 22 0.16 -1.26
## sansuu 2 21 68.95 5.24 69 68.71 5.93 60 79 19 0.46 -0.79
## rika 3 21 47.62 4.18 46 47.47 4.45 41 55 14 0.33 -1.30
## syakai 4 21 60.71 12.97 59 60.65 10.38 39 86 47 0.20 -0.94
## eigo 5 21 78.05 9.10 78 78.12 8.90 61 94 33 0.10 -0.93
## se
## kokugo 1.43
## sansuu 1.14
## rika 0.91
## syakai 2.83
## eigo 1.99
## --------------------------------------------------------
## group: 3
## var n mean sd median trimmed mad min max range skew
## kokugo 1 27 65.85 8.44 66 66.35 7.41 44 80 36 -0.55
## sansuu 2 27 77.56 3.57 77 77.35 2.97 72 86 14 0.49
## rika 3 27 53.93 4.85 53 53.65 5.93 47 66 19 0.44
## syakai 4 27 51.81 10.68 51 51.96 11.86 23 73 50 -0.29
## eigo 5 27 60.70 11.26 63 60.96 11.86 38 82 44 -0.38
## kurtosis se
## kokugo -0.22 1.62
## sansuu -0.48 0.69
## rika -0.45 0.93
## syakai 0.22 2.05
## eigo -0.87 2.17
## --------------------------------------------------------
## group: 4
## var n mean sd median trimmed mad min max range skew
## kokugo 1 30 61.57 7.00 61.0 61.58 8.15 46 75 29 -0.07
## sansuu 2 30 66.33 4.09 68.0 66.54 3.71 58 75 17 -0.36
## rika 3 30 47.30 5.36 48.5 47.62 5.93 34 55 21 -0.53
## syakai 4 31 46.48 10.40 47.0 45.96 13.34 30 71 41 0.36
## eigo 5 31 57.10 9.50 56.0 56.48 7.41 41 80 39 0.61
## kurtosis se
## kokugo -0.68 1.28
## sansuu -0.51 0.75
## rika -0.71 0.98
## syakai -0.72 1.87
## eigo -0.06 1.71
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 = 0.1645, df = 3, p-value = 0.9831
# クラスターごとに検定したりして。
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.1.
## It was executed on Mon Dec 17 14:47:35 2012.
##
##
## << DESCRIPTIVE STATISTICS >>
##
## --------------------------------------------
## A N Mean S.D.
## --------------------------------------------
## 1 21 49.4286 7.0112
## 3 27 65.8519 8.4385
## 2 21 81.9524 6.5458
## 4 30 61.5667 7.0009
## --------------------------------------------
##
##
## << 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 11473.8580 3 3824.6193 71.0636 0.0000 *** 0.6917
## Error 5112.8693 95 53.8197
## ----------------------------------------------------------------------------------------------------
## 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.
## --------------------------------------------
## 1 21 49.4286 7.0112
## 3 27 65.8519 8.4385
## 2 21 81.9524 6.5458
## 4 30 61.5667 7.0009
## --------------------------------------------
##
## ---------------------------------------------------------------------------
## Pair Interval t-value df p adj.p
## ---------------------------------------------------------------------------
## 1-2 -32.5238 14.3657 95 0.0000 0.0000 1 < 2 *
## 2-4 20.3857 9.7665 95 0.0000 0.0000 2 > 4 *
## 1-3 -16.4233 7.6941 95 0.0000 0.0000 1 < 3 *
## 3-2 -16.1005 7.5429 95 0.0000 0.0000 3 < 2 *
## 1-4 -12.1381 5.8152 95 0.0000 0.0000 1 < 4 *
## 3-4 4.2852 2.2019 95 0.0301 0.0301 3 > 4 *
## ---------------------------------------------------------------------------
##
##
## output is over --------------------///
課題は，http://www1.doshisha.ac.jp/~mjin/data/ にある職業適性能力検査データをつかって，因子分析をし（抽出方法の指定，回転の指定もする，因子得点も算出）クラスター分析してください。