TwitterAPI v1.1 に対応

　Ruby+TwitterAPIを使用して日本やインドネシアのトレンド(ホットワード)を取得して関係者に配信していたが、6/12の朝、使用していたAPIのバージョンv1のサービスがついに終了してしまった。前から知っていたことではあるのだが...　使っていたサービスはトレンド取得だけ、OAuth認証も既に組み込んでいたので、エンドポイントであるURLだけ変更すれば簡単に終わるかとおもったら... やっぱり落とし穴？っていうか自分のテク不足だけなんだろう。google先生に聞いてみると、DevelopersサイトにSingle-user OAuth with Examplesページがあった。Rubyは一番下に書いています。今までいろいろなライブラリ駆使してコードを書いている感があったけど、oauthライブラリだけでいけそうだ。あとはJSONをパースする為のjsonライブラリくらいか...　ということでさっそく試す。
URLに引数でID(WOEID)を指定することにより対象地域のホットワードを得ることができる。あとはGETしたJSONのtrendsをぐるぐる回しながら単語を出力。環境はWin7 32bit、Ruby1.8.7

require 'rubygems'
require "oauth"
require 'json'

WOEID = ARGV[0] # 引数でホットワード取得場所のIDを得る
#WOEID = '23424856' #日本
#WOEID = '1117034' #千葉

# Exchange your oauth_token and oauth_token_secret for an AccessToken instance.
def prepare_access_token(oauth_token, oauth_token_secret)
  consumer = OAuth::Consumer.new("ここはConsumer Key", "ここはConsumer secret",
    { :site => "https://api.twitter.com",
      :scheme => :header
    })
  # now create the access token object from passed values
  token_hash = { :oauth_token => oauth_token,
                 :oauth_token_secret => oauth_token_secret
               }
  access_token = OAuth::AccessToken.from_hash(consumer, token_hash )
  return access_token
end
 
# Exchange our oauth_token and oauth_token secret for the AccessToken instance.
access_token = prepare_access_token("ここはAccess token", "ここはAccess token secret")
response = access_token.request(:get, "https://api.twitter.com/1.1/trends/place.json?id=" + WOEID + "&exclude=hashtags") # id=WOEID ハッシュタグを除く
t = Time.now # 時間取得
uploadDate = t.strftime("%Y/%m/%d %H:%M:%S.#{t.usec.to_s[0, 3]}") # 時間フォーマット

if response.code == "200" then
  hash = JSON.parse(response.body) # JSON を parse
  seq = 0
  hash[0]["trends"].each do |trends| # treds のループ
    seq += 1
    name = trends["name"].to_s.kconv(Kconv::SJIS, Kconv::UTF8) # ホットワード
    url  = trends["url"].to_s.kconv(Kconv::SJIS, Kconv::UTF8) # URL
    rowdata = [uploadDate, WOEID, seq, name, url].join("\t") # タブ区切り結合
    print(rowdata,"\r\n") # 出力
end

無事成功！　結果の中に”父の日”があった。そうか明日は父の日だった。どうりで子供たちがこそこそしているわけだ。

Windows環境でjsonライブラリをインストールする場合

　Windows環境にjsonライブラリを普通にインストールするとエラーになる。

>gem install json
ERROR:  http://rubygems.org/ does not appear to be a repository
ERROR:  Could not find a valid gem 'json' (>= 0) in any repository

こちらのサイトにあるjson-1.1.9-x86-mswin32.gemをダウンロードしてlocalファイルをインストールする必要がある。1年もたつともう忘れている。

>gem install json --local
Successfully installed json-1.1.9-x86-mswin32
1 gem installed
Installing ri documentation for json-1.1.9-x86-mswin32...
Installing RDoc documentation for json-1.1.9-x86-mswin32...

これでOK!

対応している地域が増えていた

　こんな便利なサイトを使うと開発が楽になる。twitSandbox1.1　認証できていれば確認したいAPIの結果をこのサイトで見ることができる。trends/availableでは現在ホットワードが提供されている地域を確認することが出来る。見ると千葉があった。前はなかったはず...　バージョンアップで増えたのかな...　また、trends/closestでは、緯度と経度を指定すると最寄のWOEIDを知ることもできる。住んでいる船橋市の緯度経度で試すと結果は千葉の1117034だった。あとは会社のproxy越えのコードを追加するだけか。とりあえず週末のお仕事1つ終了。

マルコフ連鎖とPageRank （追加メモ）

　先日【Rメモ】マルコフ連鎖とPageRankといったメモを書いた。マルコフ連鎖の収束値とPageRankの値が一致するといった点についてふれたが、一致するのはdanmping=1の場合。danmping factorの違いによるPageRankの値を確認したのでちょっとだけ追加メモ。

　wikiのページランンクを見ると、作為的にページランクを上げようとする者に対しては、より小さい値に設定される。といった記述があるが、定義式を見るとdanmping factorが何をしているかすぐに分かる。定義式はwiki(English)が解りやすい。

d=0なら、1/N になるので均等な値になるってことです。Rのigraphパッケージのpage.rank()関数で確認。データは前回と同じこんな遷移確率の行列を準備。

	V1	V2	V3	V4	V5	V6
