Rのコンソールフォント

デフォルトはCourieｒ　Newだけど何がよいのかなぁと考えていた。

間違えにくいのは、Consolas。O（オー）と0（ゼロ）、l（エル）と1の区別がつきやすい。ゼロの中に斜線が入っているから。

しばらくConsolasで使ってみよう。

version 3

とうとう R も version 3 になる。リリースは4月3日だそうだ。プレリリースは3月6日。その前に2月19日に、version 2.15.3 がプレリリース、3月1日にリリースされる。

覚えよう

ggplot2 の使い方を覚えよう。今までは普通のplotとlattice。

Rのグラフ描画って、latticeみたいな古いものはオリジナルを引きずっているけど、ggplot2になるともはや別言語なんじゃないの？というくらい違う。latticeで感覚的にやれたことがほとんど違う言語でやることになる。

でもその分、つくり込めるらしい。

私はXpose4を激しく使うので、latticeの勉強さえしていればよかったのだが、やはり新しいものへということでggplot2にもてを出すことにした。

ちなみに、大学のRのセミナーでの格言。

「どれか一つのパッケージを選んで深く理解して使うようにするのがベストです」

そのとおりだなぁ（笑）。これ、統計ソフトにも言える。RでもS-PLUSでもSASでもいい。使い込んだソフトを持てと。要はよく理解していないとデフォルトのオプションが何やってるかわからないで、解析しちゃまずいでしょということ。

私は事実上、統計はSASしか使えないのだ。Rは集計と描画だけ。データセット作成や集計についてのみ、2ソフトできるだけ。非線形混合効果についてだけNONMEMが加わってる。

あれれ。その禁をついにおかして、ggplot2だろうか。

動作しない原因

セミナー中、どうしてもxposeでstepwise GAMができなかった。動かなかった。

xpose4の他の機能は使えたのに。なぜそんなことが起こったんだろう。

コントロールファイルのほうのTABLEの指定がおかしいのかな、と思ったんだけどどこも間違っていない。

そして、家に帰ったら。動いた。。。 orz

家に帰って、ちょっと仕事をしようと思ったのだけど挫折。

今日のデータ解析の復習とかしたくなってしまう。

Rprofile.site

R を使うとき、基本的にXpose4を使わないということがない。
よって、起動時に自動的に読み込むようにしている。

その場合、Rのインストールフォルダの、/etc/Rprofile.site というファイルをテキストエディタで書き変える。

.First <- function(){
library(xpose4)
library(latticeExtra)
cat("\n\nWelcome to Pharmacometrics Research Group\n\n", date(), "\n\n")
}

上の.Firstというのが起動時実行関数である。
libraryとして、xpose4とlatticeExtraを読み込むように設定した。

cat関数でなんとなくメッセージが出るようにもしたが、これはご愛嬌である。

なお、パッケージ読み込みのときのメッセージが冗長に感じるときは、suppressMessages(library(xpose4))とすればよい。
また、パッケージが多いときは、そこだけソースファイルにしてもよい。
遊びだけど。

R にどっぷりつかる

Control File を読み込み、一部を(例えば$SIMのSEEDとか、Xpose関係のTable番号とか)を書き直して再出力するためのRプログラムを昨日、けっこうな時間をかけて考えた。

SAS ならできるのだが、今の環境ではRを用いないといけない。この手の便利機能は、実はNONMEMの解析支援ソフトには普通に実装されていて、何も自前のプログラムを作る必要なない。例えば、前回のRunにおける推定値を次のRunの初期値にするには、PsNのupdate_initでできる。また初期値の変更のみならずXposeのための番号つけかえはPiranaで簡単に行える。ただ、R をベースにNONMEMを繰り返し動かすことは今後もありえるので、作成した。

重要な関数は、readLines、strsplit、grep、unlist、paste、catだ。これらが使えればなんとかなる。

また、グラフ描画についても粘着質に勉強中である。というか必要に迫られて描きまくっている。Pharmacometricsにおいて使いやすいRのグラフ機能は、base、lattice、latticeExtra、ggplot2である。これらをかたっぱしから使う。

ggplot2は最近売出し中（？）のようで、書籍もでている。たしかにgeom、aesなどでシステマティックにグラフが描ける。PMx ではConditional plotが多いので、ggplot2も使える。が、グラフの雰囲気がなんとなくlatticeのほうがなじむ気がする。S-PLUS, Rの基本機能であるとともに、xposeで長年使われてきたからであろう。

Data manipulationにおいては、plyrに使いやすい関数が多い。が、だいたいのところはbaseの関数でなんとかなる。NONMEMのデータセットであっても数字をcharacterとして扱ったほうがいいケースがあり、なるほどなと思う。そのため、characterの扱いを復習。

R につかり過ぎて、ちょっとくらくらする。

R の Pharmacometrics関連のパッケージ

いくつある?

