Rによるオープン・データの可視化（2）：実践！ Rで学ぶ統計解析の基礎（4）（3/3 ページ）

[柏野雄太，＠IT]

Wikileaksのアフガン戦争ダイアリー

　最近のオープン・データ的イベントとして大きかったものは、今年の7月25日にWikileaksが公開した、アフガン戦争ダイアリー（Afgan War Diary）だと思います。これは、いままで数々の機密ファイルや機密映像の暴露を行ってきた民間団体Wikileaksが、米軍からの内部情報としてまとまったレポートを公表したものです。このレポートは、アフガニスタンの駐留米軍の戦闘をイベントごとに記録したもので、2004年初から2009年末の6年間に渡る7万6900件のレコードから構成されています。

　米国政府の発表では、この「ダイアリー」自体には軍事的・諜報的な価値がないとされていますが、しかしこれだけ大量の生データが、現在も戦闘及び駐留が継続される作戦から流出したことは前代未聞です。そして、1つ1つの情報の断片だけでは見えないことも、ある切り口からまとまった量のデータをみたときには、新たな情報が見えてこないとは限りません。ただし、今回は単にアフガン戦争ダイアリーのデータをRに取り込み、それれを図示・可視化し、簡単な統計的な解析をするというだけに留めたいと思います。つまり、比較的大きなオープン・データを入手し、取り込み、図示し、解析するという過程をステップバイステップで示したいと思います。

WikiLeaksは機密情報の暴露を行う民間団体だ

　データの解析方法はニューヨーク大学の政治学の大学院生であるDrew Conway氏の方針を採用することにします。実は、筆者もこのWikileaksのデータ公開があった直後に、データをダウンロードしてRによる解析を試みていましたが、Conway氏は筆者より早く解析結果とコードを自身のブログで公開しました。そしてその後もConway氏は精力的に解析活動を続け、以下のような一連の記事を生み出しています。

●Zero Intelligence Agents

• Wikileaks Afghanistan Data
• Benford’s Law Tests for Wikileaks Data
• Wikileaks Attack Data by Year and Type Projected on Afghanistan Regional Map
• Animated Heatmap of WikiLeaks Report Intensity in Afghanistan
• In Search of Power-laws: WikiLeaks Edition

　ここでは、Conway氏の解析をベースに、Wikileaksのデータを図示・可視化するとともに、筆者も同時期にやっていたベンフォードの法則への割り当てとカイ二乗適合度検定をしてみたいと思います。利用するコードの前半のデータ操作の部分は、筆者のものより洗練されているConway氏のコードを使用させていただき、後半の解析部分は筆者のコードをConway氏のコードにマージする形で示したいと思います。この記事を書く際に、Conway氏とはメールで議論させてもらいました。Conway氏に感謝いたします。

　何よりもまずデータを読み込まなければなりません。この「アフガン戦争ダイアリー」は以下にあり、HTML、CSV、SQL、KMLファイルなどのファイルフォーマットで、誰でもダウンロードが可能な状態になっています。

Wikileaks Afghan War Diary, 2004-2010 http://www.wikileaks.org/wiki/Afghan_War_Diary,_2004-2010

　これから先以降は、この中のCSVデータをダウンロードして、圧縮を伸長した状態であると仮定します。伸長するとafg.csvというファイルができるはずです。サイズは73.9MBと、CSVファイルとしてはかなり大きなものです。そしてそのファイルが“/home/my/”ディレクトリにあると仮定します。

　setwd関数でアフガン戦争ダイアリーのCSVファイルがあるディレクトリをワーキングディレクトリとし、データを読み込みます。74MBという大きなサイズのCSVですので、読み込みに多少時間がかかるかもしれません。

> setwd("/home/my/")
> afg<-read.csv("afg.csv",stringsAsFactors=FALSE)
> head(afg)
 以下にエラー do.call("data.frame", rval) :  変数名は256バイトが上限です

　最後の行は変数名に256バイトを超えるオブジェクトがあるので、head関数が失敗していました。そこでデータの構造をみるためにstr関数を利用すると、「ダイアリー」のデータは32カラムで7万6910行のデータフレームになっていることが分かります。

> str(afg)
'data.frame':	76910 obs. of  32 variables:
 $ D92871CA.D217.4124.B8FB.89B9A2CFFCB4                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    : chr  "C592135C-1BFF-4AEC-B469-0A495FDA78D9" "D50F59F0-6F32-4E63-BC02-DB2B8422DE6E" "E3F22EFB-F0CA-4821-9322-CC2250C05C8A" "4D0E1E60-9535-4D58-A374-74367F058788" ...

