機械学習手習い: クラスタリング

「入門 機械学習」手習い、9日目。「9章 MDS:米国上院議員の類似度の視覚的な調査」です。

www.amazon.co.jp

観測値をクラスタリングするための、多次元尺度構成法(MDS:multidimensional scaling)を学び、得票に基づいて米国上院議員をクラスタリングしてみます。

# 前準備
> setwd("09-MDS/")

多次元尺度構成法でのクラスタリング

最初に、顧客の製品評価データを使って顧客を分類する例を試し、多次元尺度構成法の考え方を学びます。

まずは、製品評価データを作成します。

> set.seed(851982)
> ex.matrix <- matrix(sample(c(-1, 0, 1), 24, replace = TRUE), nrow = 4, ncol = 6)
> row.names(ex.matrix) <- c('A', 'B', 'C', 'D')
> colnames(ex.matrix) <- c('P1', 'P2', 'P3', 'P4', 'P5', 'P6')
> ex.matrix
  P1 P2 P3 P4 P5 P6
A  0 -1  0 -1  0  0
B -1  0  1  1  1  0
C  0  0  0  1 -1  1
D  1  0  1 -1  0  0

列(P1,P2..)が製品、行(A,B)が顧客を示します。値が評価値で、1が良い、-1が悪い、0が未評価を意味します。 評価値はsampleを使ってランダムに生成しています。

このデータから顧客同士の関係を取り出すため、これを関係を要約する数値に変換する必要があります。 具体的には、以下のような手順で、顧客同士の意見の一致度を示す行列にします。

• 1) 各製品に対するAさんの評価、Bさんの評価を比較し、算出された数値を合計
  • 両方とも1 or -1 なら、意見が一致(+1)
  • 片方が1で、もう一方が-1の場合、意見は不一致(-1)
  • いずれか一方が0の場合、未評価(0)
• 2) この評価を全ユーザーの組み合わせで行い、行列を作る。

例) AさんとBさんの一致度

一致度 = 製品1の一致度 + 製品2の一致度 + 製品3の一致度 + 製品4の一致度 + 製品5の一致度 + 製品6の一致度
       = (0 * -1) + (-1 * 0) + (0 * 1) + (-1 * 1) + (0 * 1) + (0 * 0)
       = -1

この操作は、行列を転置して掛け合わせることでさくっと計算できます。

> ex.mult <- ex.matrix %*% t(ex.matrix)
> ex.mult
   A  B  C  D
A  2 -1 -1  1
B -1  4  0 -1
C -1  0  3 -1
D  1 -1 -1  3

各値が、↑で説明した意見の一致度を示す数値になります。 ちなみに、AxAやBxBの場合、意見は完全一致なので、高い値になります。

次に、このデータが示す顧客間の距離を、ユークリッド距離を計算することで数値化します。 ユークリッド距離、↑のデータを4次元座標にマッピングしたときの、各座標間の距離です。(という認識であってるのかな?) ユークリッド距離は、 dist 関数で計算できます。

> ex.dist <- dist(ex.mult)
> ex.dist
         A        B        C
B 6.244998                  
C 5.477226 5.000000         
D 2.236068 6.782330 6.082763

A,Dは距離が近く、意見が似ているとわかります。

最後に、この距離行列を多次元尺度構成法で可視化します。

今回使う多次元尺度構成法では、データ中のすべての点間の距離を示す距離行列を使い、その2点間の距離を近似する座標を返します。

これを使って、顧客の関係を2次元の図にしてみます。

> png("plot1.png", width = 100, height = 100)
> ex.mds <- cmdscale(ex.dist)
> plot(ex.mds, type = 'n')
> text(ex.mds, c('A', 'B', 'C', 'D'))
> dev.off()

散布図から、A,Dは類似性があることがわかります。

米国上院議員のクラスタリング

多次元尺度構成法の概要がつかめたところで、本題に取り組みます。 法案への賛成/反対データを使って、米国上院議員を分類してみます。

まずは、データの読み込みとクリーニングから。

> library('foreign')
> library('ggplot2')

> data.dir <- file.path("data", "roll_call")
> data.files <- list.files(data.dir)

# 全データを読み込む
> rollcall.data <- lapply(data.files, function(f) {
  read.dta(file.path(data.dir, f), convert.factors = FALSE)
})

今回は、議員の名前、所属政党、賛否のみ使うので、これを取り出します。

