2012-12-10

R Advent Calendar 2012

+はじめに
久々の日記。
今年もR Advent Calendar 2012の季節です。仕事が立て込んでいて、すっかり日付変わるどころか朝になってしまい申し訳ありません。今回は、Japan.R #3でも発表したRevolution R Enterprise のお話をしたいと思います（使い回しとも言うのですが、ちょっと立て込んでいまして。。。）。

Revolution R Enterprise はR をベースにした商用のソフトウェアです。何が、Rと違うのかというと、データをHDDに格納してそれにモデルを適用できるので、巨大サイズのデータを扱えるという点があります。

取りあえず詳細は以下のスライドシェアを見て頂くとして。
http://www.slideshare.net/Hiro_macchan/japan-r3
http://www.slideshare.net/Hiro_macchan/tokyo-r25-hiromacchan
今日は、上記のスライドで触れられなかった点について触れたいと思います。

+rxTextToXdf()でFactorデータを取得した際のラベリングについて
rxTextToXdf()は、引数にテキストデータのパスと、HDD上に保存するXDFファイルパスを指定する事でRevoscaleR 上でデータを取り扱える様にしてくれます。値の型も結構色々使えるので便利で扱いやすいのですが、どうもこの関数、因子が出現した順番でlevelを設定するみたいで、このままGlmとかかけるとリファレンスになるlevelが変な感じになります。Xdfファイルの編集はrxDataStepXdf()関数に引数、transformsをリストで指定して行います。
今回の例では、

transforms = list(DayofMonth = relevel(DayofMonth, "1"), DayOfWeek = relevel(DayOfWeek, "1"))

で編集できます。relevelとか、普通の関数が使えて良かった。。。

+決定木の利用について
Revolution R Enterprise 6.1 から決定木がサポートされました。一回irisで使ってみましょう。~~そういえば、生まれて初めて決定木を使いますね。私。~~

Tree.a <- rxDTree(Species ~ Sepal.Length+Sepal.Width+Petal.Length+Petal.Width,  data = iris) 
Tree.a 
||<　
 結果はこんな感じです。

>|r|
File:   
Number of valid observations:  150 
Number of missing observations:  0 

Tree representation: 
n= 150 

node), split, n, loss, yval, (yprob)
      * denotes terminal node

1) root 150 100 setosa (0.33333333 0.33333333 0.33333333)  
  2) Petal.Length< 2.45 50   0 setosa (1.00000000 0.00000000 0.00000000) *
  3) Petal.Length>=2.45 100  50 versicolor (0.00000000 0.50000000 0.50000000)  
    6) Petal.Width< 1.75 54   5 versicolor (0.00000000 0.90740741 0.09259259) *
    7) Petal.Width>=1.75 46   1 virginica (0.00000000 0.02173913 0.97826087) *

うん、何となくうまく切れているんでしょうか？データサイズを変えたい場合、dataの引数をいじれば大丈夫です。便利ですね。

ちなみに、airlineデータでやってみると。。。

rxDTree(Cancelled ~ Month + DayOfWeek + Origin +AirTime, data = "air2008.xdf")

1) root 7009728 137434 0 (0.98039382 0.01960618) *

あらら。これは、キャンセルするかどうかをがっつり決めている因子が見つからなかったという認識で良かったんだろうか。。。ここら辺の、判別器による特性とかは抑えないと行かんなぁ。。。

+HPC Serverを用いた分散処理について
Revolution R Enterprise 6.1 では、Windows HPC Svr を利用した分散処理をサポートしている。ちょっと、手一杯なのでこれに関してはまた後日（年末くらいに）。

+最後に
今年一年、色々ありましたが少しずつ成長できた一年でした。それもこれも、色んな出会いを通じて刺激を受けて来れたおかげです。統計や、機械学習関係の出会いを授けてくれたR に感謝して、、、
皆様、よいクリスマスを。。

2012-03-23

【R】Tokyo.R#21で発表してきた。

　去る3月10日、都内某所で開催されたTokyo.R#21にて、発表をしてきました。テーマは「Rで学ぶ観察データでの因果推定」です。
　私の現在の研究テーマである医療のアウトカムリサーチを行う上で、バイアスの影響を除いた効果推定を行う事は極めて重要です。そこで利用される手法に因果推定が挙げられ、今回発表したのは、その一例です。一般的回帰分析と傾向スコアでの分析、そして、傾向スコアと一般的回帰分析が前提においている「強く無視できる割り付け」仮説が崩れた場合の対応としての操作変数を用いた推計を記載しました。