V1	0.700	0.360	0.070	0.040	0.100	0.040
V2	0.100	0.240	0.090	0.080	0.050	0.010
V3	0.080	0.300	0.700	0.400	0.180	0.070
V4	0.020	0.050	0.050	0.250	0.170	0.060
V5	0.030	0.030	0.020	0.130	0.270	0.130
V6	0.070	0.020	0.070	0.100	0.230	0.690
	1.000	1.000	1.000	1.000	1.000	1.000

※実際にRにくわせるデータには行と列のヘッダー、最終行の1.000は不要です。

library(igraph)
library(reshape)
library(ggplot2)

d1 <- as.matrix(read.table("遷移確率.txt", header=F)) # 遷移確率をmatrix形式で読込む

# 遷移確率行列からグラフ形式データを準備
d3 <- d1 # 遷移確率をコピー
rownames(d3) <- colnames(d3) # rownamesの設定
d3 <- as.data.frame(t(d3)) # from to 形式にする為の転置
d3 <- cbind(from=rownames(d3), d3) # 行名を列として設ける
d3 <- melt(d3, id="from") # from to weight 形式に整形
head(d3)

> head(d3)
  from variable value
1   V1       V1  0.70
2   V2       V1  0.36
3   V3       V1  0.07
4   V4       V1  0.04
5   V5       V1  0.10
6   V6       V1  0.04

g <- graph.data.frame(d3[1:2], directed=T) # graph形式に変換
E(g)$weight <- d3[[3]] # weightは遷移確率

df <- c(0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0) # 計算するdamping factor
pr.all <- c()
for(i in 1:length(df)){
  pr.all <- rbind(pr.all, page.rank(g, directed=T, damping=df[i])$vector) # PageRank計算
}

# ここからは描画
pr.all <- melt(data.frame(df=df, pr.all), id="df") # ggplot2用のデータ形式に変換
gg2 <- ggplot(pr.all, aes(x=reorder(df, df), y=value, fill=variable)) # 積上げ棒グラフ reorderで横軸順序を指定
gg2 <- gg2 + geom_bar(colour=I("black"), alpha=0.8, linetype=1, size=0.7) # ちょっとキレイにする
plot(gg2) # 描画

　棒グラフ右はdamping=1なのでマルコフ連鎖の収束値、左の棒グラフはdamping=0なので各系列が均等になっている様子、そしてdampingの違いによるPageRankの変化を確認できました。

参考
・Googleを支える技術 ?巨大システムの内側の世界 (WEB+DB PRESSプラスシリーズ)
・Google の秘密 - PageRank 徹底解説

【Rメモ】マルコフ連鎖とPageRank

　マーケティングではSTP戦略という手法がよく使われます。マーケティングの第一人者フィリップ・コトラー氏が提唱した代表的な手法。効果的に市場を開拓する為に 1.市場を顧客のニーズごとにセグメント化　2.参入すべきセグメントを選定しターゲットを明確にする 3.セグメントの中で顧客に対するベネフィット（利益）を検討し自らのポジションを確立する といった流れ。この戦略のS(セグメント)を議論する時、私はたたき台としてマルコフ連鎖をよく使っています。シンプルな手法ですが変化を可視化して共有すると議論がもりあがります。マーケティング担当者の受けも上々。今まではExcelで分析していたがRで完結させたくなりコードを書いたのでφ(..)メモメモ。

Rでマルコフ連鎖

　R+マルコフ連鎖については、teramonagiさんの資料 がよく分かります。（Rコードも記載されてます）。先ほどのマーケティングの例なら、現在のセグメント構成比が初期状態で、そのセグメント間の遷移確率が分かれば次の状態（セグメント構成比）を計算できるといったものです。計算は隣接行列×状態の構成比といった行列の演算、ExcelだとMMULT関数で一発。当然Rも行列の掛け算は %*% で行うことが出来るので一発。あとはその演算の繰り返し。自分で書くつもりでしたが、先ほどのteramonagiさんの資料にその関数が記載されていたのでそれを拝借。以下実際のデータとRコード。

初期構成比のデータ（セグメントの構成比）

V1	0.080
V2	0.060
V3	0.300
V4	0.100
V5	0.100
V6	0.360

遷移確率のデータ（各セグメント間で遷移する比率 ※V2->V3の確率は0.3）

	V1	V2	V3	V4	V5	V6
V1	0.700	0.360	0.070	0.040	0.100	0.040
V2	0.100	0.240	0.090	0.080	0.050	0.010
V3	0.080	0.300	0.700	0.400	0.180	0.070
V4	0.020	0.050	0.050	0.250	0.170	0.060
V5	0.030	0.030	0.020	0.130	0.270	0.130
V6	0.070	0.020	0.070	0.100	0.230	0.690
	1.000	1.000	1.000	1.000	1.000	1.000

setwd("d:/R/マルコフ連鎖") # 作業ディレクトリの変更

d1 <- as.matrix(read.table("遷移確率.txt", header=F)) # 遷移確率をmatrix形式で読込む
d2 <- as.matrix(read.table("初期構成比.txt", header=F)) # 初期構成比をmatrix形式で読込む

