【第283回】 ベーコンライスの朝から（2013年3月1日）
たしか昨夜は妻からの依頼に応えてRで作図をしてメール送信した後，1:30過ぎに布団に潜り込んだけれども，ショパンのピアノコンチェルトを流していたら妙に寝付けなかったので，たぶん眠ったのは2:00過ぎだと思うが，よく覚えていない。しかし残酷にも朝はやってくる。6:15には，いつものようにquasimodeの「Slow Motion feat. 土岐麻子」のイントロが軽快なビートを刻むので，否応なしに起きてしまう。
朝食をちゃんと作る元気が無かったので，ちょうどいい具合に残っていた白いご飯の上に，暫く買うまいと思っていたのに業務スーパーで見かけて，日持ちがするから念のためと思って買ってしまったベーコンを3枚載せ，わさび漬けがないのでキムチを付け合わせて済ませた。それだけだとビタミンCが足りないので，はるみを1個食べた。うん，元気が出てきたぞ。
出発準備ができた時点で，ちょうど2番のバスが出てしまった頃だったので，3番までは30分待ちだ。というわけで，この時間を使って，厚生労働省の平成22年都道府県別生命表の概況の「図表データのダウンロード」から，昨夜ダウンロードしておいたExcelファイルを加工して，Rで扱えるようにしてみた。このデータを使って，都道府県別平均寿命の推移を示す折れ線グラフと，死因別損失余命の都道府県別プロファイルを示すレーダーチャートを，男女別に作ってみようというわけだ。今から出発までと，バスの中くらいで，完成させたい。できたら以下に図とコードとデータを載せる予定。松原みきベストをBGMにして作業開始。
では，まず折れ線グラフから。厚生労働省サイトからダウンロードしたExcelファイルを加工して作ったタブ区切りテキストデータはpref-e0-changes.txt（OSによっては都道府県名の2バイトコードを適当に変換しなくてはいけないかもしれない）。これをRの作業ディレクトリに置いて以下を実行すると，その下のグラフが画面に表示される（Windows環境で，直接pngファイルを出力させるコードは，change-e0.R）。
x <- read.delim("pref-e0-changes.txt")
males <- t(x[, 2:11])
colnames(males) <- x$PREF
females <- t(x[, 12:21])
colnames(females) <- x$PREF
COL <- ifelse(x$PREF=="長野", "blue", ifelse(x$PREF=="沖縄", "pink", "lightgrey"))
LWD <- ifelse(x$PREF=="長野", 2, ifelse(x$PREF=="沖縄", 2, 1))
LTY <- ifelse(x$PREF=="長野", 1, ifelse(x$PREF=="沖縄", 1, 3))
years <- 1965+0:9*5
layout(t(1:2))
matplot(years, males, type="l", col=COL, lwd=LWD, lty=LTY,
main="男性の都道府県別平均寿命の推移\n（青：長野，桃：沖縄，灰色：他）")
matplot(years, females, type="l", col=COL, lwd=LWD, lty=LTY,
main="女性の都道府県別平均寿命の推移\n（青：長野，桃：沖縄，灰色：他）")
名谷駅直通のバスは，時々1分前に行ってしまうことがあって，今日の3番がそれだったのだと思う。8:02に間に合ったのに，もう影も形もないのだった。仕方ないので，いったん帰宅した。プラゴミをちゃんと出していけるからいいか。
裏RjpWikiでご指摘あった件について補足。意図はもちろんわかるのだけれども，結果としてなされたことをみると，A，B，Cという3つの標本データだけがあって，AとB，AとCの間でそれぞれ平均値の差のt検定をして，それぞれの母集団の平均に差が無いかどうかを見ていることになり，Aが2回比較されているので，検定の多重性についての言及がなくてはいけないのではないかということ。多重性を考えなくてもいいようなデザイン，例えば，同じようなサンプリングの方法だとしても，AとB，CとDという比較であれば問題が無かったのだけれども，Aを2回比較してしまった以上は，AとBまたはAとCの比較の少なくとも一方において，本当は差が無いのに偶然有意水準より小さなp値が得られてしまう確率としての第一種の過誤（αエラー）を考えなくてはいけないのではないか？ というのが，あのコメントの主旨。そういう意味で，あのストーリー自体がわざと不適切に作られている気がしたので，最初は引っかけじゃないかと思ったのだ。
昨年末のRユーザ会で岡田さんがknitrとRStudioを駆使してプレゼンされたのを見て，いつか使ってやろうと思っているのだが，最近は世界中でこの組合せが人気らしく，パリの大学院生が共同で書いているブログでも取り上げられた。コードとプレゼンが一元管理できるのが魅力。
グラフ作成の続き。今度は，死因別のゼロ歳損失余命の都道府県プロファイルをレーダーチャートで示す（レーダーチャートを描画するには，0.3.8以降のfmsbライブラリが必要。まだミラーサイトでは0.3.7しか入っていないかもしれないので，ダメだったら数日待っていただきたい）。タブ区切りテキストデータは，pref-LLY-h22.txt（これもOSによっては都道府県名の2バイトコードを適当に変換しなくてはいけないかもしれない）。これをRの作業ディレクトリに置いて以下を実行すると，その下のグラフが画面に表示される（Windows環境で，直接pngファイルを出力させるコードは，LLY-profiles.R）。
x <- read.delim("./pref-LLY-h22.txt")
COL <- ifelse(x$PREF=="長野", "blue", ifelse(x$PREF=="沖縄", "pink", "lightgrey"))
LWD <- ifelse(x$PREF=="長野", 2, ifelse(x$PREF=="沖縄", 2, 1))
LTY <- ifelse(x$PREF=="長野", 1, ifelse(x$PREF=="沖縄", 1, 3))
VX <- c("悪性新生物","高血圧を除く\n心疾患","脳血管疾患","三大死因",
"肺炎","不慮の事故","交通事故\n（再掲）","自殺","腎不全","肝疾患",
"糖尿病","高血圧","結核")
males <- x[,2:14]
females <- x[,15:27]
require(fmsb)
layout(t(1:2))
radarchart(males, maxmin=FALSE, pcol=COL, axistype=2, pty=32, plty=LTY, plwd=LWD, vlabels=VX,
title="男性の死因別損失余命（平成22年度）\n（青：長野，桃：沖縄，灰：他都道府県）")
radarchart(females, maxmin=FALSE, pcol=COL, axistype=2, pty=32, plty=LTY, plwd=LWD, vlabels=VX,
title="女性の死因別損失余命（平成22年度）\n（青：長野，桃：沖縄，灰：他都道府県）")
意外に時間がかかってしまった。もう11:00を過ぎた。
このグラフはいろいろなことを示唆してくれる。一見してわかることは，平均寿命が男女とも最長の長野県は，男女とも，がんと肺炎による死亡が少ないということだ。一方，脳血管疾患によって失われている余命は比較的大きい。これは，長野県の人は漬け物をよく食べるため，元々塩分摂取量が多く，そのために脳卒中が多かったのを，食生活改善推進員さんが歩き回って塩分摂取量を減らし，そのおかげで脳卒中が減ったと言われているのだが，それでもまだ塩分摂取が高いということかもしれない。ただし，くも膜下出血のリスク因子としては遺伝も大きいので，塩分摂取だけが問題とは言い切れないが。なお，長野県では，男性のみ交通事故によって失われている余命が大きいが，これは子供の交通事故死だと思われる。細くて見通しが悪くて歩道が狭い道路が多いのに外遊びする子供は多いので，飛び出しによる交通事故が比較的多いのであろうことは想像に難くない。沖縄のプロファイルから目立つのは，肝疾患，糖尿病が高いことだ。たぶん飲酒が多いせいだろう。女性のみ結核による損失余命が大きかったが，これは流行があったのかもしれない。たぶん，ここで書いた推測の裏付けを取るようなデータを探して分析すれば，地域相関研究だけれども，卒論か学会発表ネタくらいにはなるかもしれない。
その後，メール送受信とか，昨日やっていた仕事の続きをしていたら，昼前に留学生のアリムさんが来室し，若干の書類作成。昼飯はどうしようかなあ。今日も自販機のパンか？
結局，今日もカップ麺にしてしまった。ダメダメ。自販機に歩く途中，あまりに戸外の風が強いのに驚いた。この辺も春一番なのか？ バスが止まらなければいいのだが。
昨日から方々で話題になっているWHOのレポートの原文はhttp://www.who.int/ionizing_radiation/pub_meet/fukushima_report/en/index.htmlからpdfファイルとしてダウンロードできる。これとは別にニュースリリースも出ているが，たぶんリリースやExecutive Summaryだけでなく，本文を読んだ方がいいと思う。172ページもあって，Lifetime Attributable Riskの計算など，考え方も簡単ではないが。
特別企画：インターネット世論調査はどうあるべきか？（津田大介の「メディアの現場」Vol.66 より）で，RjpWikiと，Rを使ったゼゼヒヒのコメントの分析が紹介されていたので，RjpWikiのトップページに書き込んできた。いい記事なんだが，「しかし現在は、R Commander [*22] やR Studio、Rz [*23] 、R Excel [*24] などといった、Rの利用環境を向上させるツールも開発されるようになっている」という部分，R StudioではなくてRStudioなので念のため。あと，ここまで書くなら，EZRにも触れて欲しかったところ。なお，この記事では岡田昌史さんの肩書きが東工大の准教授とされていて吃驚したのだが，Googleで検索してみたら，同姓同名の別人であった。RjpWikiの岡田さんの所属は違うので，お間違いなきよう（注：最初，ここには正しい所属を書いていたが，岡田さんから津田さん宛に，RjpWikiは個人でやっているものだから所属は消して欲しい旨tweetがあったので，ここからも所属は消した）。
研究業績をずっと放置していたのだが，国際協力研究科から，一覧表をwebから作っていいか問い合わせるメールが来たので，久々に更新した。ただし学会の座長とかは書いていないんだが。それくらいはいいか。
R-2.15.3がリリースされた（R-announceでPeter Dalgaardからのメールも読めるが，このR-bloggersの記事の方が見やすいので，こっちをリンクしておく）。4月になるとR-3.0.0がリリースされることになっているので，これが最後のバージョン2系のリリースとなる。細かい改良がほとんどのようだ。例えば，最初の改良点"lgamma(x) for very small x (in the denormalized range) is no longer Inf with a warning."は，これまでlgamma(1e-307)は706.8936を返すけれどもlgamma(1e-308)が警告メッセージとともにInfを返していたのが，1e-308よりゼロに近い値でも数値を返すようになったということか。Windowsユーザにとっては，Windows環境でも十分にRAMを積んでいれば，64ビット版ではMacOSやLinuxと同様，R_allocが32GBを確保できるようになったことは大きいかもしれない。バグフィックスでqt(1e-12,1.2)がNaNを返していたのを修正したとあるが，自由度1.2のt分布ってどういう場合に意味があるのか想像できないくらいレアケースだと思う。良く見つけたな，というレベルのバグ。
歯科受診のため19:44に名谷駅を出るバスで帰途に就いた。治療して貰ってから帰宅し，ベーコン野菜炒めを作って晩飯。バンコクで着ていた半袖の服を含め，洗濯物の蓄積が限界に達したので洗濯して部屋干し。
0:00からの宮嶋みぎわさんのUstライブを見ようと思ってブラウザでそのページを開いていたのだが，いつの間にか寝落ちしていて見逃してしまったのが悔やまれる。
付記：このメモからは，都道府県別生命表データ関連の記述を，「都道府県別生命表掲載データからの作図法(How to make charts by prefecture data)」として独立させた。
△Read/Write COMMENTS
▼前【282】（2月末日にradarchart()関数に機能追加（2013年2月28日） ） ▲次【284】（寒くて目が覚めた朝には二度寝を（2013年3月2日） ） ●Top