内容は下記SlideShareに載ってます。もしご興味があれば。

「Rで学ぶ観察データでの因果推定」

多分、計量経済やってる人からみると不十分なところが多いと思うけど、OPENにする事で指摘をもらえると良いなぁ。。。

2012-01-27

【医療】【研究】Janeというサービス

Jane
http://www.biosemantics.org/jane/

twitter上で流れてきた、タイトル入れるとおすすめの雑誌投稿先を教えてくれるサイト。

論文投稿する側からすると、狙う雑誌を何にするかは結構重要なポイントで、慣れてくればいいんだろうけど、僕みたいなひよっこは、結構迷う。そういう意味で、狙い所を探すのに使えそう。

あと、eigenfoctorってここで初めて知った。無料で使えるIFみたいな物かぁ。使ってみようかな。

関連する著者とかも探せるから、特定分野レビューとかに対する、使い勝手は良さそうだ。Web of scienceとかは、有料だったりでハードル高いこともあるので、無料で使えるこう言うサイトは重宝したりする。

仕組みは恐らく、pubmed APIにクエリ投げて、そこから得られた情報にいくつかの他ソースを紐づけているんだろう。検索のロジックは色々噛ませていそうだけど。。。

仕組みもいいけど、ニーズにそったアイディアが秀逸。

論文投稿する際には是非活用したい。。

2012-01-10

Google Fusion Tables

はじめに
ちょっと、使用でGoogle Document を触る事がありまして、その中にCreate Table って項目を見つけたわけです。説明を簡単に読んでみると、結構お手軽にデータの可視化・共有が出来そうな感じ。一度使ってみました。

使用データは、以前こちらで紹介した全国のDPC病院の住所地リスト（一部カラムを削除してスリムにしています。）を使用します。こちらのデータをGoogle Fusion TablesでGoogle Map上に可視化してみましょう。

データのインポート
データのインポートは、ローカルからのアップロード、Google Docsからのインポート、外部URLからの取得などが使える様です。

データの確認
インポートしたデータはTable形式で保管されるため、データのチェックなどはこちらで。

Geocoding
ジオコーディング処理は、Google Map APIがやってくれます。
カラム設定で住所をLocationカラムに設定します。

Visualize
可視化設定をMapに設定します。Configure Styleでポイントスタイルの設定が可能です。今回は病床数カラムを元にポイントを色分けしました。

完成品
https://www.google.com/fusiontables/DataSource?snapid=S355266-ngK

未完成部分
１）実は、データはすでにジオコーディング済みだったので、緯度経度をLocationカラムとして指定して利用したかったのですが、うまく行きませんでした。「Two column location:」あたりで設定できそうだったんだけど。。。生じた事象としては「全医療機関中一部の医療機関しかポイントとしてインポートされない。」という現象で、原因などは不明です。誰か、使用方法ご存じでしたら教えて下さい。
Geocodingは結構負荷が高いらしく、余りGoogleさんの手を煩わせたくないので。。。

２）どうも、Google Visualization API が後ろで動いている様子なので、RのgoogleVisパッケージとかを使用して処理をスクリプト化出来そうなんですが、詳しく調べていません。また、機会があったら調査してみます。

３）症例実績とかと紐付けしたいですね。。

終わりに
使った感じ、まだ一部うまく使いこなせない所はありましたが、データの可視化と共有という意味で結構強力なツールたりえるんではないかなと思います。みなさんも、ご興味がありましたら是非。

追記(20120112)
dichikaさんからのご指摘で、Google document の設定が英語設定になっていないと、表示されないらしいです。補足ありがとうございました。*1

*1:http://d.hatena.ne.jp/dichika/20120112

2011-12-31

もうすぐ年が明ける。

今年を振り返ってみたいと思う。

といっても、今年春は、震災の記憶しか無い。通常業務を行っていた都内の私の職場も、東京で経験したことのない揺れに襲われた。本当に立っていられないくらいの揺れは、阪神大震災の揺れを布団の中で感じた小学校以来だった（私は阪神大震災の直接の被災者ではありませんが、震源から遠く離れた私の実家でも、朝飛び起きるほどの揺れを感じました。）。院内からは、大きな被害情報は伝わってこなかったが、震源が「宮城県沖」と聞いたときは、血の気が引いた。震源があれだけ遠くて、あの揺れだったら、震源地はどうなっているんだろうと。。。