#----------------------------------------------------------
# 参考：http://www.slideshare.net/teramonagi/ss-5190440 p17
#
# 数日後の天気確率計算関数
# weather.initial_:
# transition.matrix_:お天気遷移行列
# size_:何日間の遷移を計算するのか
TranslateWeather <- function(weather.initial_, transition.matrix_, size_){
  result <- matrix(NA, nrow = size_ + 1, ncol = nrow(weather.initial_)) # 格納先
  result[1,] <- t(weather.initial_) # 初期値を結果に格納
  weather.day <- weather.initial_ # 初期値
  for(day in 1:size_){
    weather.day <- transition.matrix_ %*% weather.day # マルコフ連鎖(行列演算)
    result[day + 1,] <- t(weather.day) # 結果を転置して格納
  }
  return(result)
}
#----------------------------------------------------------
res <- TranslateWeather(d2, d1, 100) # 引数：初期値、遷移行列、回数
head(res)

          [,1]       [,2]      [,3]       [,4]       [,5]      [,6]
[1,] 0.0800000 0.06000000 0.3000000 0.10000000 0.10000000 0.3600000
[2,] 0.1270000 0.06600000 0.3176000 0.08320000 0.09700000 0.3092000
[3,] 0.1602880 0.07172200 0.3246640 0.07756200 0.08934400 0.2764200
[4,] 0.1838417 0.07589820 0.3280606 0.07418924 0.08359412 0.2544162
[5,] 0.2004804 0.07878419 0.3296510 0.07189805 0.07964253 0.2395439
[6,] 0.2121961 0.08075424 0.3302923 0.07031774 0.07696189 0.2294777

ggplot2で可視化

　先ほどのマルコフ連鎖の結果を ggplot2パッケージを使い可視化。(ggplot2で書きたくなったので)

library(ggplot2)
library(reshape)

res <- as.data.frame(TranslateWeather(d2, d1, 10)) # 引数：初期値、遷移行列、回数
res2 <- cbind(t=as.numeric(rownames(res)), res) # 行番号を列として設ける
res3 <- melt(res2, id="t") # ggplot2用のデータ形式に変換
head(res2)

  t        V1         V2        V3         V4         V5        V6
1 1 0.0800000 0.06000000 0.3000000 0.10000000 0.10000000 0.3600000
2 2 0.1270000 0.06600000 0.3176000 0.08320000 0.09700000 0.3092000
3 3 0.1602880 0.07172200 0.3246640 0.07756200 0.08934400 0.2764200
4 4 0.1838417 0.07589820 0.3280606 0.07418924 0.08359412 0.2544162
5 5 0.2004804 0.07878419 0.3296510 0.07189805 0.07964253 0.2395439
6 6 0.2121961 0.08075424 0.3302923 0.07031774 0.07696189 0.2294777

head(res3)

  t variable     value
1 1       V1 0.0800000
2 2       V1 0.1270000
3 3       V1 0.1602880
4 4       V1 0.1838417
5 5       V1 0.2004804
6 6       V1 0.2121961

gg2 <- ggplot(res3, aes(x=reorder(t, t), y=value, fill=variable)) # 積上げ棒グラフ reorderで横軸順序を指定
gg2 <- gg2 + geom_bar(colour=I("black"), alpha=0.8, linetype=1, size=0.7) # ちょっとキレイにする
plot(gg2)

結構きれい。この遷移確率が続くといった仮定は入りますが、初期状態では小さいV1と大きなV6は逆転するといった結果を知ることが出来ます。このような積上げ棒グラフを描きたい場合、x軸の順序を並べ替える必要がある場合がありますのでご注意ください。実際にハマってしまいました。

マルコフ連鎖 と PageRank

　マルコフ連鎖は時間tを無限大にすると収束します。収束値は連立方程式で解けます。（これは会社の先輩から教えてもらいました　※平衡方程式というらしい）　ご存じの方もいると思いますが、その収束値はPageRankと同じ値になります。私はたまたま知ったが。。。（1年くらい前にネットワーク図のノードの評価について駒場の学生さんとお話しする機会がありPageRankが話題に。職場に戻り算出していたマルコフ連鎖の収束値とPageRankの相関をとったら一致して気づく）　RではigraphパッケージでPageRankを計算できるので一致するか確認。

library(igraph)
library(reshape)

d3 <- d1 # 遷移確率をコピー
rownames(d3) <- colnames(d3) # rownamesの設定
d3 <- as.data.frame(t(d3)) # from to 形式にする為の転置
d3 <- cbind(from=rownames(d3), d3) # 行名を列として設ける
d3 <- melt(d3, id="from") # from to weight 形式に整形
head(d3)

  from variable value
1   V1       V1  0.70
2   V2       V1  0.36
3   V3       V1  0.07
4   V4       V1  0.04
5   V5       V1  0.10
6   V6       V1  0.04

g <- graph.data.frame(d3[1:2], directed=T) # graph形式に変換
E(g)$weight <- d3[[3]] # weightは遷移確率
pr <- page.rank(g, directed=T, damping=1.0) # dampingのデフォルトは0.85
round(pr$vector, 5) # PageRankの結果

     V1      V2      V3      V4      V5      V6 
0.23924 0.08494 0.32928 0.06684 0.07135 0.20835

res <- TranslateWeather(d2, d1, 100) # 100回計算させる
tail(round(res, 5)) # マルコフ連鎖の結果

          [,1]    [,2]    [,3]    [,4]    [,5]    [,6]
