演習の本をRでこなしていくことにした。

ファーマコキネティクス―演習による理解
杉山 雄一,山下 伸二,加藤 基浩
南山堂 / ￥ 6,300 ()
在庫あり。

設問1

β層の消失速度をもとめて、それを用いてα層の消失速度を決定する。

t <- c(1,2,3,6,13,23,30,39,49)
c <- c(68.2, 51.1, 43.2, 34.0, 18.3, 11.5, 10.7, 10.6, 6.78)

23時間からの濃度からβ層の消失速度を求める

logc <- log(c)
result <- lm(logc[6:9] ~ t[6:9])
coef(result)
 (Intercept)      t[6:9] 
 2.93770720 -0.01888941 

β層の消失速度が決まったので、次は23時間までの濃度からβ層の寄与を引いた濃度からα層の消失速度を決定する

beta <- function(t){exp(2.937707-0.01889*t)}
alogc <- c[1:5] - beta(t)[1:5]
aresult <- lm(alogc ~ t[1:5])

これでα層β層それぞれの消失速度がもとまる

result

Call:
lm(formula = logc[6:9] ~ t[6:9])

Coefficients:
(Intercept)       t[6:9]  
    2.93771     -0.01889  

aresult

Call:
lm(formula = alogc ~ t[1:5])

Coefficients:
(Intercept)       t[1:5]  
     3.9916      -0.2088  

求めるモデルは

exp(3.9916-0.2088*t)+exp(2.937707-0.01889*t)

設問2

ラフに求めたA,B,alpha,betaを初期値として最小自乗法でフィッティングを行う。1exp,2expでやってみて、どっちのあてはまりがよいか考える

t <- c(1,2,3,6,13,23,30,39,49)
conc <- c(68.2, 51.1, 43.2, 34.0, 18.3, 11.5, 10.7, 10.6, 6.78)
dat <- data.frame(conc,t)
result <- nls(conc ~ a*exp(-alpha*t),start=list(a=50,alpha=0.045))
summary(result)

Formula: conc ~ a * exp(-alpha * t)

Parameters:
      Estimate Std. Error t value Pr(>|t|)    
a     63.26990    5.39879  11.719 7.45e-06 ***
alpha  0.08359    0.01777   4.706  0.00219 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 6.608 on 7 degrees of freedom

Number of iterations to convergence: 14 
Achieved convergence tolerance: 7.74e-06 

f <- function(t){63.2699*exp(-0.0836*t)}
plot(f,0,50,xlim=c(0,50),ylim=c(0,70),xlab="Time",ylab="Conc")
par(new=T)
plot(t,conc,xlim=c(0,50),ylim=c(0,70),axes=F,xlab="",ylab="")

一層の消失モデルだと特に遅い時間の当てはめが悪い。ちなみにRだとAIC関数使えば適合度がでる。

t <- c(1,2,3,6,13,23,30,39,49)
conc <- c(68.2, 51.1, 43.2, 34.0, 18.3, 11.5, 10.7, 10.6, 6.78)
dat <- data.frame(conc,t)
result <- nls(conc ~ a*exp(-alpha*t)+b*exp(-beta*t),start=list(a=50,alpha=0.277,b=20,beta=0.0198))
summary(result)
f <- function (t){55.75*exp(-0.3769*t)+29.157*exp(-0.03189*t)}
plot(f,0,50,xlim=c(0,50),ylim=c(0,70),xlab="Time",ylab="Conc")
par(new=T)
plot(t,conc,xlim=c(0,50),ylim=c(0,70),axes=F,xlab="",ylab="")

プロットするときにpar関数で上書きするのと、xlim,ylimを明示して軸をきちんとそろえつつ、目盛やタイトルを上書きして重ねないように気をつける。

ivivcというパッケージを発見。PKfitでkelとかVdを求めておけばこれで対話的に色々（というかルーチーンワークが）できるらしい。

demoのflashも用意されている。

****************************************************************************
*                                                                          *
*                                  ivivc                                   *
*                                 v 0.1.4                                  *
*                                                                          *
* In vitro-in vivo correlation (IVIVC) is defined as the correlation       *
* between in vitro drug dissolution and in vivo drug absorption.           *
* This package is used to develop and validate an IVIVC model.             *
* The following steps wii be conducted as follows:                         *
* ->1: Input/Edit In Vivo Absorption Data: IV, Oral solution or IR drug    *
* ->2: Develop an IVIVC Model: Fitting IV, Oral solution or IR drug        *
* ->3: Input/Edit In Vitro Dissoution Data and                             *
*      In Vivo absorption Data: ER drug with Different Release Rates       *
* ->4: Develop an IVIVC Model: Model Dependent Method                      *
* ->5: Evaluate an IVIVC model: Prediction Error                           *
*                                                                          *
*           Please type 'run()' to get started.                            *
*                                                                          *
****************************************************************************

