主成分得点のデータセットへの追加

FactoMineRの“PCA”コマンドは主成分得点を散布図に描いてくれるので大変便利であある。しかし、より掘り下げた分析を進めるには、やはり各都道府県の主成分得点をデータセットに追加しておきたい。それには以下のようにする。

まず主成分分析の結果を保存している“pca”の中身を確認する。

names(pca)

## [1] "eig"  "var"  "ind"  "svd"  "call"

色々な結果が保存されていることが分かるが、散布図を描くときに参照された座標値は結果“pca＄ind＄coord”の第１列“Dim.1”と第２列“Dim.2”であることが確認できる。実際、これらをそれぞれ“PC1”、“PC2”と名前をつけてデータセット“iryou.std”に追加し、散布図を描いてみると上と同じ散布図が得られる。

iryou.std$PC1 <- pca$ind$coord[,1]
iryou.std$PC2 <- pca$ind$coord[,2]
ggplot(iryou.std, aes(x=PC1,y=PC2)) +
  geom_point() +
  geom_text(aes(label=rownames(iryou.std)))

クラスタリング結果のランダム性

上で得られた二つの主成分（PC1とPC2）に基づいて、都道府県別医療状況のクラスター分析を行う。

OBSの授業「ビジネス統計分析」でも説明したとおり、クラスター分析には非階層型クラスタリングと階層型クラスタリングの二つの方法がある。

ここでは授業を復習する意味合いから、まず非階層型クラスタリングを行う。非階層型クラスタリングで大切な点は、第１にクラスタ数を指定しなければならないこと。第２はクラスターの種として最初にランダムにとられる何個かの点の選び方によって結果として得られるクラスター構成に違いがもたらされる、つまり実行するたびに異なった結果になる可能性があるということである。

都道府県を幾つのクラスターに分けるかだが、散布図をみて自然にクラスター数が見て取れるわけではない。こんな場合は、まずは上下左右に分けてみるという意味合いからクラスター数を４、あるいは３と指定して、その後でクラスターごとの特徴の違いを調べるという手順をとると色々なことが分かることが多い。

clst<-kmeans(iryou.std[,c("PC1","PC2")],centers=4,iter.max=100)

結果の一覧は次のとおりである。

names(clst)

## [1] "cluster"      "centers"      "totss"        "withinss"    
## [5] "tot.withinss" "betweenss"    "size"         "iter"        
## [9] "ifault"

length(clst$cluster)

## [1] 47

結果の中の“cluster”がグループ分けされた後、各都道府県に割り当てられたグループ番号である。これを元のデータ表に“cluster”という名前で追加してから、四つのクラスターで層別化して都道府県の散布図を描く。

iryou.std$cluster <- clst$cluster
p.cluster <- ggplot(iryou.std,aes(x=PC1,y=PC2,color=as.factor(cluster)))
p.cluster <- p.cluster + geom_point()
p.cluster <- p.cluster + geom_text(aes(label=rownames(iryou.std),vjust=1))
p.cluster <- p.cluster+xlim(-8.5,8.5)
p.cluster <- p.cluster+theme_bw()
plot(p.cluster)

クラスター分析の結果にランダム性が残っていることを確認するため同じクラスタリングをもう一度行ってみる。

clst.again <- kmeans(iryou.std[,c("PC1","PC2")],centers=4,iter.max=100)

２番目のクラスタリングに基づいて散布図を描く。

iryou.std$cluster.again <- clst.again$cluster
p.cluster <- ggplot(iryou.std,aes(x=PC1,y=PC2,color=as.factor(cluster.again)))
p.cluster <- p.cluster + geom_point()
p.cluster <- p.cluster + geom_text(aes(label=rownames(iryou.std),vjust=1))
p.cluster <- p.cluster+xlim(-8.5,8.5)
p.cluster <- p.cluster+theme_bw()
plot(p.cluster)

１番目と２番目の散布図を見比べると、左上の区域にな東京都と大阪府が、最初のクラスタリングでは同じグループに入っているのに、２番目の散布図では別々のグループに分かれている。

各グループとクラスター番号との対応や散布図におけるクラスター番号の位置取りも二つの散布図で違いがある。

