1 非階層型クラスター分析

1.1 ファイルの読み込み

# CSVファイルの読み込み（ヘッダーと列名がある場合）
dat <- read.csv("data/data_ch11-1.csv", header = TRUE, row.names = 1)
#row.namesは行名の指定。row.names=1だと1列目が列名になる
head(dat)
Kokugo Eigo Syakai Sugaku Rika
S001 85 92 80 70 74
S002 54 60 69 99 80
S003 94 76 71 64 63
S004 76 95 80 68 61
S005 63 68 59 82 97
S006 72 78 70 95 93 
# 各科目の平均点（＝列平均）を計算
colMeans(dat)
## Kokugo   Eigo Syakai Sugaku   Rika 
##  73.20  76.60  70.45  77.20  75.95
# 各学習者の平均点（＝行平均）を計算
#rowMeans(dat)

1.2 ユークリッド距離

# ユークリッド距離の計算
dat.d <- dist(dat)
# ユークリッド距離の計算結果の確認
#dat.d

1.3 ウォード法

# ウォード法によるクラスターの作成
dat.hc <- hclust(dat.d, method = "ward.D2")
summary(dat.hc)
##             Length Class  Mode     
## merge       38     -none- numeric  
## height      19     -none- numeric  
## order       20     -none- numeric  
## labels      20     -none- character
## method       1     -none- character
## call         3     -none- call     
## dist.method  1     -none- character
dat.hc$merge#結合
##       [,1] [,2]
##  [1,]  -18  -19
##  [2,]   -2  -10
##  [3,]  -11  -12
##  [4,]   -5    1
##  [5,]   -8  -14
##  [6,]   -4   -7
##  [7,]  -15  -16
##  [8,]  -20    5
##  [9,]   -6   -9
## [10,]   -3    3
## [11,]  -17    7
## [12,]    6   11
## [13,]   -1   10
## [14,]    4    8
## [15,]  -13   13
## [16,]    9   14
## [17,]    2   16
## [18,]   12   15
## [19,]   17   18
#cophenetic(dat.hc) #コーフェン行列
#（コーフェン行列は、属するクラスター間の距離を表す）

階層的クラスターを形成する階層の情報である。マイナス符号が付いているのが個体の番号、

マイナスが付いていないのがクラスターの番号である。行番号がクラスター形成の順番である。 http://hyoccori-toukei.com/id-3/id/2019.html.

1.4 樹形図

# 樹形図による可視化
#plot(dat.hc)

# 樹形図における個々のデータを一番下にそろえて配置
plot(dat.hc, hang = -1)

# 個々のデータが含まれるクラスターの確認
cutree(dat.hc, k = 2)#（クラスター数が2の場合）
## S001 S002 S003 S004 S005 S006 S007 S008 S009 S010 S011 S012 S013 S014 S015 S016 
##    1    2    1    1    2    2    1    2    2    2    1    1    1    2    1    1 
## S017 S018 S019 S020 
##    1    2    2    2
cutree(dat.hc, k = 4)#（クラスター数が4の場合）
## S001 S002 S003 S004 S005 S006 S007 S008 S009 S010 S011 S012 S013 S014 S015 S016 
##    1    2    1    3    4    4    3    4    4    2    1    1    1    4    3    3 
## S017 S018 S019 S020 
##    3    4    4    4
# 樹形図による可視化
plot(dat.hc, hang = -1)
# 樹形図におけるクラスターを囲む四角形を表示（クラスター数が2の場合）
rect.hclust(dat.hc, k = 2)
rect.hclust(dat.hc, k = 4)

1.5 ヒートマップ

# 分析データを行列に変換
dat.2 <- as.matrix(dat)
# 分析データの形式を再確認
class(dat.2)
## [1] "matrix" "array"
dat.2
##      Kokugo Eigo Syakai Sugaku Rika
## S001     85   92     80     70   74
## S002     54   60     69     99   80
## S003     94   76     71     64   63
## S004     76   95     80     68   61
## S005     63   68     59     82   97
## S006     72   78     70     95   93
## S007     71   98     72     69   60
## S008     68   66     58     95   88
## S009     74   70     74     86   95
## S010     57   64     66    100   82
## S011     88   83     70     66   64
## S012     82   83     69     65   66
## S013     91   77     90     55   69
## S014     65   64     58     90   84
## S015     85   87     77     64   56
## S016     77   89     75     57   54
## S017     74   91     85     60   57
## S018     61   69     62     87   93
## S019     60   67     64     83   96
## S020     67   55     60     89   87
# ヒートマップつきの樹形図の作成（ユークリッド距離，ウォード法）
heatmap(dat.2, cexCol = 0.5, hclustfun = function(x){hclust(x, method = "ward.D2")})

2 非階層型クラスター分析

2.1 k-means法　