[96,]  0.23924 0.08494 0.32928 0.06684 0.07135 0.20835
[97,]  0.23924 0.08494 0.32928 0.06684 0.07135 0.20835
[98,]  0.23924 0.08494 0.32928 0.06684 0.07135 0.20835
[99,]  0.23924 0.08494 0.32928 0.06684 0.07135 0.20835
[100,] 0.23924 0.08494 0.32928 0.06684 0.07135 0.20835
[101,] 0.23924 0.08494 0.32928 0.06684 0.07135 0.20835

PageRankの結果とマルコフ連鎖の収束値が一致していることが確認できます。（100回計算しなくても既に収束しているようですが）　また、PageRankにはdampingパラメータがあります。デフォルトは0.85ですが1にしないとマルコフ連鎖の収束値にはなりません。

参考
・マルコフ連鎖モンテカルロ法入門－1
・reorder 関数を使って横軸の順序の並べ替えを行う方法

R changepoint パッケージちょっと試した

　よく見ている R-bloggers で R で changepoint パッケージを使ったコードが掲載されていた。（実際のサイトはこちら）

変化点検出は、これまでもいくつか書籍やブログ、R勉強会のTokyoRなどで紹介されている。

・データマイニングによる異常検知　ChangeFinderが掲載されている山西さんの書籍
・異常検知（変化点検出）のパッケージを作ってみた TokyoR主催者yokkunsさんの日記
・隠れマルコフモデルで異常検知（R Advent Calendar2012） teramonagiさんのブログ（コードはこちら）

changepoint パッケージは初めてだったので、ちょっとだけ試してみた。まずは inside-R にあったサンプルコードを実行。サンプルデータは正規分布の乱数です。

library(changepoint)

# change in variance
x=c(rnorm(100,0,1),rnorm(100,0,10))
ansvar=cpt.var(x)
plot(ansvar)
print(ansvar)
 
# change in mean
y=c(rnorm(100,0,1),rnorm(100,5,1))
ansmean=cpt.mean(y)
plot(ansmean,cpt.col='blue')
print(ansmean)
 
# change in mean and variance
z=c(rnorm(100,0,1),rnorm(100,2,10))
ansmeanvar=cpt.meanvar(z)
plot(ansmeanvar,cpt.width=3)
print(ansmeanvar)

マニュアルにも cpt.mean、cpt.var、cpt.meanvar を見てから始めて．．．といった説明もあるので、基本的な関数はこの3つ？　それぞれの plot はこのようになります。

change in variance

change in mean

change in mean and variance

といった感じです。具体的なデータではどのようになるの？ということで日経平均株価を試してみましたのでその結果を紹介します。データはダウンロードサイトから csv ファイルを準備しました。（R では yokkunsさんが開発した株価が取得できる RFinanceYJ というパッケージもあります）csv ファイルは5列目が終値です。また、関数は cpt.meanvar、引数の method は AMOC、PELT、SegNeigh、BinSeg の中から PELT を選択しております。

library(changepoint)

## ---日経平均日足データの読込み
setwd("d:/R/changepoint") # 作業ディレクトリの変更

d <- read.csv("日経平均日足.csv", header=T) # データフレームで読込み

head(d)
       日付    始値    高値    安値    終値
1  2012/5/8 9190.50 9207.56 9159.47 9181.65
2  2012/5/9 9112.72 9115.94 9021.20 9045.06
3 2012/5/10 9013.26 9075.63 8985.90 9009.65
4 2012/5/11 9019.40 9050.61 8944.63 8953.31
5 2012/5/14 8986.22 9031.09 8947.82 8973.84
6 2012/5/15 8910.85 8930.78 8838.78 8900.74


ansmeanvar=cpt.meanvar(d[,5], method='PELT')
plot(ansmeanvar,cpt.width=3)
print(ansmeanvar)

Class 'cpt' : Changepoint Object
       ~~   : S4 class containing 10 slots with names
              date data.set cpttype method test.stat pen.type pen.value cpts ncpts.max param.est 

Created on  : Sun May 12 17:57:17 2013 

summary(.)  :
----------
Changepoint type      : Change in mean and variance 
Method of analysis    : PELT 
Test Statistic  : Normal 
Type of penalty       : SIC with value, 11.04292 
Maximum no. of cpts   : Inf 
Number of changepoints: 42 

一応変化点が検出されてっるぽい。まだまだ勉強不足なので時間があったらもう少し調べる。

参考
・changepoint: An R Package for Changepoint Analysis
・データマイニングによる異常検知

・異常検知（変化点検出）のパッケージを作ってみた
・隠れマルコフモデルで異常検知（R Advent Calendar2012）

ネットワーク分析をもうちょっと勉強（実践編）

　前回は 竹本先生の資料 や WEBページ を拝見し R で行うネットワーク分析について新たに知ることができた関数、分析手法についてまとめてみた。さっそく Non-overlapping 手法を実務で実践したが早くも壁にぶちあたってしまった。実務で取り扱うデータを分析する際にはつきものだ。。。今回分析したネットワークはこんな感じ。ノード数=446、エッジ数=792、図ではエッジの矢印が見えないかもしれないが有向グラフ。（前回同様サイトスケープにくわせOrganicレイアウトで表示　※背景色のみ白に変更）