このように一度きりのクラスタリングで４７都道府県のグループ分けを確定するのは問題が残る。また、クラスター数をここではアドホックに４と指定したが、なぜクラスター数は４でなければならないか、その最適性にも説明が必要である。

最適なクラスター数

最適なクラスター数を決めるには、クラスター内の分散に着目する。理想的なクラスタリングは、現実には発生しないにせよ、各クラスターの内部においては全て同一の値、クラスター間にのみ違いあるという場合である。したがってクラスター内の分散が小さくなるクラスター数の方がよいと判定するのは自然である。実際、クラスター数を増やしていくと、クラスター内のデータはまとまってくる傾向がある。クラスター内の分散を合計した値が減少する速さに着目して、その減少速度が遅くなる点でクラスター数を決める方式は合理性がある。

更に、クラスター内分散に対するクラスター間分散の比率に着目してもよい。これはCalinski-Harabasz(CH) Indexである。理想的なクラスタリングにおいては、クラスター内分散がゼロ、クラスター間でのみ分散が残るのでCH指標は無限大になる。通常のケースではなるべくCH指標がピークをとるクラスター数を選べばよい。

まず所要の関数を定義しておく。

sqr_edist <- function(x, y){
  sum((x-y)^2)
}

wss.cluster <- function(clustermat){
  c0 <- apply(clustermat, 2, FUN=mean)
  sum(apply(clustermat, 1, FUN=function(row){sqr_edist(row,c0)}))
}

wss.total <- function(dmatrix, labels){
  wsstot <- 0
  k <- length(unique(labels))
  for(i in 1:k){
    wsstot <- wsstot + wss.cluster(subset(dmatrix, labels==i))
  }
  wsstot
}

totss <- function(dmatrix){
  grandmean <- apply(dmatrix,2,FUN=mean)
  sum(apply(dmatrix, 1, FUN=function(row){sqr_edist(row, grandmean)}))
}
ch_criterion <- function(dmatrix, kmax, method="kmeans"){
  if(!(method %in% c("kmeans", "hclust"))){
    stop("method must be one of c('kmeans','hclust')")
  }
  npts <- dim(dmatrix)[1] # number of rows
  totss <- totss(dmatrix)
  wss <- numeric(kmax)
  crit <- numeric(kmax)
  
  wss[1] <- (npts-1)*sum(apply(dmatrix,2,var))
  for(k in 2:kmax){
    if(method=="kmeans"){
      clustering <- kmeans(dmatrix, k, nstart=10, iter.max = 100)
      wss[k] <- clustering$tot.withinss
    } else {
      d <- dist(dmatrix, method="euclidean")
      pfit <- hclust(d, method="ward")
      labels <- cutree(pfit, k=k)
      wss[k] <- wss.total(dmatrix, labels)
    }
  }
  bss <- totss - wss
  crit.num <- bss/(0:(kmax-1))
  crit.denom <- wss/(npts - 1:kmax)
  crit <- crit.num / crit.denom
  crit[1] <- NA
  list(crit=crit, wss=wss, totss=totss)
}

クラスター数を１個（＝クラスタリングしない）から１０個までに範囲を決め、順に主成分得点“PC1”と“PC2”に基づくkmeans“を実行しよう。このとき、クラスターごとにクラスター内分散とCH Indexを保存する。そうしてから横軸にクラスター数を、縦軸にクラスター内分散、CH Indexをとって折れ線グラフを描く。

pcs <- c("PC1", "PC2")
clustcrit <- ch_criterion(iryou.std[,pcs], 10, method="kmeans")
critframe <- data.frame(k=1:10, ch=scale(clustcrit$crit), wss=scale(clustcrit$wss))
critframe2 <- melt(critframe, id.vars = c("k"), variable.name = "measure", value.name = "score")
ggplot(critframe2, aes(x=k, y=score, color=measure)) +
  geom_point(aes(shape=measure)) +
  geom_line(aes(linetype=measure)) +
  scale_x_continuous(breaks = 1:10, labels = 1:10)

