お正月に広津先生のクラスタリング手法をRで実装しました。個体ごとの時系列データをクラスタリングするのに使えます(実際は時系列ではなく一般の2-wayのデータに適用できます)。
個人的に感じる適正なサンプルサイズと時点のサイズはおよそ、10~1000個体、4~30時点程度です。これ以上時点が多い場合は、状態空間モデルなどの方がよいと思われます。
参考文献として以下の3つを挙げます。
[1] Hirotsu, C.(2009): Clustering rows and/or columns of a two-way contingency table and a related distribution theory. Computational Statistics and Data Analysis 53, 4508-4515
[2] 広津千尋 (2004) : 交互作用は相互作用?2 コレステロール低下剤Mは有効か?, ESTRELA2004年3月号
[3] 広津千尋 (1992) : 実験データの解析-分散分析を超えて-, 共立出版, 東京
[1]は今回のアルゴリズムの原論文です。 [2]は日本語で分かりやすいです(ただし(3)式の係数行列、(6)式の符号などにtypoがあるので注意)。 [3]の10章には論文の旧バージョンのアルゴリズムが載っており、また背景や計算などが詳しく書かれています(ただし図10.1及び図10.2が欠落しているので注意)。
プログラム内で使用しているデータは[2]と[3]で使われているもので、ほんとの薬を飲ませた人が12人、プラセボ(なんちゃって薬を飲ませた人)の人が11人の合計23人の体内総コレステロール量の時間変化(6時点)のデータになります。 ほんとの薬がどれほど有効か示すにはいくつか方法があると思いますが、ここでは、(1)薬とプラセボを混ぜてクラスタリング→(2)薬とプラセボを行、どのグループに属したかを列にした混合表を作ってFisherの正確確率検定で判定 という流れで評価しています。 その他、グループごとの交互作用を図示したい時などは[2]を参照してください。
さて、Rのソースコードは以下になります。
hirotsu.clustering <- function(data.m, cluster.n){ t.n <- ncol(data.m) s.n <- nrow(data.m) cache.m <- cbind(data.m[,1], sapply(2:t.n, function(j){ apply(data.m[,1:j]/j, 1, sum) })) - apply(data.m, 1, mean) get.diff.group <- function(gr1, gr2){ n1 <- length(gr1) n2 <- length(gr2) s <- sum(sapply(1:(t.n - 1), function(j) { j/(t.n - j) * (sum(cache.m[gr1, j])/n1 - sum(cache.m[gr2, j])/n2)**2 })) n1*n2/(n1+n2) * t.n * s } get.diff.groups <- function(gr.l){ sum(sapply(1:(length(gr.l)-1), function(ix){ get.diff.group(unlist(gr.l[1:ix]), gr.l[[ix+1]]) })) } get.sigstat <- function(gr.l){ j.v <- 1:(t.n-1) rho.v <- sqrt(t.n/j.v/(t.n-j.v)) chi.square.inv <- sum(((t(data.m[,-t.n] - data.m[,-1]) - (colMeans(data.m)[-t.n] - colMeans(data.m)[-1]))/rho.v)**2) get.diff.groups(gr.l)/chi.square.inv } merge.n <- s.n - cluster.n member.l <- as.list(1:s.n) alive.v <- 1:s.n dist.m <- matrix(NA, s.n + merge.n, s.n + merge.n) for (i1 in 1:(s.n-1)){ for (i2 in (i1+1):s.n){ dist.m[i1, i2] <- get.diff.group(member.l[[i1]], member.l[[i2]]) } } for (m in 1:merge.n){ new.i <- s.n + m min.ix <- which(dist.m == min(dist.m, na.rm=T), arr.ind=T) new.member.v <- sort(as.vector(min.ix[1,])) member.l[[new.i]] <- sort(unlist(sapply(new.member.v, function(x){ member.l[[x]] }))) alive.v <- alive.v[-which(alive.v %in% new.member.v)] for (i in alive.v){ dist.m[i, new.i] <- get.diff.group(member.l[[i]], member.l[[new.i]]) } alive.v <- c(alive.v, new.i) dist.m[new.member.v, ] <- NA dist.m[ ,new.member.v] <- NA } gr.l <- sapply(alive.v, function(x){ member.l[[x]] }) repeat{ try.l <- replicate(cluster.n, vector()) for (i in 1:s.n){ exclude.l <- lapply(gr.l, function(x){ x[which(x!=i)] } ) check <- sapply(exclude.l, length) if (any(check == 0)){ try.l[[which(check == 0)]] <- i next } nearest <- which.min(sapply(gr.l, function(x){ get.diff.group(i, x) })) try.l[[nearest]] <- c(try.l[[nearest]], i) } if (get.diff.groups(gr.l) > get.diff.groups(try.l) || identical(all.equal(gr.l, try.l), TRUE)) break gr.l <- try.l } list(grouping=gr.l, sigstat=get.sigstat(gr.l)) } cluster.n <- 3 data.m <- rbind(c(317, 280, 275, 270, 274, 266), c(186, 189, 190, 135, 197, 205), c(377, 395, 368, 334, 338, 334), c(229, 258, 282, 272, 264, 265), c(276, 310, 306, 309, 300, 264), c(272, 250, 250, 255, 228, 250), c(219, 210, 236, 239, 242, 221), c(260, 245, 264, 268, 317, 314), c(284, 256, 241, 242, 243, 241), c(365, 304, 294, 287, 311, 302), c(298, 321, 341, 342, 357, 335), c(274, 245, 262, 263, 235, 246), c(232, 205, 244, 197, 218, 233), c(367, 354, 358, 333, 338, 355), c(253, 256, 247, 228, 237, 235), c(230, 218, 245, 215, 230, 207), c(190, 188, 212, 201, 169, 179), c(290, 263, 291, 312, 299, 279), c(337, 337, 383, 318, 361, 341), c(283, 279, 277, 264, 269, 271), c(325, 257, 288, 326, 293, 275), c(266, 258, 253, 284, 245, 263), c(338, 343, 307, 274, 262, 309)) result <- hirotsu.clustering(data.m, cluster.n)
> result $grouping $grouping[[1]] [1] 1 3 9 10 23 $grouping[[2]] [1] 2 5 6 12 13 14 15 16 17 19 20 21 22 $grouping[[3]] [1] 4 7 8 11 18 $sigstat [1] 0.5269027 > p.value [1] 0.0029
hirotsu.clustering
関数が本体になります。NA
のないmatrix
と分けたいクラスター数を引数に取ります。result$grouping
でどの個体がどのグループに属しているか分かります。result$sigstat
は少し正確さを犠牲にして分かりやすく言うとグループ間距離を補正したような統計量です。
クラスタリングの結果の図は以下になります。
group1は減少群、group2は不変群、group3は増加群と判断できます。
この後に、ほんとの薬とプラセボの方々がどのグループに属したかの人数の分布が以下になります。
Treat | group 1 (減少群) | group 2 (不変群) | group 3 (増加群) |
---|---|---|---|
ほんとの薬 | 4 | 4 | 4 |
プラセボ | 1 | 9 | 1 |
Fisherの正確確率検定(両側検定)で調べてみるとp-value=0.090で有意ではない(ほんとの薬は有効とは言えない)という結果になります。