その後は、テレビやラジオから流れてくるニュースに一喜一憂する日々だった。１０日後には被災地にも医療支援で向かったし、募金も協力した。自衛隊や米軍が行き来する高速道路をひたすら運転して見た被災地の惨状は、今でも目に焼き付いている。でも、僕に出来たことは本当に微々たる者で、自分の力不足を感じた。

五月頃には、原発の状況もはっきりせずで世間の情勢はまだ落ち着かなかったが、自分にやれることをやる以外に何も出来ないことを知り、やるべき事を少しずつやっていくことにした。

その頃、TokyoRという、Rという言語の勉強会に参加する機会があった。今まで、全く縁の無かったプログラマーさんやエンジニアさんと話す機会は僕にとって、とても新鮮で、新しい研究を計画する上で重要なアイディアを与えてくれたし、今もまだ、いくつかネタを提供し続けてもらっている。
僕に出来ることは少ないものの、何らかの貢献はして行きたい。（このブログも、そういった活動の中で作ってみることにした。）

研究を進める上での基礎技術の習得をしていたら、夏は終わっていた気がする。たまに、自分に一切関係ない内容を勉強していく事は意外と大切だと感じた。
まぁ、夏場はサイクリングにあちこち出かけていたので、そっちに時間を食われていた感じは否めないが。。。

秋から冬にかけては、申請書や研究計画書をまとめることに心血を注いだ。僕をかつて指導して下さった先生からも、意外と好評かを頂き（基本的に、その先生からは怒られることしか無かったからうれしかったのを覚えている。）、論文化に向け舵を切った。

冬になったら、今度は研究とは別の業務が忙しくなった。研究は完全にスタッグしてしまった。人に何かをお願いして業務を回す難しさや、タスク管理の不十分さを感じた冬だった。

一年を振り返ると、春の震災を超えて、秋までは自分のやるべき事をなし、冬になるにつれて、諸々の障害が顕在化した。主に、タスク管理が不十分である事は危険なことと、人に何かをお願いする際に、他人の力量を十分把握した上で考えなければ読み損ねた分が自分に跳ね返ってくることを改めて認識した年だったと言える。

それでも、新たな人脈を構築し、新たな技術を習得し、新たな領域に足場を組む為の準備が出来つつある事は今年の進展だと思う。来年さらにそこに力を入れたい。

年始段階で考えていた、キャリア形成の事については、ここに記載する事は避ける。しかし、結局実績積まないとどうしようもないという当たり前の結論を一年かけて得ただけだった気がする。

来年は、さらに前に進めるよう、今年の失敗を踏まえて一年を頑張っていきたい。

最後に、今年一年お世話になった友人、同僚、先輩、家族、皆さんに心より御礼申し上げ、今年の総括とする。

２０１１．１２．３１　ブログ著者

2011-12-12

【R】R での行政データの利用（ＤＰＣ分科会データ）

　R は、統計方面に非常に強い環境です。私も、昔からちょこちょこ利用しています。今回、少しR を用いた記事を書く機会に恵まれましたので、R の紹介をかねて、1つアウトプットを出したいと思います。

　この記事はR Advent Calendar 2011 の12日目の記事となります。先にブログを書かれた皆さんの記事を読んでみと。。。何というかすごい。なんか、ばりばりモデル作ったり、すごいグラフが出てきたり、、、、遺伝子。。。？無理。こんなん無理。。そもそも、ブログ書くこと自体難しい。。。

　ということで、「１Ｒユーザー」として、R をどんな風に使っているかをご覧頂いて、今回の記事にしたいと思います。コードの不整合、意味わからないなどのご指摘、お待ちしております。（心が折れない程度でお願いします。）あと、これをR でやる理由がわからないというご指摘に関しては、、、察して下さい。

#以上前置き。

　皆さん、行政が提供している統計データを利用されたことってありますか？私のように、社会医学系の研究をしている人間はそういった行政データを使用する事も時々あるわけですが。。。。。これが、非常に使いにくい！