この図をみると、クラスター数を増やすに伴ってクラスター内分散は減少している。一方、CH指標はクラスター数が３の前後でピークを形成している。６個以降はCH指標が上昇しているが、これはクラスター内分散の減少速度が鈍くなる一方でクラスター数の増加によってクラスター間の分散が増えていることよる。

最初に行ったクラスタリングではクラスター数を４としたが、最適なクラスター数としては３のほうが良いと判断できる。

安定したクラスタリング結果

上ではクラスタリング結果には常にランダム性が残ることに注意を促した。実際、都道府県の医療状況に着目するとしても、現実に他の都道府県と明確な違いがあり、クラスターを設けることの根拠がある現実のクラスターがある一方で、クラスター数を４と指定したことが原因となって実際には存在しないクラスターが計算上は導かれてくる、そんな見せかけのクラスターもまた得られがちである。実態の伴わないクラスタリングで形成されたグループは安定性に欠け、クラスター分析を行うたびに所属するグループが変わったりする。

パッケージ“fpc”にある“clusterboot”は何度も反復して計算をおこない最終的に安定したクラスタリングを返してくれるので大変便利である。

ここではクラスタリングを100回反復して各グループの安定性をみてみる。

cboot <- clusterboot(iryou.std[,c("PC1","PC2")], clustermethod = kmeansCBI, krange=3:4)

最数的に確定できたグループを表示させよう。

group <- cboot$result$partition
for(i in 1:3){
  print(paste("cluster", i))
  print(iryou$Pref[group==i])
}

## [1] "cluster 1"
##  [1] Hokkaido  Toyama    Ishikawa  Wakayama  Tottori   Shimane   Hiroshima
##  [8] Yamaguchi Tokushima Kagawa    Ehime     Kochi     Fukuoka   Saga     
## [15] Nagasaki  Kumamoto  Oita      Miyazaki  Kagoshima
## 47 Levels: Aichi Akita Aomori Chiba Ehime Fukui Fukuoka Fukushima ... Yamanashi
## [1] "cluster 2"
##  [1] Miyagi   Tochigi  Saitama  Chiba    Tokyo    Kanagawa Fukui   
##  [8] Gifu     Shizuoka Aichi    Mie      Shiga    Kyoto    Osaka   
## [15] Hyogo    Nara     Okayama  Okinawa 
## 47 Levels: Aichi Akita Aomori Chiba Ehime Fukui Fukuoka Fukushima ... Yamanashi
## [1] "cluster 3"
##  [1] Aomori    Iwate     Akita     Yamagata  Fukushima Ibaraki   Gumma    
##  [8] Niigata   Yamanashi Nagano   
## 47 Levels: Aichi Akita Aomori Chiba Ehime Fukui Fukuoka Fukushima ... Yamanashi

これらのグループはクラスター計算のたびにどの程度安定して同じ構成員と組み合わされたのだろうか。

cboot$bootmean

## [1] 0.8165997 0.7997885 0.8132880

表示された数値は、いわゆる“Jaccard Coefficient”（ジャッカード係数)を100回にわたるクラスタリング計算を通して平均した値であって、ある回のクラスター計算による構成メンバーが別の回の構成メンバーとどの程度まで共通しているかという指標の平均値である。当然ながら、クラスター計算によらず同じ構成員が含まれるグループは安定性が高いわけである。上の数値をみると、どれも７５％を超えており、十分な安定性を有している。クラスターの安定性を示す平均ジャッカード係数が0.5未満になる場合は、そのクラスターに実質的な意味はないと判断される。

上の結果に基づいて層別化された散布図にしてみる。

iryou.std$cluster.final <- cboot$result$partition
p.cluster <- ggplot(iryou.std,aes(x=PC1,y=PC2,color=as.factor(cluster.final)))
p.cluster <- p.cluster + geom_point()
p.cluster <- p.cluster + geom_text(aes(label=rownames(iryou.std),vjust=1))
p.cluster <- p.cluster+xlim(-8.5,8.5)
p.cluster <- p.cluster+theme_bw()
plot(p.cluster)

医療施設、医師、看護師等の数

• D1: 一般病院数（10万人当たり）

iryou %>% group_by(cluster) %>% summarise(avgD1=mean(D1))

