機械学習手習い: クラスタリング
「入門 機械学習」手習い、9日目。「9章 MDS:米国上院議員の類似度の視覚的な調査」です。
観測値をクラスタリングするための、多次元尺度構成法(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')