runしてdemoするとなんか一瞬で終わる。

Rでここまでのワークフローを対話的にやらすパッケージってのは自由度が少ない気もするが、あとでまた。

とりあえずPKfitをいじろうと思ってるんだけど、pdfのマニュアルがいい加減すぎるんだよね。

chemoinformatics的には最もlateなフェーズである、QSPkRとかモデル動物からヒト外挿あたりにブームがやってきたので、色々と調べたりいじったりする日々を。

どっちかっていうとPK/PD,PBPKのようなモデルに興味があるのでそういったパッケージがないかと、CRANを漁った。

CRAN Task View: Analysis of Pharmacokinetic DataによるとPK/PD modelはnlmeODEパッケージを使えばいいらしい。

PBPKに関してはR-ForgeにRDynamicってのを見つけた。

ちょっとnlmeODEのサンプルを

data(Theoph)
TheophODE <- Theoph
TheophODE$Dose[TheophODE$Time!=0] <- 0
TheophODE$Cmt <- rep(1,dim(TheophODE)[1])
OneComp <- list(DiffEq=list(
dy1dt = ~ -ka*y1 ,dy2dt = ~ ka*y1-ke*y2),
ObsEq=list(
    c1 = ~ 0,
    c2 = ~ y2/CL*ke),
    Parms=c("ka","ke","CL"),
    States=c("y1","y2"),
    Init=list(0,0))
    TheophModel <- nlmeODE(OneComp,TheophODE)

Theoph.nlme <- nlme(conc ~ TheophModel(ka,ke,CL,Time,Subject),
data = TheophODE, fixed=ka+ke+CL~1, random = pdDiag(ka+CL~1),
start=c(ka=0.5,ke=-2.5,CL=-3.2),
control=list(returnObject=TRUE,msVerbose=TRUE),
verbose=TRUE)

plot(augPred(Theoph.nlme,level=0:1))

よくよく考えてみるに、既に3年くらい前にはこういうことに必要性を認識していたのだけどなぁ。

Computer Applications in Pharmaceutical Research and Development (Wiley Series in Drug Discovery and Development)

Wiley-Interscience / 17107円 ( 2006-06-30 )

シリーズ物だったのね、知らんかった。

Analysis of Integrated and Cointegrated Time Series with R (Use R!)
Bernhard Pfaff
Springer New York / 7631円 ( 2008-08-11 )

Data Manipulation with R (Use R!)
Phil Spector
Springer / 7861円 ( 2008-03-19 )

Bayesian Computation with R (Use R)
Jim Albert
Springer-Verlag / ?円 ( 2007-07 )

Applied Spatial Data Analysis with R (Use R!)
Roger S. Bivand
Springer / 8586円 ( 2008-08-14 )

Analysis of Integrated And Cointegrated Time Series With R (Springer Series in Statistics)
Bernhard Pfaff
Springer-Verlag / ?円 ( 2005-12-22 )

Applied Statistical Genetics with R: For Population-based Association Studies (Use R!)
Andrea S. Foulkes
Springer New York / 5834円 ( 2009-04-17 )

Analysis of Phylogenetics And Evolution With R (Use R)
Emmanuel Paradis
Springer New York / 6473円 ( 2006-07-28 )

Nonlinear Regression with R (Use R!)
Christian Ritz
Springer / 7631円 ( 2008-11-21 )

Wavelet Methods in Statistics with R (Use R!)
Guy Nason
Springer / 8089円 ( 2008-08-11 )

soaplibにはCherryPyが組み込まれているので簡単にウェブサービスがつくれる。さらにrpy2を使えばRのウェブサービスがお手軽に構築できて便利。

インストールの手順はwikiを参照

以下のコードをsoapr.pyというファイル名で保存する。これは与えた数字の数だけRがrnormして返すという簡単なサービス。

from soaplib.wsgi_soap import SimpleWSGISoapApp
from soaplib.service import soapmethod
from soaplib.serializers.primitive import String, Integer, Array, Float
import rpy2.robjects as robjects