## # A tibble: 3 x 2
##   cluster avgD1
##     <int> <dbl>
## 1       1  9.61
## 2       2  5.01
## 3       3  5.58

• D2: 一般診療所数（10万人当たり）

iryou %>% group_by(cluster) %>% summarise(avgD2=mean(D2))

## # A tibble: 3 x 2
##   cluster avgD2
##     <int> <dbl>
## 1       1  86.2
## 2       2  76.6
## 3       3  72.7

• D3: 一般病院病床数（10万人当たり）

iryou %>% group_by(cluster) %>% summarise(avgD3=mean(D3))

## # A tibble: 3 x 2
##   cluster avgD3
##     <int> <dbl>
## 1       1 1434.
## 2       2  949.
## 3       3 1059.

• D4: 医療施設に従事する医師数（10万人当たり）

iryou %>% group_by(cluster) %>% summarise(avgD4=mean(D4))

## # A tibble: 3 x 2
##   cluster avgD4
##     <int> <dbl>
## 1       1  258.
## 2       2  220 
## 3       3  196.

• D5: 医療施設に従事する（准）看護師数（10万人当たり）

iryou %>% group_by(cluster) %>% summarise(avgD5=mean(D5))

## # A tibble: 3 x 2
##   cluster avgD5
##     <int> <dbl>
## 1       1 1259.
## 2       2  829.
## 3       3  920.

• D6: 救急自動車数（10万人当たり）

iryou %>% group_by(cluster) %>% summarise(avgD6=mean(D6))

## # A tibble: 3 x 2
##   cluster avgD6
##     <int> <dbl>
## 1       1  6.51
## 2       2  4.68
## 3       3  6.9

• D7: 年間救急出動件数（千人当たり）

iryou %>% group_by(cluster) %>% summarise(avgD7=mean(D7))

## # A tibble: 3 x 2
##   cluster avgD7
##     <int> <dbl>
## 1       1  42.6
## 2       2  44.7
## 3       3  39.5

• D8: 薬局数（10万人当たり）

iryou %>% group_by(cluster) %>% summarise(avgD8=mean(D8))

## # A tibble: 3 x 2
##   cluster avgD8
##     <int> <dbl>
## 1       1  48.5
## 2       2  40.7
## 3       3  45.0

患者数、利用率等

• D9: 一般病院外来患者数（常勤医師1人1日当たり）

iryou %>% group_by(cluster) %>% summarise(avgD9=mean(D9))

## # A tibble: 3 x 2
##   cluster avgD9
##     <int> <dbl>
## 1       1  8.86
## 2       2  8.35
## 3       3  9.79

• D10: 一般病院在院患者数（常勤医師1人1日当たり）

iryou %>% group_by(cluster) %>% summarise(avgD10=mean(D10))

## # A tibble: 3 x 2
##   cluster avgD10
##     <int>  <dbl>
## 1       1   8.39
## 2       2   6.16
## 3       3   7.58

• D11: 一般病院在院患者数（看護師・准看護士1人1日当たり）

iryou %>% group_by(cluster) %>% summarise(avgD11=mean(D11))

## # A tibble: 3 x 2
##   cluster avgD11
##     <int>  <dbl>
## 1       1   1.35
## 2       2   1.24
## 3       3   1.29

• D12: 一般病院病床利用率

iryou %>% group_by(cluster) %>% summarise(avgD12=mean(D12))

## # A tibble: 3 x 2
##   cluster avgD12
##     <int>  <dbl>
## 1       1   82.1
## 2       2   79.1
## 3       3   77.6

• D13: 一般病院平均在院日数

iryou %>% group_by(cluster) %>% summarise(avgD13=mean(D13))

## # A tibble: 3 x 2
##   cluster avgD13
##     <int>  <dbl>
## 1       1   31.5
## 2       2   24.5
## 3       3   25.3

• D14: 一般病院年間新入院患者数

iryou %>% group_by(cluster) %>% summarise(avgD14=mean(D14))

## # A tibble: 3 x 2
##   cluster avgD14
##     <int>  <dbl>
## 1       1 13678.
## 2       2 11178.
## 3       3 11889.

