授業用に作ったものを、こちらにも載せておく。

ちょっとした説明

今井. 2018. 『社会科学のためのデータ分析入門(上)』第3章に記載されている、DW-NOMINATEの1次元目及び2次元目のスコアを元に連邦下院議員を分類するタスクを扱う。分析の対象は第80議会 (1947-48)、クラスターの数は $k = 4$ に設定する。

ここでは、アルゴリズムの挙動がよくわかるように、Rの標準関数であるkmeans() は用いず、各ステップが終了するごとにその時点での分類結果を図示するようにしている。

Rコード

# Load package
require(tidyverse) # for alpha() function

# Load & subset data
url <- "https://raw.githubusercontent.com/kosukeimai/qss/master/MEASUREMENT/congress.csv"
congress <- read.csv(url)
congress.80 <- subset(congress, congress == 80)

# k-means settings
k <- 4 # number of clusters
max.iter <- 100 # maximum number of iterations
set.seed(1234)

# Initialize 
congress.80$class <- sample(1:4, nrow(congress.80), replace = TRUE)
congress.80$class.new <- NA
plot(dwnom2 ~ dwnom1, data = congress.80, pch = 16, col = alpha(class + 1, 0.2),
     xlab = "DW-NOMINATE Score (1st dim)", ylab = "DW-NOMINATE Score (2nd dim)",
     main = "Iteration: 0")
x.mean <- tapply(congress.80$dwnom1, congress.80$class, mean) # compute centroids
y.mean <- tapply(congress.80$dwnom2, congress.80$class, mean)
points(x.mean, y.mean, pch = 15, col = 2:5, cex = 1.5)
for (i in 1:max.iter){
  # Reassign each point to the nearest cluster
  for (j in 1:nrow(congress.80)){ # loop over legislators
    dist <- sqrt((congress.80$dwnom1[j] - x.mean)^2 + (congress.80$dwnom2[j] - y.mean)^2)
    congress.80$class.new[j] <- which(dist == min(dist))
  }
  # Plot the current class assignments
  plot(dwnom2 ~ dwnom1, data = congress.80, pch = 16, col = alpha(class.new + 1, 0.2),
       xlab = "DW-NOMINATE Score (1st dim)", ylab = "DW-NOMINATE Score (2nd dim)",
       main = paste0("Iteration: ", i))
  x.mean <- tapply(congress.80$dwnom1, congress.80$class.new, mean) # recompute centroids
  y.mean <- tapply(congress.80$dwnom2, congress.80$class.new, mean)
  points(x.mean, y.mean, pch = 15, col = 2:5, cex = 1.5)
  # break if class assignment doesn't change
  if (sum(congress.80$class != congress.80$class.new)== 0) break 
  # update class assignment
  congress.80$class <- congress.80$class.new 
  Sys.sleep(2)
}