　このデータフレームのカラムにカラム名を追加します。そのために「ダイアリー」のデータ構造について調べる必要があります。このために、http://wardiary.wikileaks.org/にあるボールド体による“the structure of the report”を参考にしましょう。

　これを眺めてみると、まず俯瞰的な描像を手に入れたいと思うはずです。そのために、「6年間に渡って毎月どのくらいのイベントがどこで生じたのか」ということを図示することにします。それに必要なカラムは、

• イベントが起きた日付データのDateOccurred
• 場所を表すRegion
• イベントの種類（友軍への攻撃、敵への攻撃など）を表すAttackOn（これはAffilicationでも代替できそうです）

くらいでしょうか。

　ちなみに、この「ダイアリー」の省略コードにある、KIAは“Killed in Action”, WIAは“Wounded in Action”だそうです。そして、FriendlyWIAは友軍が負傷、CivilianKIAは一般市民が死亡人数を表すそうです。これで、攻撃の種類による負傷者数や死者数を知ることができそうです。その他の略字は、先ほどのリンク先の説明を参考にしてください。

　さて、「ダイアリー」のデータ構造が分かったところで、読み込んだデータフレームに、上記のサイトにあるカラム名を追加します。

> colnames(afg)<-c("ReportKey", "DateOccurred", "Type", "Category", "TrackingNumber", "Title", "Summary", "Region", "AttackOn", "ComplexAttack", "ReportingUnit", "UnitName", "TypeOfUnit", "FriendlyWIA", "FriendlyKIA", "HostNationWIA", "HostNationKIA", "CivilianWIA", "CivilianKIA", "EnemyWIA", "EnemyKIA", "EnemyDetained", "MGRS", "Latitude", "Longitude", "OriginatorGroup",  "UpdatedByGroup", "CCIR", "Sigact", "Affiliation", "DColor", "Classification")
> head(afg)
                             ReportKey        DateOccurred                Type            Category
1 C592135C-1BFF-4AEC-B469-0A495FDA78D9 2004-01-01 00:00:00     Friendly Action Cache Found/Cleared
2 D50F59F0-6F32-4E63-BC02-DB2B8422DE6E 2004-01-01 00:00:00    Non-Combat Event          Propaganda
3 E3F22EFB-F0CA-4821-9322-CC2250C05C8A 2004-01-01 00:00:00        Enemy Action         Direct Fire
4 4D0E1E60-9535-4D58-A374-74367F058788 2004-01-01 00:00:00     Friendly Action Cache Found/Cleared
5 84A769E7-73B4-4B30-AD36-06497F766F15 2004-01-02 00:00:00 Suspicious Incident        Surveillance
6 B81BE975-9081-4D19-A8FE-C2DB0C18A095 2004-01-02 00:00:00 Suspicious Incident        Surveillance
        TrackingNumber                      Title
1 2007-033-004738-0185 CACHE FOUND/CLEARED  Other
2 2007-033-010818-0798          PROPAGANDA  Other
3 2007-033-004042-0850         DIRECT FIRE  Other
4 2007-033-004738-0279 CACHE FOUND/CLEARED  Other
5 2007-033-010701-0111                SURV  Other
6 2007-033-010701-0220                SURV  Other
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               Summary
1                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       USSF FINDS CACHE IN VILLAGE OF ....

　コードではcolnamesを利用してカラム名を付け、headで最初の6行を取り出しました。ここでは長くなりますので、Summaryの先頭部分以降は省略しています。Summaryに長い文字のレポートが入っていることが分かります。どうやら、最初のhead関数がミスをしたのはここをタイトルとして読み込もうとしたことが理由のようです。

　さて、ここで「6年間に渡って毎月どのくらいのイベントがどこで生じているか」という図示をするために、日付データや地域データを整形します。日付データはas.Dataとしてdateオブジェクトに変換します。幸い、このCSVのDateOccurredにあるデータはUnix形式の日付データなので、そのまま素直にはいります。また、地域を表すRegionは因子として扱いたいのでas.factorでfactorオブジェクトに変換します。

> afg$DateOccurred <- as.Date(afg$DateOccurred)
> afg$Region<-as.factor(afg$Region)