> rollcall.simplified <- function(df) {
  no.pres <- subset(df, state < 99) # state 99は副大統領。今回はこれは除外する。
  
  # 投票方法にはいくつかあるので、賛成(1)/反対(-1)/非投票(0)に丸める
  for(i in 10:ncol(no.pres)) {
    no.pres[,i] <- ifelse(no.pres[,i] > 6, 0, no.pres[,i])
    no.pres[,i] <- ifelse(no.pres[,i] > 0 & no.pres[,i] < 4, 1, no.pres[,i])
    no.pres[,i] <- ifelse(no.pres[,i] > 1, -1, no.pres[,i])
  }
  
  return(as.matrix(no.pres[,10:ncol(no.pres)]))
}

> rollcall.simple <- lapply(rollcall.data, rollcall.simplified)

元データができたので、距離行列を作成。

> rollcall.dist <- lapply(rollcall.simple, function(m) dist(m %*% t(m)))

MDSを計算。

# 図にした時の配置を直観に合わせるため、-1を掛けて反転させてます。
> rollcall.mds <- lapply(rollcall.dist,
  function(d) as.data.frame((cmdscale(d, k = 2)) * -1))
> head(rollcall.mds[[1]])
          V1       V2
2  -11.44068 293.0001
3  283.82580 132.4369
4  885.85564 430.3451
5 1714.21327 185.5262
6 -843.58421 220.1038
7 1594.50998 225.8166

さらに、 rollcall.data から、議員名と所属政党を取り出して追加します。

> congresses <- 101:111
> for(i in 1:length(rollcall.mds)) {
  names(rollcall.mds[[i]]) <- c("x", "y")
  
  congress <- subset(rollcall.data[[i]], state < 99) # 副大統領は除外
  # 議員名を取り出す。姓名が含まれている場合があるので、名だけを取り出す。
  congress.names <- sapply(as.character(congress$name),
    function(n) strsplit(n, "[, ]")[[1]][1])
  
  rollcall.mds[[i]] <- transform(rollcall.mds[[i]],
                                 name = congress.names,
                                 party = as.factor(congress$party),
                                 congress = congresses[i])
}
> head(rollcall.mds[[1]])
           x        y      name party congress
2  -11.44068 293.0001    SHELBY   100      101
3  283.82580 132.4369    HEFLIN   100      101
4  885.85564 430.3451   STEVENS   200      101
5 1714.21327 185.5262 MURKOWSKI   200      101
6 -843.58421 220.1038 DECONCINI   100      101
7 1594.50998 225.8166    MCCAIN   200      101

データができたので、グラフにしてみます。 まずは、第110上院議会の散布図。

> cong.110 <- rollcall.mds[[9]]
> base.110 <- ggplot(cong.110, aes(x = x, y = y)) +
  scale_size(range = c(2,2), guide = 'none') +
  scale_alpha(guide = 'none') +
  theme_bw() +
  theme(axis.ticks = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        panel.grid.major = element_blank()) +
  ggtitle("Roll Call Vote MDS Clustering for 110th U.S. Senate") +
  xlab("") +
  ylab("") +
  scale_shape(name = "Party", breaks = c("100", "200", "328"),
              labels = c("Dem.", "Rep.", "Ind."), solid = FALSE) +
  scale_color_manual(name = "Party", values = c("100" = "black",
                                                "200" = "dimgray",
                                                "328"="grey"),
                     breaks = c("100", "200", "328"),
                     labels = c("Dem.", "Rep.", "Ind.")) +
  geom_point(aes(shape = party, alpha = 0.75, size = 2))
> ggsave(plot = base.110, filename = 'plot2.png')

所属政党でざっくり二分される結果に。まぁ、そうですよね。

最後に、残りの議会の図も作成して終わり。

> all.mds <- do.call(rbind, rollcall.mds)
> all.plot <- ggplot(all.mds, aes(x = x, y = y)) +
  geom_point(aes(shape = party, alpha = 0.75, size = 2)) +
  scale_size(range = c(2, 2), guide = 'none') +
  scale_alpha(guide = 'none') +
  theme_bw() +
  theme(axis.ticks = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        panel.grid.major = element_blank()) +
  ggtitle("Roll Call Vote MDS Clustering for U.S. Senate (101st - 111th Congress)") +
       xlab("") +
       ylab("") +
       scale_shape(name = "Party",
                   breaks = c("100", "200", "328"),
                   labels = c("Dem.", "Rep.", "Ind."),
                   solid = FALSE) +
      facet_wrap(~ congress)
> ggsave(plot = all.plot, filename = 'plot3.png')