さっそくエラー多発

　何も気にせずにNon-overlapping手法のコミュニティ抽出関数を実行したらエラーが多発。

> # おまじない
> ls(all=T)
character(0)
> rm(list=ls(all=TRUE))
> 
> # パッケージの読み込み
> library(igraph)
> 
> # ワーキングディレクトリの変更
> setwd("d:/R/igraph") # 作業ディレクトリの変更
> 
> # ネットワークの読込み
> d <- read.table("エッジデータ.txt", header=T) # データフレームで読込み
> g <- graph.data.frame(d[,2:3], directed = TRUE) # データフレームからグラフオブジェクトに変換
> g <- simplify(g, remove.multiple = TRUE, remove.loops = TRUE) # 多重エッジや自己ループの削除
> 
> eb <- edge.betweenness.community(g) # 辺媒介性 媒介性の高い辺を削除して分割
> rw <- walktrap.community(g) # ランダムウォークに基くアルゴリズム
> fg <- fastgreedy.community(g) # 貪欲アルゴリズム Q値が高くなるいようにペアをまとめていく
 以下にエラー fastgreedy.community(g) : 
  At fast_community.c:537 : fast greedy community detection works for undirected graphs only, Unimplemented function call
> le <- leading.eigenvector.community(g) # スペクトラル最適化法 グラフラプラシアンを使ってQ値が最大となるような分割を探す
 警告メッセージ： 
In leading.eigenvector.community(g) :
  At community.c:1357 :This method was developed for undirected graphs
> ml <- multilevel.community(g) # 多段階最適アルゴリズム
 以下にエラー multilevel.community(g) : 
  At community.c:2421 : multi-level community detection works for undirected graphs only, Unimplemented function call
> sp <- spinglass.community(g) # 焼きなまし法 グラフラプラシアンを使ってQ値が最大となるような分割を探す
 以下にエラー spinglass.community(g) : 
  At clustertool.cpp:285 : Cannot work with unconnected graph, Invalid value
> op <- optimal.community(g) # 全探索（時間がかかるため小さなネットワークのみ有効）
 以下にエラー optimal.community(g) : 
  At glpenv05.c:70 : glp_malloc: no memory available
, CPU time exceeded
>

うひょ～　なんだ？　辺媒介性とランダムウォークは問題なし、貪欲アルゴリズムと多段階最適アルゴリズムは無向グラフのみか。（スペクトラル最適化法は警告）全探索はマシンのスペック不足のようだ。最も試したかった焼きなまし法はネットワークが分断されているということか。上の図を見ても確認できるがネットワークは既にいくつかのコミュニティに分かれているってことだ。前に igraph を使っていたときは clusters 関数を使い連結成分のみのネットワークに分割し、その単位でコミュニティ抽出していたことを思い出した。前処理、前処理。

ネットワークの連結成分を確認する

　igraph の clusters 関数は連結成分を確認することができる。（Rで学ぶデータサイエンスの#8ネットワーク分析にも記載されている）

> g.cl <- clusters(g) #連結成分を確認
> g.cl$membership # 連結成分
  [1] 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 3 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 4 1 1 5 1
 [39] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 6 1 1 1 1 1 7 1 1 1 1 1 1 1 4 1 1 1 1

:(中略)

[419] 1 1 1 1 8 8 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
> g.cl$csize # 連結成分のサイズ
[1] 400   9   5  12   4   2  11   3
> g.cl$no # 連結成分の数
[1] 8

そうそう、この関数で分断されているネットワークがいくつあるか、またノードはどの群に属しているか確認できる。前はこの関数を使い Excel でネットワークを分割していたが今回は R で何とかしてみるか。。。

decompose.graph 関数

　分断されたネットワークを連結成分単位に分割するコードを書こうかと考えたが、ひょっとしたら 既に関数があるかもと思い約300ページある igraph のマニュアルを確認してみることにした。まずは目次から。すると decompose.graph 関数というものがあった。decompose = 分解 ということは？？？　ひょっとして発見した！？　ということでさっそく試す。（dで始まる関数だったの結構早く見つかって良かった）

> dcg <- decompose.graph(g) # 連結成分でグラフを分割する modeは"weak"みたい
> length(dcg) # グラフの連結数
[1] 8
> dcg
[[1]]
IGRAPH DN-- 400 752 -- 
+ attr: name (v/c)

[[2]]
IGRAPH DN-- 9 8 -- 
+ attr: name (v/c)

：（中略）

[[8]]
IGRAPH DN-- 3 2 -- 
+ attr: name (v/c)

ビンゴ！　ネットワークは8つに分割されていて、各ネットワークはそれぞれ配列っぽく格納されているじゃないか。これは使える！

非連結ネットワークのコミュニティを抽出する

　ということで、非連結ネットワークに対してコミュニティ抽出するコードを書いてみた。私のＲ経験値ではこんなコード。（この例では焼きなまし法だが他の手法でも同様）

# おまじない
ls(all=T)
rm(list=ls(all=TRUE))

# パッケージの読み込み
library(igraph)

# ワーキングディレクトリの変更
setwd("d:/R/igraph") # 作業ディレクトリの変更