　さて、ここまででafgオブジェクトには必要なデータがすべてそろっていることになります。それでは改めて問いますが、「6年間に渡って毎月どのくらいのイベントがどこで生じているか」を図示するためにはどうすればいいでしょう？ まず、読み込むデータフレームはafgであり、日付のx軸はDataOccrredである。そして、これを30日ごとのヒストグラムとして表す。しかし、ヒストグラムはイベントの種類AttackOnもしくは Affiliationによって色分けした積み上げグラフにし、しかも、これらすべてについて地域を表すRegionごとに書き出したい。これを実行するのは大変なようです。しかし、ggplot2を使用し、このまま「グラフィック文法」に沿ってコーディングすれば、この作業は、いとも簡単にできてしまいます。

　以下にそのコードを示します。

> library(ggplot2)
> event <- ggplot( afg, aes(x = DateOccurred) ) + stat_bin( aes( y = ..count.., fill = AttackOn ), binwidth=30 ) + facet_wrap( ~Region ) + 
    opts( title = "Monthly Event Report per Region and Target" ) + xlab( "Date") + ylab( "Report Counts" ) + 
    scale_fill_manual( values = c( "darkred", "darkblue", "darkgreen", "orange" ), name = "AttackOn" ) + scale_x_date( major = "years", minor = "months" )
> ggsave("wikileaks_event.png", plot = event, width=10, height=5)

　ggplotでデータフレームとx軸を指定し、stat_binでヒストグラムを描きます。stat_binのbinwidthは30日ごとのヒストグラムにしたいので30としています。また、stat_binのaes関数には、y軸はカウントデータを表す...cont...を、色分けはAttackOnで行うので、それをfillの引数としました。最後のRegionごとに書き出すのはggplot2の「十八番」であるファセットグラフを利用します。facet_wrap関数がそれです。

　こうして書き出したグラフは次のようになります。

　これを見ると分かるのは、戦闘はすべての地域で行われるのではなく、タリバン勢力の支配地域となっている東および南に多いということです。また、各々のファセットグラフの時間変動をみると季節性の変動があることが分かります。アフガニスタンは冬が厳しいためです。さらに、時間を経るごとに戦闘数が多くなっていることが分かります。つまりアフガニスタン戦争が終結した後も、内戦状態がますます激しくなっている状況が想像できます。先ほど、2011年7月から米軍は撤退を開始するという報道がありましたが、果たしてこの状態で米軍は来年から撤退を開始できるでしょうか？ この解析からは疑問符が付くように思えますが、それは杞憂でしょうか？

　まだ、解析は始まったばかりなのですが、紙面の都合上、次に進みたいと思います。次に検証を試みたいのは、果たしてこのデータは本物なのかどうか、ということです。もちろん、すでに米国政府はこのデータは正式なものであると認めていますので、それを信じるならばこれ以上の解析は必要ありません。しかし、ここではその判断を一旦留保しまして、データが人為的に作成された可能性があるかという点について調べたいと思います。

ベンフォードの法則で検証