PK, PKfit, PKtools, PKreport, PKGraph, PKmodelfinder, PopPK.　このあたりはPK、PopPKのためのもの。

MIfuns, npde, Xpose4.　これらはわかりやすくPharmacometricsゆかりのもの。

drc, drfit のように用量反応性の解析のためのパッケージも?

nlmeODE のように事実上、Pharmacometricsのためといっていいものもある。Tornoeだし。

Xpose4 以外はほとんど使っていない。最近、人数も増えてきたし、このあたりの機能開拓をしてもいいのかもしれない。

AUCの計算

RでAUCの計算をするということは誰でもやっていることだろうが。。。

AUCcalc <- function(t,Cp){
tmp <- c()
for (i in 1:(length(t)-1)) {
tmp[i] <- (Cp[i]+Cp[i+1])*(t[i+1]-t[i])/2
}
return(sum(tmp))
}

こんなのでできる。AUC計算が教科書通りなので、わかりやすい。

しかし同じことをするにも。

AUCcalc <- function(t,Cp) sum(diff(t) * (Cp[-1] + Cp[-length(Cp)])) / 2

このほうが R らしいプログラミングだと思う。

他にROCパッケージにもAUCという関数があるらしいが、それは使ったことがない。

# AUClast(Linear Trapezoidal)計算
# 1名分のAUClast計算
AUCcalc <- function(x,y){
if(x[1]!=0){
x<-c(0,x)
y<-c(0,y)
}
sum(diff(x)*(y[-1]+y[-length(y)]))/2
}

# 各被験者のAUClastを計算し、resultへ格納する
RepAUC <-function(dat){
result <- c(length(unique(dat$ID)))
for(i in unique(dat$ID)){
tmp <- dat[dat$ID==i,]
result[i]<-AUCcalc(tmp$TIME,tmp$DV)
}
return(result)
}

WinNonLin Phoenix を使えという話もあるが、NONMEMでモデリングするときに各被験者のAUCを算出してグラフにしたかったりとかそういうことがあって。他のソフトから結果を持ってくるよりも、R/S-PLUS内で計算したかった。

という風に、自分にほしい道具をちまちま作りつつ作業をしていると、RにはPKというパッケージがあることを見つける。

このPKというパッケージは、non-compartmental analysis のためのパッケージのようだ。上に書いたような関数は、auc.complete(conc,time,group=...)としてできそう。

他にも半減期計算したりする関数がある。マニュアルは2011.9に改訂されているし、今のバージョンは2011年のものだから、定期的にメンテナンスされているのかもしれない。

半減期はプロットをみながら点を選択したりするが、WinNonLinのようにAdj-Rsqから自動選択することなぞできれば便利だろうな。

NONMEMのモデリングのときに計算する可能性があるのは、AUCとCmax。そのAUCがinfiniteである必要は必ずしもないから、AUClastがでれば十分だろう。

ループのProgressを外部ファイルへ

Progress カウンタを簡単につくってループに入れ込む。R内部で反復回数がよくわからないからだ。
if(irep==1){ PROGRESS[1,] <- c(i,Sys.time()) }
  else{PROGRESS <- rbind(PROGRESS,c(i,Sys.time()))}
write.table(file="Progress.tab",sep=",",row.names=F,PROGRESS)
既存のファイルに上書きするやりかたをしている。最後の行だけを追加するというやり方が思いつかなかった。というか、できるのだけどなんか効率が悪かった。

Simulationの律速段階は全然違うところにあるから、これでも速度には影響あるまい。

文字列を扱う

MEM in S などを勉強したりしているので、頭がRモードになりやすい。

さて、Rの場合、文字列を扱う関数も強力である。

cat(), deparse(), formatC(), grep(), match(), pmatch(), nchar(), parse(), paste(), strsplit(), sub(),gsub(), substring(), toupper(), tolower()

これだけの関数がリゲルの Programmieren mit R には載っている。知らないものもあって面白い。

NONMEMの出力を読んだり、コントロールファイルをRの中で書き換えたりすることは基本テクニックであるが、これら関数に習熟しておいたほうがよさそうだ。

一番の基本はreadLines(ファイル名)として、1行1行テキストファイルを読むことである。やってみるとよい。

そこからコントロールファイルの一部を、例えばシミュレーションのseedを書き換える場合、substr(ｵﾌﾞｼﾞｪｸﾄ, 1, 4)=="$SIM"のように検索すれば、書き換える位置がわかる。そしてpaste関数で書き換える文字列を構成してあげるといい。

また、scan関数というのもある。これは文字列のひとまとまりを単位としてひたすら読み込む。行列に格納してもいいが、こうやるときりが悪いことが多いので、ベクトルにしてしまうほかはない。

数値を扱っているときは、行数、列数に規則性があることが多いが、文字列をハンドリングすると不規則なところも多く、なんというか独特のコーディングになるようだ。

ちなみに、文字列として加工しきったあと、数値としてほしいところはas.numeric関数で数値にもどしてやるとよい。