　なぜ使いにくいかを列挙していくと
・良くてExcel、悪いとPDF、最悪が報告書をスキャンした画像ファイル（PDF)
・データというよりも、集計表として公開しているため、すぐに統計処理をかけられない。
・データを提供する場所が、統一されておらず、わかりにくい。
　等々。

　行政側からするとデータを出すにしても色々なところとの調整やら、面倒なクレームやらが多いんで、なんとか出せるギリギリが上記の（使い勝手の悪い）データとなっていることは重々承知なんですが、出来れば将来的には API叩いたらJSONやらXMLやらで正規化されたデータがポンと返ってくるようになると良いなと思ってます。
　頑張って！！特に厚労省と総務省！！
　ちなみに、近年は統計法改正によって、オーダーメイドでのRaw データ切り出しに応じてくれるようになりましたね。すばらしいと思います。

#以上まだ前置き

本題：今回の目的は
　さて、今回はR の記事ですので、厚生労働省が公開しているExcel ベースでのデータをきれいにして統計処理をある程度しやすい形のデータに持って行くことを目的にR をいじってみました。

★使用するデータ

H22年度ＤＰＣデータ（分科会公開データ）

　ＤＰＣデータとは、全国の包括支払制度で入院料を算定している（一部は準備期間中の）病院から厚労省に提出される患者の情報です。H22時点で、1649病院（病院がH22年10月段階で8670施設なので、19％）がデータを厚労省に提出しており、そのデータを元に厚労省サイドで集計を行い、支払制度の調整や改正を実施しています*1。
　この際の集計結果は、中央社会保障審議会ＤＰＣ評価分科会*2にて報告され、そのデータはHPを介して誰でも確認する事が出来ます。

　こちらのデータの内
（７）疾患別・手術別集計
（９）疾患別・手術有無別・処置１有無別集計
（１０）疾患別・手術有無別・処置２有無別集計
　あたりを使うことで、全国のＤＰＣ病院の患者分布をDPCコードと呼ばれる疾患（と処置）を用いた分類毎に出す事が出来ます。特定地域で特定の疾患に強い病院は何処なのかといった情報を可視化可能なので、便利です。専門的に、そういった情報を配信しているサイト*3も存在したりします。

今回は、（７）疾患別・手術別集計を利用して、データの整形を実施します。

★使用するパッケージ

XLConnect

　dichikaさんの資料*4
　phosphor_mさんの資料*5
　などを参考にしています。

reshape2

　a_bickyさんの資料（reshape2に関連する部分）*6を参考にしています。

１．まずはダウンロード
　データを全てダウンロードします。（リストは、下記のURLに保存。）
https://docs.google.com/spreadsheet/pub?hl=en_US&hl=en_US&key=0AmNL5KN9Q1ZQdHdYb3NyeW52dWNWYWxMNElyZ2l4RlE&single=true&gid=0&output=csv

　分科会のページソースをgoogle refine使ってクレンジングしてURLを抽出、その後少し整形しています。リストのURLを適当なファイル名（リスト2列目）で保存します。取りあえず、R でExcelなどをバイナリ形式で落として、保存するソースは下記の通り。
　download.file()でmode="wb"にすると良いらしい。データを落としてくる間で、エラーが起きると危険かなぁ。RCUrlのgetURLオプションで、curlのretry件数指定が出来ると良いんだけど、やり方わからなかった。