# ネットワークの読込み
d <- read.table("エッジデータ.txt", header=T) # データフレームで読込み
g <- graph.data.frame(d[,2:3], directed = TRUE) # データフレームからグラフオブジェクトに変換
g <- simplify(g, remove.multiple = TRUE, remove.loops = TRUE) # 多重エッジや自己ループの削除

# 連結成分でネットワーク
dcg <- decompose.graph(g)

# 連結ネットワーク単位にコミュニティーと中心性を計算する
sp.df.all <- c() # コミュニティ分割結果データフレーム
for (i in 1:length(dcg)){
  set.seed(1) # シードを固定
  sp <- spinglass.community(dcg[[i]]) # 焼きなまし法 グラフラプラシアンを使ってQ値が最大となるような分割を探す
  sp.df <- cbind(i, as.data.frame(sp$names), as.data.frame(sp$membership)) # データフレーム形式に変換
  sp.df.all <- rbind(sp.df.all, sp.df) # コミュニティーに結果を追加
}
colnames(sp.df.all) <- c("i", "NodeName", "Community") # 列名変更
write.table(sp.df.all, "Community.txt", row.name=F, col.names=T, sep="\t", quote=F, append=F) # ファイル出力

焼きなまし法によるコミュニティ抽出結果

　上のコードを使い焼きなまし法によるコミュニティ抽出を行い冒頭のネットワーク図のノードを色分けするとこのようになった。（描画はcytoscapeにノード属性として結果ファイルを食わせて色分け）

でたぁっ　初回にしては結構いい感じ～～　ここまできたら中心性指標から特徴的なノードも見てみるか。

媒介中心性から特徴的なノードを確認

　とり急ぎ媒介中心性を使い特徴となりそうなノードも確認。媒介中心性の指標値が大きいノードを赤色にしてグラデーションをかけてみた。（描画はコミュニティ同様cytoscapeにノード属性として結果ファイルを食わせて色分け）

betweenness(g, directed=TRUE, weights=NULL) # 媒介中心性

中心性指標は分断されているネットワークでも計算できる（前のバージョンでは計算できなかった関数があったような気もするが。。。）が、中心性指標によってはコミュニティ抽出同様、連結成分ごとに計算しなければいけない場合もあるかも。。。その時はコミュニティ抽出の for ループの中で計算しちゃえば良いってことだ。

中心性指標を比較してみる

　分断されているネットワークの場合には指標値が異なる結果を返す中心性指標がある為、分断されていない最も大きなネットワークで比較。とりあえず相関。

i <- 1 # 分断されていない最も大きなネットワーク
cl <- closeness(dcg[[i]], mode="out", weights=NULL, normalized=TRUE) # 近接中心性
dg <- degree(dcg[[i]], mode="out") # 次数中心性
ev <- evcent(dcg[[i]], directed=TRUE, weights=NULL)$vector # 固有ベクトル中心性
pr <- page.rank(dcg[[i]], directed=TRUE, damping=0.85, weights=NULL)$vector # PageRank
bp <- bonpow(dcg[[i]], exponent=0.2) # ボナチッチのパワー中心性
bw <- betweenness(dcg[[i]], directed=TRUE, weights=NULL) # 媒介中心性
nc <- cbind(cl, dg, ev, pr, bp, bw) # 中心性指標まとめ

library(psych)
pairs.panels(as.data.frame(nc), smooth=FALSE, density=FALSE, ellipses=FALSE, scale=TRUE) # 散布図

やはり固有ベクトル中心性とPageRankには強い相関がある。他は。。。ふ～ん。。。

コミュニティ抽出 焼きなまし法と辺媒介性の比較

　最期にコミュニティ抽出で比較した焼きなまし法と辺媒介性の2つを比較した結果。辺媒介性ではコミュニティが138もある。（対象ノード数は400）　dendPlot 関数で確認してみると、げげ～　多くのノードがまとまった比較的大きなコミュニティーが数個と、ほぼ１つで形成されているだろうノードが数多く存在する結果となっているようであった。今回分析したネットワークの特徴？　有向グラフだから？　ということで今回は焼きなまし法を採用した。パフォーマンスも良いということなので。

> i <- 1 # 分断されていない最も大きなネットワーク
> sp <- spinglass.community(dcg[[i]]) # 焼きなまし法 グラフラプラシアンを使ってQ値が最大となるような分割を探す
> eb <- edge.betweenness.community(g) # 辺媒介性 媒介性の高い辺を削除して分割
> max(sp$membership) # 焼きなまし法のコミュニティ数
[1] 16
> max(eb$membership) # 辺媒介性のコミュニティ数
[1] 138
> dendPlot(eb)

まずはNon-overlapping手法を実践。今度はOverlapping手法とか実践してみながら今のデータ分析に有効か確認しよう。

参考
・ネットワーク分析をもうちょっと勉強 - でたぁっ 感動と失敗の備忘録

ファイル出力でちょっとハマった

　データフレームをファイルに出力する際にちょっとだけハマってしまった。
今までは前に作ったRコードをコピペしていたので無意識に、write.csv 関数を使っていた。それに気づかず、write.table 関数を書いていた。結果を見るとどうもおかしい。コンソールから help(write.table) でマニュアルを見ながら引数を加えた。結果こんな感じで個人的によく利用する形式になった。こんな時にしかマニュアル見ないので今のうちにメモ。