病因別死亡率等

• D15: 標準化死亡率（千人当たり）

iryou %>% group_by(cluster) %>% summarise(avgD15=mean(D15))

## # A tibble: 3 x 2
##   cluster avgD15
##     <int>  <dbl>
## 1       1   1.90
## 2       2   1.83
## 3       3   1.94

• D16: 生活習慣病による死亡者数（10万人当たり）

iryou %>% group_by(cluster) %>% summarise(avgD16=mean(D16))

## # A tibble: 3 x 2
##   cluster avgD16
##     <int>  <dbl>
## 1       1   648.
## 2       2   531.
## 3       3   676.

• D17: 悪性新生物による死亡者数（10万人当たり）

iryou %>% group_by(cluster) %>% summarise(avgD17=mean(D17))

## # A tibble: 3 x 2
##   cluster avgD17
##     <int>  <dbl>
## 1       1   327.
## 2       2   276.
## 3       3   329.

• D18: 糖尿病による死亡者数（10万人当たり）

iryou %>% group_by(cluster) %>% summarise(avgD18=mean(D18))

## # A tibble: 3 x 2
##   cluster avgD18
##     <int>  <dbl>
## 1       1   12.6
## 2       2   10.5
## 3       3   13.9

• D19: 高血圧疾患による死亡者数（10万人当たり）

iryou %>% group_by(cluster) %>% summarise(avgD19=mean(D19))

## # A tibble: 3 x 2
##   cluster avgD19
##     <int>  <dbl>
## 1       1   6.18
## 2       2   5.29
## 3       3   6.75

• D20: 心疾患（高血圧を除く）による死亡者数（10万人当たり）

iryou %>% group_by(cluster) %>% summarise(avgD20=mean(D20))

## # A tibble: 3 x 2
##   cluster avgD20
##     <int>  <dbl>
## 1       1   190.
## 2       2   151.
## 3       3   188.

• D21: 脳血管疾患による死亡者数（10万人当たり）

iryou %>% group_by(cluster) %>% summarise(avgD21=mean(D21))

## # A tibble: 3 x 2
##   cluster avgD21
##     <int>  <dbl>
## 1       1  112. 
## 2       2   87.9
## 3       3  138.

• D22: 自殺者数（10万人当たり）

iryou %>% group_by(cluster) %>% summarise(avgD22=mean(D22))

## # A tibble: 3 x 2
##   cluster avgD22
##     <int>  <dbl>
## 1       1   21.1
## 2       2   19.8
## 3       3   23.8

distance <- dist(iryou.std[,vars],method="euclidean")

fit <- hclust(distance, method="ward")
plot(fit)

group <- cutree(fit, k=3)
for(i in 1:3){
  print(paste("cluster", i))
  print(iryou$Pref[group==i])
}

## [1] "cluster 1"
##  [1] Hokkaido  Toyama    Ishikawa  Fukui     Kyoto     Wakayama  Tottori  
##  [8] Shimane   Okayama   Hiroshima Yamaguchi Tokushima Kagawa    Ehime    
## [15] Kochi     Fukuoka   Saga      Nagasaki  Kumamoto  Oita      Miyazaki 
## [22] Kagoshima
## 47 Levels: Aichi Akita Aomori Chiba Ehime Fukui Fukuoka Fukushima ... Yamanashi
## [1] "cluster 2"
##  [1] Aomori    Iwate     Akita     Yamagata  Fukushima Ibaraki   Gumma    
##  [8] Niigata   Yamanashi Nagano   
## 47 Levels: Aichi Akita Aomori Chiba Ehime Fukui Fukuoka Fukushima ... Yamanashi
## [1] "cluster 3"
##  [1] Miyagi   Tochigi  Saitama  Chiba    Tokyo    Kanagawa Gifu    
##  [8] Shizuoka Aichi    Mie      Shiga    Osaka    Hyogo    Nara    
## [15] Okinawa 
## 47 Levels: Aichi Akita Aomori Chiba Ehime Fukui Fukuoka Fukushima ... Yamanashi

Zumel, N. and J. Mount “Practical Science with R”, Ch.8, 2014