url.list <- read.csv("DPC_URL_list.csv",head=TRUE)
dir.create("C:\\dpc_temp")
setwd("C:\\dpc_temp")
for(i in c(1:length(url.list[ , 1]))){
   url <- as.character(url.list[i, ]
   download.file(url, paste(,url.list[i,2],".xls",sep=""), mode="wb")
   }

２．さて、データを見てみると・・・・
　公開されているデータを見てみると、いくつか面倒な事象にぶつかります。
　-ヘッダ部分が結合セルになっている。
　　これが、意外と面倒でデータを簡単に整理して統計にかけられない要因の1つとなっています。また、これのせいかはわかりませんが、XLConnectで読み込みすると、列の右の方が切れることがあります。
　-ヘッダの行数がファイルによってまちまち。
　　これも、データを整理し難くする要因です。
　-何か、おしりに脚注が入っている。
　　これも。。。（略）

３．どうしようか。。
　　取りあえず、XLConnect 使ってデータをR に取り込みます。
　右の列が切れてしまう事象については、取りあえず256列目まで読んで回避します。真ん中あたりで、NAが続いたらcolの終わりと考えてそれ以降のデータを削除するというはなはだ危ういロジックを使っています。また、おしりの脚注を除くために読み込み行数を指定します。読み込み行数は先ほどのcsvの３列目に記載。ここら辺の引数はreadWorksheetないの endCol, endRowで指定できます。

for(i in 1:nrow(url.list)){
   wb <- loadWorkbook(paste(url.list[ i, 2],".xls",sep=""), create = TRUE)
   x <- readWorksheet(wb,sheet = 1,startCol = 1, startRow = 1, endCol = 256, endRow = url.list.shisetsu[ i, 3], head = FALSE)
   x <- x[ ,!is.na(x[100 ,])]
   eval(parse(text=paste(url.list.shisetsu[ i, 2],"<-x",sep="")))
}

４．どうしようか。。。２
　　さらに、問題なのはヘッダ中に登場する結合セルです。XLConnectでは、結合セルは左上端のセルに値を入れて、その他はNAで返す仕様となっています。今回のファイルは行をまたぐ結合はないので、header部分でNAを左端の値にそろえて値のreplaceをかけます。

#data <- 対象とするdata.frame
#header <- headerの行数
header.replace <- function(data, header){
	for(z in 1:header){
	x <- vector(mode = "character", length = length(data[z, ]))
	a <- (1:length(data[z, ]))[!is.na(data[z, ])]
if(length(a)>1){
		for(i in 1:(length(a)-1)){
	            b <- c(a[i], a[i+1]-1)
	            x[b[1]:b[2]] <- data[z, b[1]]
 	            }
	   b <- a[length(a)]
	   x[b:length(data[z, ])] <-  data[z, b[1]]
	   data[z, ] <- x
	}else{
	   b <- a[length(a)]
	   x[b:length(data[z, ])] <-  data[z, b[1]]
	   data[z, ] <- x
}
}
return(data)
}
#MDC_OP(手術分岐があるMDC06コード）を処理する。
data.list <- ls()
data.list <- data.list[grep("MDC_OP",data.list) ]
header.replace.4 <- function(data){header.replace(data,4)}
for(i in 1:length(data.list)){
	data <- eval(parse(text=paste(data.list[i])))
	eval(parse(text=paste(data.list[i]," <- header.replace.4(data)",sep="")))
}

５．どうしようか。。。３
　最後にヘッダ部分を１行にまとめましょう。

#change.name関数は、複数行で欠かれたHeaderを元に、カラム名を作成する。
#data <- 必要となるdata.frame
#header <- headerの数
change.name <- function(data, head){
	col.name.vec <- vector(mode="character", length(ncol(data)))
	for(i in 1:ncol(data)){
	col.name.vec[i] <- paste(data[1:head, i],collapse="_")
	}
	colnames(data) <- col.name.vec
	data <- data[-1:-head,]
	return(data)
}
for(i in 1:length(data.list)){
	data <- eval(parse(text=paste(data.list[i])))
	eval(parse(text=paste(data.list[i]," <- change.name.4(data)",sep="")))
}

６．ヘッダ部分を１行にしたデータを元にlong形式データを作成
　　reshape2は、データを縦持ちデータに変換するときもとても便利です。melt()関数で、IDとなる列を指定すれば、ぱぱっとlong形式データにしてくれます。ばらばらのデータフレームを一気に1つのデータフレームにまとめます。

data.OPE.full <- NA
for(data.name in data.MDC.OPE){
   data.a <- eval(parse(text = data.name ))
   data.OPE <- melt(data.a, id = c(1,2))
   data.OPE <- cbind(data.OPE ,colsplit(data.OPE[ , 3], "_" ,names = c("DPC06","DPCname","value_type","OPE","PRO")))
   data.OPE$source <- "OPE"
   data.OPE <- data.OPE[ , c(1,2,5:10,4)]
   data.OPE.full <- rbind(data.OPE.full,data.OPE)
}
#最後に少し、カラムを修正
data.OPE.full[ , 1] <- as.integer(data.OPE.full[ ,1 ])
data.OPE.full[ , 3] <- paste("0",data.OPE.full[ , 3], sep="")
data.OPE.full[ , 3] <- substring(data.OPE.full[ , 3],nchar(data.OPE.full[ , 3])-5,nchar(data.OPE.full[ , 3]))