write.table(d, "write.txt", quote=F, sep="\t", row.names=F, col.names=T)

利用できる引数はこちら

quote	TRUEの場合は二重引用符で囲まれる
sep	デリミタを指定 タブは"\t"
row.names	行ヘッダーの出力有無 TRUEで出力
col.name	列ヘッダーの出力有無 TRUEで出力
na	欠損値に対する文字列
eol	行末の改行コード
dec	小数点に使用する文字
append	TRUEは追記される FALSEは上書き
qmethod	文字列にある二重引用符に対処する方法を指定する文字列
fileEncoding	エンコーディングを指定
ファイル名	"" 空白にするとコンソールに出力

以外に、ファイル名出力で空白にしとくとコンソールに出力されるので出力前の確認に便利かも。

ちなみに、今までよく使っていたcsv出力はこんな感じ。

write.csv(d, "write.csv")
write.csv(d, "write.csv", quote=F, row.names=F)

write.csv2 という関数もあるけど何が違うのだろう。。。
MASSパッケージには、write.matrix 関数もあるみたい。

ネットワーク分析をもうちょっと勉強

　マーケティング施策や商品開発の糸口を発見する為に購買履歴データでネットワーク分析を行っているのだが、今年に入りノードやエッジの数が多い隣接行列を取り扱うようになってきた。といってもまだまだビッグデータとはいいがたいが。。。cytoscapeにくわせOrganicレイアウトで表示するとこんな感じです。（他の設定はデフォルト）

今後はもっとデータ量が多くなりそうなので、ネットワーク分析をもうちょっと勉強。

　これまで参考にしていたのはRで学ぶデータサイエンスの#8ネットワーク分析。R+igraphで中心性指標を計算したり、コミュニティの抽出などはさくさくできる。会社の先輩に話をすると、もっと勉強になる資料があるよと九州工業大学の竹本先生の資料を教えてくれた。2013/3/8に開催された第2回 Rでつなぐ次世代オミックス情報統合解析研究会で発表された資料らしい。タイトルは「R+igraphではじめる生物ネットワーク解析」です。拝見させて頂くとRで学ぶデータサイエンスシリーズには記載されていないigraphの使い方、linkcommというパッケージを知ることができた。個人的に新しく知ることができた内容はこんな感じ。

・データフレームからigraphで扱うグラフに変換できる（逆もできる）
・くわせたデータの多重エッジや自己ループを削除できる
・最短パスを見つけ描画できる
・コミュニティ抽出にはスペクトル最適化法や焼きなまし法など他にも手法がある
・ネットワークの比較（和、積、差）を行うことができる
・コミュニティ抽出にはオーバーラッピングコミュニティーがある（linkcommパッケージ）

これはさっそく試してみるしかない！

個人的に新たに知ったigraphパッケージの関数

　まずはデータフレーム形式で読込みから。データは竹本先生のページにあった演習用データをダウンロードさせて頂きました。URLはこちらです。R seminar on igraph - supplementary information - Kazuhiro Takemoto

library(igraph)
d <- read.table("eco_EM+TCA.txt") # データフレームで読込み
g <- graph.data.frame(d, directed = FALSE) # データフレームからグラフに変換
g <- simplify(g, remove.multiple = TRUE, remove.loops = TRUE) # 多重エッジや自己ループの削除

df <- get.data.frame(g) # グラフオブジェクトからデータフレームに変換

感動！　これでRDBから抽出したデータをすぐ食わせられる。RODBC使っちゃおう。また、その逆もできそうだ。cytoscapeにくわせるエッジリストも簡単じゃん。ちなみに今までだとこんな感じでRにくわせていた。（RDBからデータを抽出しExcelとか使いながら整形）

g <- graph(c(1,2,1,3,1,4,1,5,1,9,2,3,2,4,3,4,5,6,5,7,5,9,6,7,6,8,7,8),n = 9,directed = FALSE)

次は最短のパスを見つけ描画してみよう。最短パスを知る関数はRで学ぶデータサイエンスにもあったが描画の工夫はなかったと思う。
例は2とのノード"D-Glucose"、"2-Oxoglutarate"間の最短パスみたいだ。

sum <- get.all.shortest.paths(g, "D-Glucose", "2-Oxoglutarate", mode="out")
V(g)[sum$res[[1]]]
E(g)$color <- "grey"
E(g, path=sum$res[[1]])$color <- "red" # ※最短経路が複数ある場合は1つ目
E(g, path=sum$res[[1]])$width <- 3
plot(g, vertex.label=V(g)$name, vertex.size=12, layout=layout.fruchterman.reingold)

いい感じ。（この例ではsum$resには2つの結果があったので1つめの結果を描画）

次はコミュニティ抽出法。はじめの2つは知っていたけど他にもこんなにありました。

eb <- edge.betweenness.community(g) # 辺媒介性 媒介性の高い辺を削除して分割
rw <- walktrap.community(g) # ランダムウォークに基くアルゴリズム
fg <- fastgreedy.community(g) # 貪欲アルゴリズム Q値が高くなるいようにペアをまとめていく
le <- leading.eigenvector.community(g) # スペクトラル最適化法 グラフラプラシアンを使ってQ値が最大となるような分割を探す
ml <- multilevel.community(g) # 多段階最適アルゴリズム
sp <- spinglass.community(g) # 焼きなまし法 グラフラプラシアンを使ってQ値が最大となるような分割を探す
op <- optimal.community(g) # 全探索（時間がかかるため小さなネットワークのみ有効）

