RでAffinity propagation - fukuitのブログ

とりあえず、必要があってRでaffinity propagationを書いてみた。MATLABやCのコードは公開されているけれど、Rのコードは公開されてない。必要なら作るしかない。

Rでちゃんと関数を書く方法がよく分からないので、その辺は適当に。

まずは、Rでdata.frame間の距離を測ってmatrixにする - fukuitの日記の続き。入力された、データポイント同士の距離を求めて、affinity propagationで言うところの、similarityの行列を作る。この距離はユークリッド距離に限定しちゃってるけど、本来は引数で指定できるようにすべきなんだろうな。そういうのの書き方が分からない。

ちなみに、データポイント同士のユークリッド距離を求めて、ソレを二乗した物に-1をかけている。さらに、そうすると、自身に対する距離は0になってしまって、後々過剰に評価されてしまうので、自身に対するsimilarityは全てのデータポイント同士の距離のmedianを設定する。なんで、medianなんだろう。

sim <- function( df ){      #### make similarity table
  if( !is.data.frame(df) ){
    df <- as.data.frame(df)
  }
  len <- nrow(df)
  s <- as.matrix(dist(df, method="euclidean"))
  c <- c()  # datapoints (except i==j)
  for(i in 1:len)
    for(j in 1 :len){
      # convert to negative euclidean distance
      s[i,j] <- -1*s[i,j]^2 
      if ( i != j )
        c<-c(c,s[i,j])
    }
  # preference <- median value ... preference means self similarity score
  m <- median(c) 
  for(i in 1:len)
    s[i,i] <- m     
  return(s)
}

次は、responsibilityを計算する。

responsibility <- function(s,  # similarity table
                           a,  # availability table
                           r   # responsibility table
                           ){
  len <- nrow(s)  
  for(i in 1:len){
    for(k in 1:len){
      tmp<-c()
      if (i != k){
        for(kk in 1:len)
          if(k != kk)
            tmp <- c(tmp, a[i,kk]+s[i,kk]) 
      }else{   # for r(k,k) ... self responsibility
        for(kk in 1:len)
          if(k != kk)
            tmp <- c(tmp, s[i,kk]) 
      }
      r[i,k] <- s[i,k] - max(tmp)
    }
  }
  return(r)
}

コレを計算するためには、availability tableが必要だが、availabilityの初期値は0なので、一発目の計算はソレを代入する。
ここで求めたresponsibilityの値を使って、availabilityを計算する。

availavility<-function(  s,    # similarity table
                         r,    # responsibility table
                         a     # availablity table           
                       ){
  len <- nrow(s)
  for(i in 1:len){
    for(k in 1:len){
      tmp <- 0
      for(ii in 1:len)
        if ((ii != k) && (ii != i) && (r[ii,k] > 0))
          tmp <- tmp + r[ii,k]
      if (i != k)
        a[i,k] <- min(0,r[k,k]+tmp) # for self availability
      else
        a[i,k] <- tmp
    }
  }
  return(a)
}

多分、もっと良い方法があるんだろうな、と、思いながら書いてるワケだけれど、そういうことを相談できる相手が身近に欲しいところだ。
responsibilityとavailabilityを互いに互いを使って、値が収束するまで計算を繰り返す。このiterationが収束するのは、引数として設定する上限回数に達した時か、あるいはavailabilityとresponsibilityを加算して求めたcluster centerが一定回数以上変化しなかった時だ。

apclust <- function(df,
                lambda=0.5,
                maxits=100,
                convits=10)
{
  len<-nrow(df)             # count of datapoint
  a <- diag(0, len)         # availability
  r <- diag(0, len)         # responsibility
  e <- diag(0, len)         # exemplar
  kc <- c()                 # clustercenter
  kcPrev <- c()             # clustercenter of previous iteration
  kcnt <- 0                 # iteration count

  s<-sim(df)               # similarity table  
  # remove degeneracies
  s <- s+(1e-12)*randn(len)*(max(s)-min(s))
  # iteration
  for (iter in 1:maxits){
    cat("*")
    aPrev<-a
    rPrev<-r
    kcPrev<-kc
 
    r <- responsibility(s,a,r)       ## get responsibility
    r <- (1-lambda)*r+lambda*rPrev   ## dumping factor lambda
    a <- availavility(s,r,a)         ## get availability
    a <- (1-lambda)*a+lambda*aPrev   ## dumping factor lambda
    e <- r+a                         ## calc exemplar
    
    kc <- c()  # kc for cluster center
    for(i in 1:len)
      if (e[i,i] > 0)
        kc <- c(kc,i)  
     
    if( (!is.null(kc)) && (!is.null(kcPrev)) ){
      if ( (length(kc) == length(kcPrev)) &&
          (dist(rbind(kc, kcPrev), method="euclidean") == 0) ){  
        kcnt <- kcnt+1
        if (kcnt > convits)
          break
      }
    }else{
      kcnt <- 0
    }
  }
 
  ### print output
  print(e)    # exemplar table
  print(iter) # iteration count
  cl <- c()
  for(i in 1:len){
    if ( max(e[i,]) > 0 )
      cl <- c(cl, which.max(e[i,])) # cluster
    else
      cl <- c(cl, -1*which.max(e[i,])) # not cluster center
  }
  samplename <- rownames(df)
  result <- data.frame(SAMPLE=samplename, EXEMPLAR=cl)
  return(list(exemplar=e, center=kc, cluster=result))
}

これで、代表値が求められるはず。とりあえず、Affinity Propagation Web Applicationに載ってるToy Problemと同じ結果は得られたので、効率等はさておき、それなりに動作しているんだろうとは思う。

動作確認

40個の乳癌細胞と7個の正常細胞の発現データを使って、試してみる。つまり、答えが既に分かっているデータを使って、正しく乳癌細胞と正常細胞を見分けられるかどうかを確認する。
もう雑に、esetを取ってきて、特に正規化もしないまま発現量のsimilarityから、分類する。

 > library(GEOquery)
 > gse3744.eset <- getGEO('GSE3744')[[1]]
 > apclust(t(exprs(gse3744.eset)))

すると、とりあえず乳癌細胞はいくつかのクラスタに分かれたが、正常細胞に関しては、一つの正常細胞のみからなるクラスタが形成された。乳癌細胞もそれぞれのクラスタごとの特徴みたいなものをつかめると、これはこれで面白い研究ネタに繋がるかもしれない(繋がらないかもしれない)。

という具合に、オレが書いた程度のRの関数でも、とりあえず一瞬でクラスタリングできるという意味で、affinity propagationはすごい技術なのかもよ。