　ここで「ベンフォードの法則」が登場します。「ベンフォードの法則」は、データが人為的に作成されたかどうかを検証するための有用なツールです。その利用の概要を一言で述べますと、自然界に存在する多くの数字や数値に対して最初の1桁だけを取り出して数えた場合、それがある法則（log10((1+x)/x）に従っているということを利用するのです。

　有名なのが、経理データで不正経理をしているとおぼしき数値の先頭桁だけを取り出しヒストグラムにした場合、ベンフォードの法則に従っていない場合は、人がデッチ上げをした可能性がありそうだ、というものです。実際に、不正な会計操作を発見した実例が以下のエッセイに書いてあります。英語が億劫でなければ読んでみてください。

I've Got Your Number

　その他の詳細は、ウィキペディアを参考にするか、筆者も昔個人ブログに記述したことがあったので、それを参考にしてください。

• ウィキペディア: ベンフォードの法則
• wrong, rogue and log: Benfordの法則 -その1-
• wrong, rogue and log: Benfordの法則 -その2-

　さて、このベンフォード則の「ダイアリー」への適用ですが、先ほどの図示した月ベースではデータ数が少なすぎ、日ベースではデータの種類が少なすぎるので、都合がよくなさそうです。ここでは、週ごとにデータをまとめることにします。実は、これはすごく簡単で、すでにRのdateオブジェクトになっているので、format.Date関数を利用してフォーマットを指定するだけです。

> weeklyCount<-cbind(table(cbind(format.Date(afg$DateOccurred,"%Y %W"))))
> head(weeklyCount)
        [,1]
2004 00   12
2004 01   27
2004 02   43
2004 03   31
2004 04   27
2004 05   16

　最終的にデータフレームで利用するために縦ベクトルであって欲しいので、cbind、table、cbindという姑息な手段を利用しました。さて、このweeklyCountの先頭桁を取り出すには、工夫が必要そうです。幸いC言語の数値フォーマットに変換する関数formatCを利用すれば、例えば、

> formatC(12, format="e")
[1] "1.2000e+01"

のようにできるので、これを文字列の部分文字列取り出し関数substringを利用すればよさそうです。以下にコードを示しました。

> pullLeaddigit<-function(x) { as.numeric(substring(formatC(x, format = "e"), 1, 1)) } 
> head(pullLeaddigit(weeklyCount))
[1] 1 2 4 3 2 1

　実際に適用すると、良好なようです。さて、このpullLeaddigit関数をつかって、weeklyCountの先頭桁を取り出し、それを散布図プロットをして、ベンフォードの法則の理論曲線と比べてみましょう。

> digits <- table(pullLeaddigit(weeklyCount))
> digits <- cbind(digits)
> digits <- as.data.frame(digits)
> colnames(digits) <- "digitCount"
> digits$digitCount <- digits$digitCount/sum(digits$digitCount)
>
> dbenford <- function(x) { return(log10((1+x)/x)) }
> 
> benford <- ggplot(digits, aes(x=1:9, y=digitCount)) + geom_path(aes(colour = "observation")) + geom_point(aes(colour = "observation")) + stat_function(fun = dbenford, aes(colour = "theory")) + scale_colour_manual( values = c("observation" = "orange", "theory" = "blue")) + xlab("Digit") + ylab("Count Density") + opts(title="Applying Benford's law to Wikileaks Data") 
> ggsave("wikileaks_benford.png", plot=benford, width=7, height=7)

　以下が出力されるグラフです。

　オレンジがデータ、青がベンフォードの法則の理論曲線です。これをみると分かるように、1が理論曲線に比べて少ない気はしますし、2は多すぎるような気がしますが、大局的には良く当たっているように思います。

　さて、それが本当に統計学的に言えるかどうかを調べるには、どうすればよいでしょうか？ そうです、第1回、第2回で取り上げた統計的検定のカイ二乗適合度検定を用いればいいのです。この場合の帰無仮説は、「アフガン戦争ダイアリーの週データの先頭桁カウントは、ベンフォード則の分布と同じである」です。chisq.testを利用して早速調べてみましょう。

> chisq.test(digits$digitCount, dbenford(1:9))
	Pearson's Chi-squared test
data:  digits$digitCount and dbenford(1:9) 
X-squared = 72, df = 64, p-value = 0.2303
 警告メッセージ： 
In chisq.test(x = digits$digitCount, y = dbenford(1:9)) :
   カイ自乗近似は不正確かもしれません

　結果を見ると分かるように、p-value = 0.2303である、つまり有意水準5％でも帰無仮説は棄却できないということです。これは言い換えると、「アフガン戦争ダイアリーの先頭桁カウントは、ベンフォード則の分布と同じである」ことを否定するのは統計的に言えなさそうだということです。以上、ベンフォード則とのグラフをみても、カイ二乗適合度検定をみても、アフガン戦争ダイアリーは人が捏造したものでなさそうだ、ということです。

　以上、かなり重たい内容になってしまいましたが、このようにデータさえ手に入れれば、データ解析を自分の手でできることが感覚的に分かったかと思います。

次回について

　今回は、ggplot2の紹介、Webサービス経由でデータリポジトリを利用する方法、そして、それなりに大きなデータをどのように扱って可視化し、解析するか、ということについて実践しました。次回はオープン・データを図示・可視化する際の注意点をまとめると共に、統計学的な推測や予測についての初歩についてお話したいと思います。それでは、2週間後にまたお目にかかりましょう。

Index

Rによるオープン・データの可視化（2）

Page1
ggplot2という「デフォルト」グラフィックパッケージ

Page2
世界銀行のData Catalog

Page3
Wikileaksのアフガン戦争ダイアリー
ベンフォードの法則で検証
次回について

「実践！ Rで学ぶ統計解析の基礎」バックナンバー

• 大相撲のアノーマリー （2）
• 大相撲のアノーマリー （1）
• 夏の異常気象をオープン・データで確認
• インターリュード: TwitterとR
• Rによるオープン・データの可視化（2）
• Rによるオープン・データの可視化（1）
• Rで実践する統計的検定の初歩
• Rは統計解析のブッシュナイフだ