【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]))

★データの解説
　出来たデータは、

• １列目　告示番号

　ＤＰＣ病院の告示番号です。

• ２列目　施設名

　ＤＰＣ病院の施設名称です。

• ３列目　DPC06

　ＤＰＣ上６桁です。疾患の分類を示します。

• ４列目　DPCname

　ＤＰＣ上６桁の名称です。大まかな疾患名を示します。

• ５列目　value_type

　数値の種類（在院日数か患者数か）を示します。

• ６列目　OPE

　手術あり症例か、なし症例かを示します。

• ７列目　PRO

　処置分岐の有無を示します。

• ８列目　source

　今回のデータは（７）疾患別・手術別集計から取得しましたが、
　他のソースから取得した場合、値が異なる（端数処理の関係）ため
　ソースがわかるようにしています。

• ９列目　value

　値です。

★これを使って
　先ほどＵＰした、ＤＰＣ病院一覧と紐づけてSPファイルにすることができます。
　結果として例えば、「関東近県の脳腫瘍：ＤＰＣ上六桁010010（手術あり）症例の分布は？」みたいな処理を、ある程度気軽に実施できます。

皆さんも、お時間があればやってみて下さい。
以上、興味ない人には一切興味ない誰徳？な記事でした。

追記：みなさんのブログを拝見すると、コードがきれいに書かれている訳ですが、あれはどうなっているんだろう？どうも、今ひとつな記載になってしまう。
2012-12-12-21:04
誰徳←誰得でした。
　
　

*1:参考：DPCはやわかりマニュアル　2010年度版（2010年7月発行）
http://di.mt-pharma.co.jp/medical/dpc_manual/

*2: 平成23年度第9回診療報酬調査専門組織・ＤＰＣ評価分科会
　http://www.mhlw.go.jp/stf/shingi/2r9852000001u23a.html

*3:病院情報局
http://hospia.jp/

*4:快適なエクセルライフのご提案
http://www.slideshare.net/dichika/tokyor18
[R]googleVisとXLConnectで医療費を可視化
http://d.hatena.ne.jp/dichika/20110617

*5:重回帰職人の朝は早い
http://www.slideshare.net/m884/ss-10334382

*6:RではじめるTwitter解析
http://www.slideshare.net/abicky/rtwitter