$membershipの最大数（コミュニティ数）を比較してみるとこんな感じ。

eb	rw	fg	le	ml	sp	op
4	5	5	5	4	6	6

手法が違うのでコミュニティ数も変わるのは当然だが、資料によると計算時間とパフォーマンスはトレードオフになっているとのこと。ふむふむ。焼きなまし法（アニーリング）なんかいいのかな？

焼きなまし法で抽出したコミュニティ毎に色を変え描画するとこんな感じになった。大きなネットワークで行ったらどうなるのだろう。美しさならcytoscapeなどのソフトに食わせるしかないのかな？それは今度。

V(g)$color <- sp$membership
plot(g, vertex.size=12, vertex.label=V(g)$name, layout=layout.fruchterman.reingold)

コミュニティの状況はデンドログラムでも確認できるのだが、dendPlot関数を使うこともできるようだ。dendPlotを使うと綺麗な罫線で囲んでくれるから見やすい。fastgreedy.communityの結果はこんな感じ

dend <- as.dendrogram(fg)
plot(dend) # 罫線は出ません

dendPlot(fg) # 罫線が出るデンドログラム

※焼きなまし法だとNot a fully hierarchical community structureといったエラーになる。完全な階層になっていないの? 階層構造になっていないので表示できない。

ネットワークの比較

　竹本先生のページR seminar on igraph - supplementary information - Kazuhiro Takemotoを最期まで読むと、最期にネットワーク比較なるものが。。。igraphでできたの～って感じです。
g1、g2　といったグラフに対して行う場合にはこのように使います。

graph.union(g1, g2) # 和
graph.intersection(g1, g2) # 積
graph.difference(g1, g2) # 差 ※ただし |g1| > |g2|

igraphのマニュアルには、graph-operators に記載されていました。目次には union とかで記載されていません。300ページ近くもあるからただでさえ把握困難なのに。。。
ノード名で行いたい場合は～～.by.name関数がありました。

graph.union.by.name(g1, g2) # 和
graph.intersection.by.name(g1, g2) # 積
graph.difference.by.name(g1, g2) # 差

差の結果はこんな感じです。データは適当に作っております。
plot(g1)

plot(g2)

plot(graph.difference.by.name(g1, g2))

こんな図では確認しにくいですが g1 - g2 ができていました。g1からg2にあるエッジが消える感じです。隣接行列の演算をしているのでしょう。これは時系列変化などでの活躍間違いなしだ。

オーバーラッピングによるコミュニティー抽出

　最後にコミュニティーはオーバーラップしているとい考え方を取り入れた分析手法。linkcomパッケージを使うのか。これは初めてのパッケージ。4種類ある手法の2種類が実装されているみたいだ。さっそくlinkcomパッケージをインストールして確認。

まずは、LinkCommという階層クラスタリング手法

library(linkcomm)
g <- read.table("eco_EM+TCA.txt")
lc <- getLinkCommunities(g)

デンドログラムがでた。

階層構造と分割密度なるものが確認でき、lc$pdmax を見ると最大は0.833のようだ。
任意の密度でネットワーク図を描画するにはこんな感じ

lc <- newLinkCommsAt(lc, cutat=0.84) # cutstに任意の分割を指定
lc$clusters # クラスタリング結果を確認
lc$nodeclusters # ノード名でクラスタリング結果を確認
plot(lc, type="graph") # ネットワーク図を描画

階層クラスタリング手法の結果はこんな感じでcytoscape用エッジ属性をファイルに出力できるようだ。cytoscapeも勉強不足。

linkcomm2cytoscape(lc, interaction="pp", ea="edge_atrib.ea") # cytoscape用エッジ属性の出力

オーバーラッピングで使える関数は１つOCG(Overlapping Cluster Generator)という手法。これは貧欲アルゴリズムでLinkCommよりも良いパフォーマンスを示すとのこと。

ocg <- getOCG.clusters(g)
plot(ocg, type="graph") # ネットワーク図を描画

でたぁっ。階層クラスタリング手法でも同様だがノードの色がすごい。いくつのコミュニティーに属するかでノードが円グラフみたいになっている。これはcytoscapeでもできるのかな。。。今はそんなことはいいか。後で考えよう。

また、複数コミュニティーに属するノードはテーブルで描画することもできる

plot(ocg, type="members") # 複数コミュニティーに属するノードを見る

すげ～～～～～～～～～。竹本先生どうもありがとうございました。m(__)m　本当に感謝です！
先生の資料を教えてくれた先輩もありがとう！
自分なり今までよりネットワーク分析ができそうな気になってきた。バックグラウンドや関数の特徴などを理解しながらさっそく実務でいろいろ使ってみよう。CRANにあるigraphのマニュアルって300ページもあるから後回しにしていたのも反省。やっぱりマニュアル読まないとダメだ。

参考
・Rで学ぶデータサイエンスの#8ネットワーク分析

・R+igraphではじめる生物ネットワーク解析

・竹本先生のページ
R seminar on igraph - supplementary information - Kazuhiro Takemoto