Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 12-31-2011

大晦日だ。

そば殻?

もちろんそば殻なんか使われていない。

やぎ。あまりにも良いのでアイフォンの待ち受け画面になった。

世界地図を見ていたら、「択捉島・1945年からロシアによって占領されている。日本が領土権を主張」と書いてあった。いかにも。

中米系スーパーでみつけたジーザスろうそく。使い道はよくわからないけれど、とりあえずインテリアなんだと思う。アップで見ると リアル でかなり怖い。

韓国系スーパーでみつけた十七茶。日本には負けたくないんだろうなぁ。

平成171819202122232425262728293031323310111210111213141516171819202122232425262728293031日 土曜日|どうしても分類不可コメント(2)

Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 12-21-2011

ところで、最近の日本のニュースと言えばTPP。環太平洋・・・協定。日本にとってプラスかマイナスか?とにかく一大ニュースだ。でも太平洋を渡ってみると、アメリカでは誰もそんなこと話題にしていない。ま、普通に生活していて耳にするニュースではない。

つい先日のロイターアメリカ版のニュースに載っていたけれど、その見出しが The most important trade deal you've never heard of だ。つまり『誰も聞いたことがないけれど最重要の貿易の話題』 アメリカでは誰も聞いたことがないニュース呼ばわりされている。

ここまで前置き。

僕はアメリカに住んでいるのでアメリカ産のお米を食べている。こしひかり品種。来年で創業100年のとある貿易会社が販売している。他にもキッコーマングループの貿易会社の名前もよく目にする。どちらも日本の会社だ。

最近、とあるアジア系食料品店で少し変なお米を見た。初めて見た。happyrice1.jpg「こうふくまい」なら分からないでもないけどねぇ。どうやら韓国の会社らしい。すぐ隣に韓国語が書かれたお米も売っているんだけど、日本語で名前が書いてある方が売れるんだろうか?

日本人がいろいろと英語を使うのと同じ感覚なんだろうか?ちょっと違う気がする。

この会社のページに他のこうふくごめもあった。

平成171819202122232425262728293031323310111210111213141516171819202122232425262728293031日 水曜日|強いて言えばアメリカコメント(0)

Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 12-17-2011

Kaplan-Meier estimates をグラフにするプログラムを書いた。Rだ。

仕事上、何度となくグラフにしていて、その都度せくせくコードを書いていたのだけど、今回、昔書いたプログラムを大幅に改善した。生存曲線を書いて、下のマージンに numbers at risk を書く。

例。セットアップから。

library(survival)
kma <- survfit( Surv(time, status) ~ rx + adhere, data=colon )

この kmaの中にグラフにするべき数字が隠されている。

僕が書いた function はkmplot()という名前で ここ にあるので、

source('http://biostat.mc.vanderbilt.edu/wiki/pub/Main/TatsukiRcode/TatsukiRcodeKMplot.r')

で読み込める。

例の続き:

シンプルにkmplot(kma)で、それなりのグラフは得られるんだけど、細かにいろいろと設定することが可能。

kmplot(kma, mark='',
   xaxis.at=c(0,.5,1:9)*365, xaxis.lab=c(0,.5,1:9),
   lty.surv=c(1,2), lwd.surv=1, col.surv=c(1,1,2,2,4,4),
   col.ci=0,
   group.names=c('Obs ','Obs tumor adh','Lev','Lev tumor adh',
    'Lev+5FU ','Lev+5FU tumor adh'),
   group.order=c(5,3,1,6,4,2),
   extra.left.margin=6, label.n.at.risk=FALSE, draw.lines=TRUE,
   cex.axis=.8, xlab='Years', ylab='Survival Probability',
   grid=TRUE, lty.grid=1, lwd.grid=1, col.grid=grey(.9),
   legend=TRUE, loc.legend='bottomleft',
   cex.lab=.8, xaxs='r', bty='L', las=1, tcl=-.2
)

これで、この図 が得られる。

唯一絶対必要なインプットはsurvfit()で作ったオブジェクト(この例ではkma)だけ。それ以降のインプットは・・・

まずmarkは、打ち切りデータをどう表示するか。特に気にしなくて良い。いらない場合はmark=''

次にxaxis.at。これはn at riskを計算して表示する x軸の値(つまり時)を指定する。何もしていないと適当な所になる。xaxis.labxaxis.atで表示するラベル。例えば、時間は日にちで指定されているけれど、年を表示したいという時は、xaxis.at=(0:5)*365, xaxis.lab=0:5とすれば、0日、365日、730日のところに、0,1,2と表示することが出来る。

lty.surv, lwd.surv, col.survは、それぞれ線の種類、太さ、色だ。グループが4つあれば、4つ指定することもできるけれど、4つ未満の場合は繰り返して使われる。Rではいつものことだ。lty.surv=1:2, lty.lwd=1, lty.col=1:4とか。