class RService(SimpleWSGISoapApp):

    @soapmethod(Integer,_returns=Array(Float))
    def rnorm(self,num):
        r = robjects.r        
        y = r.rnorm(num)
        return [i for i in y]

if __name__=='__main__':
    from cherrypy.wsgiserver import CherryPyWSGIServer
    server = CherryPyWSGIServer(('localhost',7789),RService())
    server.start()

実行するとサーバが立ち上がりhttp://localhost:7789/soapr.wsdlにアクセスするとwsdlを出力する。

テスト用のクライアント

>>> from soaplib.client import make_service_client
>>> from soapr import RService
>>> client = make_service_client('http://localhost:7789/',RService())
>>> print client.rnorm(5)
[-0.40538156526399999, 0.45889158855099998, -0.67243766066699995, 2.6881017757299999, 0.113257266309]
>>> print client.rnorm(3)
[-0.16610789340500001, 0.66591497139099998, -1.74812066455]

参考

• Documentation for rpy2
• soaplib

少し読んでみたけど、グラフ分割問題の評価関数を固有値問題に落とすことで最適なクラスタを見つけるのは最初読んで大雑把につかんでたんだけど、あんまネットやらないのにWeb系の研究やってるっていうwwwのエントリとそこから辿れる論文でもうちょっと深く理解した。

ある種の化合物ライブラリのクラスタリング問題にも使えそうな感じなので、agglomerative clusteringでないやり方も覚えておくとよいですな。

以下写経。

>source("sc.R")
>x <- iris[1:4]
>d <- as.matrix(dist(x))
>d[d == 0] <- 0.1
>w <- 1/d
>ans <- mcut(w,3)
> ans[1:50]
 [1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[39] 2 2 2 2 2 2 2 2 2 2 2 2
> ans[51:100]
 [1] 1 3 1 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 1 1 3 3 3 3 3 1 3 3 1 3
[39] 3 3 3 3 3 3 3 3 3 3 3 3
> ans[101:150]
 [1] 1 1 1 1 1 1 3 1 1 1 1 1 1 1 1 1 1 1 1 3 1 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[39] 3 1 1 1 1 1 1 1 1 1 1 1

この本は集合知プログラミング読んで、クラスタリングに興味が湧いたら読んでみるといい本だと思う。

Rで学ぶクラスタ解析
新納 浩幸
オーム社 / 3360円 ( 2007-11 )

集合知プログラミング
Toby Segaran
オライリージャパン / 3570円 ( 2008-07-25 )

R newsにRandomForestの解説が載っていると教えてもらったので、早速読んだ。

Classification and Regression by randomForest

Some notes for practical useは一回読んでおくと吉。

昨日気になったので、とりあえずどんな感じか確かめようと。

> library("fastICA")
> library("mclust")
> data("iris")
> b <- fastICA(iris[1:4],2)
> mc <- Mclust(b$S,G=3,modelNames="VVV")
> plot(mc,b$S)

混合分布はクラスタリングの本にも載っているが、PRMLの下巻のほうがK-meansからEM,混合モデルの導入の流れがスムーズで分かりやすかったかな。

Rで学ぶクラスタ解析
新納 浩幸
オーム社 / 3360円 ( 2007-11 )

パターン認識と機械学習 下 - ベイズ理論による統計的予測
C. M. ビショップ
シュプリンガー・ジャパン株式会社 / ?円 ( 2008-07-11 )

出典がちょっとどこだかわすれたが。

• PCA 分散の大きな成分を抽出
• ICA 非正規性を最大にする成分を抽出

独立って、直交しないといけないんじゃないかなんて思っていたがこの資料みて納得。

その上で主成分分析、独立成分分析を読んだらさくっと理解できた。

> data("iris")
> a <- prcomp(iris[1:4])
>biplot()

> library("fastICA")
> data("iris")
> b <- fastICA(iris[1:4],2)
> plot(b$S)

ある種のデータを二次元のプロットとして見たい場合、PCAだと変に鎖上になってしまうものもICAだともうちょっと広がって見えてくれるかな。単にマップして見たいってだけだったらこういうやり方でもいいか。

あとPCAとかICAで次元を圧縮してから混合ガウス分布みたいなのでクラスタリングするってのはダメなのだろうか?

Rで学ぶクラスタ解析に書いてあったので、これを見るべし的。

Rで学ぶクラスタ解析
新納 浩幸
オーム社 / 3360円 ( 2007-11 )

8章が次元縮約の章で特異値分解,pLSI,NMFのコードが載ってる。

やはり、この本面白いですねぇ。