# 乱数の種を固定
set.seed(1)
# k-means法（クラスター数は2）
dat.km <- kmeans(dat, center = 2)
dat.km
## K-means clustering with 2 clusters of sizes 10, 10
## 
## Cluster means:
##   Kokugo Eigo Syakai Sugaku Rika
## 1   82.3 87.1   76.9   63.8 62.4
## 2   64.1 66.1   64.0   90.6 89.5
## 
## Clustering vector:
## S001 S002 S003 S004 S005 S006 S007 S008 S009 S010 S011 S012 S013 S014 S015 S016 
##    1    2    1    1    2    2    1    2    2    2    1    1    1    2    1    1 
## S017 S018 S019 S020 
##    1    2    2    2 
## 
## Within cluster sum of squares by cluster:
## [1] 2005.9 1690.7
##  (between_SS / total_SS =  76.4 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"
# k-means法によるクラスタリング結果の確認
dat.km$cluster
## S001 S002 S003 S004 S005 S006 S007 S008 S009 S010 S011 S012 S013 S014 S015 S016 
##    1    2    1    1    2    2    1    2    2    2    1    1    1    2    1    1 
## S017 S018 S019 S020 
##    1    2    2    2
# パッケージの読み込み
library("cluster")
# k-means法によるクラスタリング結果の可視化
clusplot(dat, dat.km$cluster)

# k-means法によるクラスタリング結果の可視化（個々の学習者のIDを表示）
clusplot(dat, dat.km$cluster, labels = 2)

# Gap統計量の計算
set.seed(1)
dat.Gap <- clusGap(dat, kmeans, K.max = 5)
# 計算されたGap統計量の確認
dat.Gap
## Clustering Gap statistic ["clusGap"] from call:
## clusGap(x = dat, FUNcluster = kmeans, K.max = 5)
## B=100 simulated reference sets, k = 1..5; spaceH0="scaledPCA"
##  --> Number of clusters (method 'firstSEmax', SE.factor=1): 2
##          logW   E.logW         gap     SE.sim
## [1,] 5.157829 4.929641 -0.22818770 0.06221962
## [2,] 4.460194 4.586987  0.12679330 0.05167384
## [3,] 4.307199 4.421919  0.11471922 0.05583342
## [4,] 4.198610 4.267757  0.06914638 0.05672929
## [5,] 3.985888 4.134093  0.14820484 0.06046819
# Gap統計量とクラスター数の関係の可視化
plot(dat.Gap)

Gap統計量とは、ランダムな点を用いた残差平方和から、実際にクラスタリングを行った結果の残差平方和を引いたものなので、 「Gap統計量が大きい＝正しい

Gap 統計量とは, クラスタ数kのクラスタリング結果の評価として, クラスタのコンパクトト 性を表す指標 W𝑘を用い, 与えられたデータセットPと, 一様乱数分布の(つまり, クラスタを持たない) B 個のデータセットPb(b = 1,…,B)の W𝑘を比較することで 最適なクラス タ数を推定する手法である. Bはブートストラップ回数といい, 後の計算において, このB個のデータセットの平均をとることで,与えられたデータセットとの対照としてより精度の高い比較を実現する. それゆえ, この値が高いほど精度が向上する. 応用にもよるが, 一般的には 100 前後で用いられ, 500 程度で極めて高い正確性をもつ

https://proceedings-of-deim.github.io/DEIM2020/papers/D1-5.pdf

3 決定木

# CSVファイルの読み込み（ヘッダーがある場合）
dat <- read.csv("data/data_ch11-2.csv", header = TRUE)
# 読み込んだデータの冒頭の確認
head(dat)
Listening Reading Writing Speaking Grade
62 88 61 79 pass
67 60 90 67 pass
73 57 79 47 pass
72 54 57 78 fail
69 67 78 81 pass
91 66 60 65 pass 
# 読み込んだデータ（全体）の概要の確認
summary(dat)
##    Listening        Reading         Writing         Speaking    
##  Min.   :52.00   Min.   :50.00   Min.   :56.00   Min.   :40.00  
##  1st Qu.:59.00   1st Qu.:60.00   1st Qu.:62.75   1st Qu.:49.75  
##  Median :67.00   Median :68.00   Median :70.00   Median :62.50  
##  Mean   :68.03   Mean   :69.74   Mean   :70.92   Mean   :62.22  
##  3rd Qu.:73.00   3rd Qu.:74.50   3rd Qu.:77.25   3rd Qu.:73.25  
##  Max.   :99.00   Max.   :96.00   Max.   :90.00   Max.   :85.00  
##     Grade          
##  Length:100        
##  Class :character  
##  Mode  :character  
##                    
##                    
## 
# テストの合格者数と不合格者数を確認
table(dat$Grade)
## 
## fail pass 
##   20   80
library("rpart")
# 決定木分析（Gradeを結果変数，それ以外を予測変数として指定）
rpart.result <- rpart(Grade ~ ., data = dat)
# 分析結果の確認
rpart.result
## n= 100 
## 
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
## 
## 1) root 100 20 pass (0.20000000 0.80000000)  
##   2) Reading< 61.5 28 13 fail (0.53571429 0.46428571)  
##     4) Writing< 75.5 16  1 fail (0.93750000 0.06250000) *
##     5) Writing>=75.5 12  0 pass (0.00000000 1.00000000) *
##   3) Reading>=61.5 72  5 pass (0.06944444 0.93055556) *
#install.packages("partykit", dependencies = TRUE)
# パッケージの読み込み
library("partykit")
## Loading required package: grid
## Loading required package: libcoin
## Loading required package: mvtnorm
# 決定木の可視化
plot(as.party(rpart.result))