lty.ci, lwd.ci, col.ciは confidence interval を表示する時に使う。表示したくなければ、lty.ci=0あるいはcol.ci=0とする。

group.namesは別に指定しなくていいんだけど、その場合、kmaオブジェクトで使われているものをそのまま引用する。それだと、たいていの場合うざい。
group.names=c('trt1','trt2','trt3','trt4')とかに書き換えてしまえば良い。 group.orderを指定しないと、ま、そのままグラフになるのだけど、n at riskの表示と、上のグラフがだいたい同じ順序だと見やすい。そこで、表示の順番を変えられる。例えば、group.order=c(2,4,1,3)とすると2番目のグループを1番上に書いて、その次が4番目のグループ、その下に、1番目のグループ、3番目のグループとなる。ま、例を見てください。

group.namesを短く指定してもそれでも見切れる場合はextra.left.marginを指定する。label.n.at.risk=TRUEとすると"Numbers at risk"と表示する。必要ないかな。

grid=TRUEで座標を表示する。細かな設定はcol.gridなどで。

legend=TRUEにすると、グループの名前と線の色、種類なんかを表示できる。場所はloc.legendで変えられる。

その他のインプットはplot()に行くので、cex.labelbtyなども指定できる。

図だけで、下の数字はいらないや、という場合はsimple=TRUEで消せる。

他にも ggplot を使って似たようなことをした人は何人かいるけれど、ggplot なんか使わないで済むならそれにこしたことはない、と思う。

平成171819202122232425262728293031323310111210111213141516171819202122232425262728293031日 土曜日|統計学コメント(2)

Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 12-11-2011

アイフォン4Sを買った。今までの使っていた携帯電話は、日本では10年くらい前に使っていたのと似た感じだったのでいきなりかなりのレベルアップ。テキストメッセージさえ送受信出来なかったんだから。写真は撮れるけれど撮った写真をどうしようもない。(ま、電話としてしか使わないプランだったからだけどね)

いきなりアイフォンだ。出先でメールチェックできるのは便利だね。ま、でもそんなに急ぎの用事なんかはめったにないのだけど、「あれ?ミーティングの場所どこだっけ?」ということがなくなる。理論上では。

ま、いいや。

そこで悲しいのがアイポッド。存在意義が無くなった。役立たず。無用。無駄。時代遅れ。全く存在意義が無いわけでもないか。今のところ思いつかないけれど。悲しさで言ったら(前にも使った比喩だけど)アニマルキングダムができたあとのジャングルクルーズだ。

平成171819202122232425262728293031323310111210111213141516171819202122232425262728293031日 日曜日|強いて言えばアメリカコメント(0)

Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 12-10-2011

10月にヴァンダービルト大学がとあるギネス記録を更新した。僕も貢献した。

「1日で何人にインフルエンザ予防接種を出来るでしょうか」記録。

記録は 12,850人。確か8時間だったか。今までの記録が 6,000人強だったので、一気に二倍。

つい先週、ヴァンダービルトの記録が正式にギネス記録として登録されたというニュースがあった。

平均で1分に30人弱に注射をしていたことになるので、かなりの規模だ。たしか常時40人が注射をしていた、とかそんなんだったと思う。

緊急時に大学関係者に速やかに予防接種をするプランというのがあるのだけど、それのテストを兼ねてのイベントだったらしい。

平成171819202122232425262728293031323310111210111213141516171819202122232425262728293031日 土曜日|仕事コメント(0)

Sunday, Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, 12-02-2011

例えばナッシュビル動物園に来る人は何人グループで来るんだろうか?という疑問を持ったとする。一人で来る人もいるだろうし、恋人同士もいるだろうし、家族4人とか、家族+親戚9人とか、学校の行事で120人ご一行様、とかもあるだろう。

ま、とりあえずグループの平均サイズを知りたいとしよう。理由はなんであれ、ね。

チケットを自動改札機に通す時に(本当はナッシュビル動物園にそんなハイカラなものはないのだけど)、ランダムにピピッとなるようにして、それがなった人に「あなたのグループは何人ですか?」と訊いてデータを集める。

ランダム(無作為)にサンプルを選んでいるので、このサンプルの平均を求めればきっとバイアス(偏り)の無い推定ができるに違いない。

違いない?

★☆☆★

例:とある小さな動物園に昨日やってきたグループは全てで10グループでした。小さな動物園だし、昨日は雨だったから10グループしか来なかったのだ。

グループのサイズは 1, 1, 2, 2, 2, 3, 4, 4, 4, 10。一人で来た写真家、雨の中のデート、雨天決行の強行ツアーグループ10人など。

平均グループサイズは (1+1+2+2+2+3+4+4+4+10)/10 = 3.3 で良い?

もし、平均グループサイズをこう定義すると、前述の「ランダムにピピッ」というサンプルのとり方では、使えるデータは取れない。というのも、グループの中で誰かひとりでもサンプルに入ればそのグループの人数がデータとして記録されるので、大きなグループほどサンプルに入りやすくなってしまうからだ。10人のグループはサンプルに入る確率が1人グループの10倍ある。

この弱小動物園の例で計算してみると、この日の入場者数は 33 人だ。その一人ひとりのグループサイズを考慮すると、グループサイズ1が2人、2が6人、3が3人、4が12人、10が10人ということになる。平均だと( 1 + 1 + 2+2 + 2+2 + 2+2 + 3+3+3 + 4+4+4+4 + 4+4+4+4 + 4+4+4+4 + 10+10+10+10+10+10+10+10+10+10 ) / 33 = 5.2。

平均グループサイズは 3.3? 5.2? どっち?

平成171819202122232425262728293031323310111210111213141516171819202122232425262728293031日 金曜日|統計学コメント(2)

ブログ検索

カレンダー

JanFebMar AprMay JunJulAugSepOctNovDec 2005200620072008200920102011
- - - - 1 2 3
4 5 6 7 8 9 10
11 12 13 14 15 16 17
18 19 20 21 22 23 24
25 26 27 28 29 30 31

最近のコメント

履歴

2017年05月 (1)

2017年04月 (1)

2017年03月 (1)

2017年02月 (1)

2017年01月 (1)

2016年12月 (1)

2016年11月 (1)

2016年09月 (1)

2016年08月 (1)

2016年07月 (1)

2016年05月 (1)

2016年04月 (1)

2016年03月 (1)

2016年02月 (1)

2016年01月 (1)

2015年12月 (1)

2015年07月 (1)

2015年06月 (2)

2015年05月 (2)

2015年04月 (3)

2015年03月 (2)

2015年02月 (1)

2015年01月 (3)

2014年12月 (1)

2014年11月 (1)

2014年10月 (2)

2014年09月 (2)

2014年08月 (2)

2014年07月 (1)

2014年06月 (1)

2014年05月 (2)

2014年04月 (1)

2014年03月 (2)

2014年02月 (1)

2014年01月 (2)

2013年12月 (1)

2013年11月 (3)

2013年10月 (2)

2013年09月 (1)

2013年08月 (4)

2013年07月 (1)

2013年06月 (2)

2013年05月 (2)

2013年04月 (3)

2013年03月 (1)

2013年02月 (1)

2013年01月 (4)

2012年12月 (1)

2012年11月 (3)

2012年10月 (1)

2012年09月 (1)

2012年08月 (3)

2012年07月 (3)

2012年06月 (2)

2012年05月 (6)

2012年04月 (2)

2012年03月 (8)

2012年02月 (2)

2012年01月 (1)

2011年12月 (6)

2011年11月 (5)

2011年10月 (4)

2011年09月 (6)

2011年08月 (9)

2011年07月 (5)

2011年06月 (5)

2011年05月 (5)

2011年04月 (6)

2011年03月 (17)

2011年02月 (6)

2011年01月 (10)

2010年12月 (10)

2010年11月 (4)

2010年10月 (6)

2010年09月 (5)

2010年08月 (11)

2010年07月 (8)

2010年06月 (8)

2010年05月 (3)

2010年04月 (8)

2010年03月 (11)

2010年02月 (4)

2010年01月 (8)

2009年12月 (6)

2009年11月 (6)

2009年10月 (6)

2009年09月 (7)

2009年08月 (6)

2009年07月 (10)

2009年06月 (10)

2009年05月 (10)

2009年04月 (6)

2009年03月 (7)

2009年02月 (9)

2009年01月 (12)

2008年12月 (6)

2008年11月 (10)

2008年10月 (8)

2008年09月 (9)

2008年08月 (12)

2008年07月 (8)

2008年06月 (12)

2008年05月 (12)

2008年04月 (12)

2008年03月 (11)

2008年02月 (10)

2008年01月 (10)

2007年12月 (12)

2007年11月 (14)

2007年10月 (13)

2007年09月 (11)

2007年08月 (16)

2007年07月 (10)

2007年06月 (10)

2007年05月 (6)

2007年04月 (10)

2007年03月 (13)

2007年02月 (10)

2007年01月 (8)

2006年12月 (13)

2006年11月 (15)

2006年10月 (9)

2006年09月 (8)

2006年08月 (18)

2006年07月 (14)

2006年06月 (16)

2006年05月 (23)

2006年04月 (20)

2006年03月 (12)

2006年02月 (14)

2006年01月 (21)

2005年12月 (20)

2005年11月 (17)

2005年10月 (18)

2005年09月 (16)

2005年08月 (10)

2005年07月 (6)