Drkcore

20092009 R MCMC JAGS

JAGSでメタアナリシス

マルコフ連鎖モンテカルロ法の5-2をJAGSで

マルコフ連鎖モンテカルロ法 (統計ライブラリー)
豊田秀樹
朝倉書店 / ￥ 4,410 ()
在庫あり。

モデル

model
{
  for( i in 1 : Num ) {
    rc[i] ~ dbin(pc[i], nc[i]);
    rt[i] ~ dbin(pt[i], nt[i]);
    logit(pc[i]) <- mu[i];
    logit(pt[i]) <- mu[i] + delta[i];
    mu[i] ~ dnorm(0.0,1.0E-5);
    delta[i] ~ dnorm(d, tau);
  }
  d ~ dnorm(0.0,1.0E-6);
  tau ~ dgamma(0.001,0.001);
  delta.new ~ dnorm(d, tau);
  sigma <- 1 / sqrt(tau);
}

Rのコード

library(rjags)

data <- list(rt = c(49, 44, 102, 32, 85, 246, 1570),
     nt = c(615, 758, 832, 317, 810, 2267, 8587),
     rc = c(67, 64, 126, 38, 52, 219, 1720),
     nc = c(624, 771, 850, 309, 406, 2257, 8600),
    Num = 7)

inits <- list(d = 0, delta.new = 0, tau=1,
     mu = c(0, 0, 0, 0, 0, 0, 0),
  delta = c(0, 0, 0, 0, 0, 0, 0))

m <- jags.model(
  file = "model.txt",
  data = data,  
  inits = list(inits),
  n.chains = 1
)


update(m,5000)

x <- coda.samples(
    m,
    c("d","sigma"),
    thin = 1, n.iter = 15000
)

summary(x)

plot(x)

結果

> summary(x)

Iterations = 6001:21000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 15000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

         Mean      SD  Naive SE Time-series SE
d     -0.1400 0.08028 0.0006555       0.001885
sigma  0.1285 0.08188 0.0006686       0.002196

2. Quantiles for each variable:

          2.5%      25%     50%     75%    97.5%
d     -0.31570 -0.18244 -0.1338 -0.0908 0.003366
sigma  0.02855  0.06916  0.1110  0.1672 0.331475

mcmc_5_2

19092009 R MCMC JAGS

JAGSでロジスティック回帰モデル

macでwinbugsを頑張って動かすのもアレなので、潔くJAGSで。

マルコフ連鎖モンテカルロ法の5-1

マルコフ連鎖モンテカルロ法 (統計ライブラリー)
豊田秀樹
朝倉書店 / ￥ 4,410 ()
在庫あり。

生態学のデータ解析 - R で JAGSによるとjagsの場合にはBUGS言語の文（式）の終わりにはセミコロンが必要。

model
{
  for(i in 1:N){
    r[i] ~ dbin(p[i],n[i]);
    e[i] ~ dnorm(0.0,tau);
    logit(p[i]) <- beta0+beta1*x1[i]+beta21*x21[i]+beta22*x22[i]
    +beta23*x23[i]+beta24*x24[i]+beta25*x25[i]
    +beta26*x26[i]+beta27*x27[i]+e[i];
  }
  beta0 ~ dnorm(0.0,1.0E-6);
  beta1 ~ dnorm(0.0,1.0E-6);
  beta21 ~ dnorm(0.0,1.0E-6);
  beta22 ~ dnorm(0.0,1.0E-6);
  beta23 ~ dnorm(0.0,1.0E-6);
  beta24 <- 0;
  beta25 ~ dnorm(0.0,1.0E-6);
  beta26 ~ dnorm(0.0,1.0E-6);
  beta27 ~ dnorm(0.0,1.0E-6);
  tau ~ dgamma(0.001,0.001);
  sigma<-1/sqrt(tau);
}

Rのコード

library(rjags)

data <- list(
r = c(1,13,44,155,92,100,18,0,0,2,1,0,2,0),
n = c(13,70,115,301,153,141,26,6,16,25,32,8,9,4),
x1 = c(0,0,0,0,0,0,0,1,1,1,1,1,1,1),
x21 = c(1,0,0,0,0,0,0,1,0,0,0,0,0,0),
x22 = c(0,1,0,0,0,0,0,0,1,0,0,0,0,0),
x23 = c(0,0,1,0,0,0,0,0,0,1,0,0,0,0),
x24 = c(0,0,0,1,0,0,0,0,0,0,1,0,0,0),
x25 = c(0,0,0,0,1,0,0,0,0,0,0,1,0,0),
x26 = c(0,0,0,0,0,1,0,0,0,0,0,0,1,0),
x27 = c(0,0,0,0,0,0,1,0,0,0,0,0,0,1),
N = 14)

inits <- list(beta0 = 0, beta1 = 0, beta21 = 0,  beta22 = 0, beta23 = 0, beta24 = 0, beta25 = 0, beta26 = 0,  beta27 = 0, tau = 1)

m <- jags.model(
  file = "model.txt",
  data = data,  
  inits = list(inits),
  n.chains = 1
)

update(m,10000)

x <- coda.samples(
    m,
    c("beta0","beta1","beta21","beta22","beta23","beta24","beta25","beta26","beta27", "sigma"),
    thin = 1, n.iter = 40000
)

summary(x)

plot(x)

サマリーの表示

> summary(x)

Iterations = 11001:51000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 40000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

           Mean     SD Naive SE Time-series SE
beta0  -0.02576 0.4188 0.002094        0.03425
beta1  -3.11738 0.6059 0.003030        0.01815
beta21 -3.03168 1.4556 0.007278        0.03269
beta22 -1.52877 0.6757 0.003379        0.02993
beta23 -0.28549 0.6355 0.003177        0.04374
beta24  0.00000 0.0000 0.000000        0.00000
beta25  0.35897 0.5926 0.002963        0.02957
beta26  1.05996 0.6389 0.003194        0.04301
beta27  0.70287 0.7677 0.003838        0.03026
sigma   0.34110 0.3980 0.001990        0.03152

2. Quantiles for each variable:

           2.5%      25%      50%      75%   97.5%
beta0  -1.11529 -0.15412  0.02137  0.16865  0.7406
beta1  -4.48731 -3.46881 -3.06188 -2.70570 -2.0882
beta21 -6.39669 -3.83976 -2.86884 -2.05688 -0.6382
beta22 -2.90198 -1.87085 -1.54216 -1.21238 -0.0738
beta23 -1.27301 -0.63488 -0.39022 -0.08097  1.4261
beta24  0.00000  0.00000  0.00000  0.00000  0.0000
beta25 -0.87361  0.10421  0.34604  0.59078  1.7201
beta26  0.08425  0.72545  0.96155  1.25666  2.7809
beta27 -0.80720  0.33750  0.71085  1.10396  2.1284
sigma   0.02816  0.08266  0.19387  0.45606  1.3892

5章は事例がたくさん載っているので、この勢いで一通りやっておこう。

5.1　ロジスティック回帰モデル　
5.2　メタ分析　
5.3　多項ロジットモデル　
5.4　対数線形モデル　
5.5　ポアソン回帰　
5.6　2値データに対する回帰分析　
5.7　トービット回帰モデル　
5.8　変曲点のある回帰分析　
5.9　生存時間分析(ワイブル回帰)　
5.10　生存時間分析(コックス回帰)　
5.11　時系列モデル　
5.12　分散分析　
5.13　分散成分分析　
5.14　分散分析(枝分かれ配置)　
5.15　一般化可能性理論　
5.16　反復測定データの分散分析　
5.17　階層線形モデル　
5.18　項目反応理論(2母数2値モデル)　
5.19　項目反応理論(段階反応モデル)　
5.20　項目反応理論(名義反応モデル)　
5.21　項目反応理論(部分採点・評定尺度モデル)
5.22　項目反応理論(連続反応モデル)
5.23　多次元IRT
5.24　項目反応理論(混合名義反応モデル)　
5.25　項目反応理論における特異項目機能(DIF) の分析　
5.26　正規混合モデル　
5.27　潜在クラス分析　
5.28　成長曲線モデル　
5.29　非線形成長曲線モデル　
5.30　因子分析　
5.31　多母集団分析　
5.32　非線形SEM

02092009 chemoinformatics R

openbabelのフィンガープリントを使ってSVMの予測モデルを作る

久々にケモインフォクックブックを更新した。

で、このフィンガープリント使って予測モデルを作ってみる。

Benchmark Data Set for in Silico Prediction of Ames Mutagenicity のSupporting InformationにAmes試験のデータセットがあるのでこれを使ってcsvを用意する。

babel -ismi ci900161g_si_001/smiles_cas_N6512.smi -ofpt ci900161g.fpt -xh -xfFP2

これでfptファイルを用意してクックブックので変換。

f2b.py ci900161g.fpt > fingerprint.txt

さらにフィンガープリントをcsvにするのは下のやっつけスクリプトで。

file = "ci900161g_si_001/smiles_cas_N6512.smi"

dic = {"0":"negative","1":"positive"}
act = {}
for l in open(file,"r"):
   smi,id,num = l.split()
   act[id] = dic[num]


fingerprints_file = "fingerprint.txt"

header = "ID,"

for i in range(1,1025):
   header += "bit" + str(i) + ","

header += "act"
print header

for l in open(fingerprints_file, "r"):
   col = ""
   id, fp = l.split()
   col += "\"" +id + "\","
   for c in fp:
       col += c + ","
   col += "\"" + act[id] + "\""
   print col

実行

python fconv.py > ci900161g.csv

でRで解析してみる。

ames <- read.csv("/Users/kzfm/ci900161g.csv")
set.seed(50)
tr.num <- sample(6512,2500)
ames.train <- ames[tr.num,]
ames.test <- ames[-tr.num,]
ames.svm <- ksvm(act ~.,data=ames.train[,-1])
ames.pre <- predict(ames.svm, ames.test[,c(-1,-1026)])
ames.tab <- table(ames.test[,1026],ames.pre)
sum(diag(ames.tab))/sum(ames.tab)

結果

>     sum(diag(ames.tab))/sum(ames.tab)
[1] 0.7771685
> ames.tab
          ames.pre
           negative positive
  negative     1424      425
  positive      469     1694

まぁまあかなぁ。フィンガープリントとかカーネルを考えればもうちょっと精度が上がる気もするが。

参考

[連載]フリーソフトによるデータ解析・マイニング第 31 回. Rとカーネル法・サポートベクターマシン.

マシンラーニング (Rで学ぶデータサイエンス 6)
辻谷将明,竹澤邦夫
共立出版 / ￥ 3,675 ()
通常1～2か月以内に発送

18082009 R PRML MCMC

SIR(sampling-importance-resampling)

水産資源学におけるベイズ統計の応用ワークショップのスライドの22,23枚目

g(x)からM個のサンプルx1*, …, xM*を発生
w*= f(x*)/g(x*)を計算
x1*, …, xM*からw1*, …, wM*に比例する確率でサンプルをm個とってくる
得られたm個のサンプルはf(x)からの事後サンプル

となっているのだが、23枚目のRのコードは

n2 <-sample(length(x1), 1000, replace=T, prob=dlnorm(x1,1.1,0.6))

とwを求めずにいきなりサンプリングしているので？となったが、一様分布だから重みを計算する必要ないのね。

わざわざ書くならこうか。

x1 <-runif(100000,0,60)
w <- dlnorm(x1,1.1,0.6) / (dunif(x1,0,60))
n2 <-sample(length(x1), 1000, replace=T, prob=w)
x2 <-x1[n2]
plot(density(x2))

パターン認識と機械学習下 - ベイズ理論による統計的予測
C. M. ビショップ
シュプリンガー・ジャパン株式会社 / ￥ 8,190 ()
在庫あり。

13082009 chemoinformatics R PRML

カーネル主成分分析

普通のPCAとちがうのは非線形主成分を扱えること。でも、カーネル関数をうまく選ばないといけないことを意味する。

Rとカーネル法を参考に。

> library(kernlab)
> x <- as.matrix(iris[,1:4])
> iris.kpc1 <- kpca(x,kernel="rbfdot",feature=2,kpar=list(sigma=0.1))
> plot(pcv(iris.kpc1),col=as.integer(iris[,5]))

kpca

あータニモトカーネルを使ったkernel PCAとかなんかいいかもしれない。

パターン認識と機械学習下 - ベイズ理論による統計的予測
C. M. ビショップ
シュプリンガー・ジャパン株式会社 / ￥ 8,190 ()
在庫あり。

06082009 R PRML

Rでベータ分布

PRMLも再読してるが、手動かしながら読むほうがしっくりくる。 p.70のベータ分布の逐次学習ってのが面白そうなのでやってみた。

パターン認識と機械学習上 - ベイズ理論による統計的予測
C. M. ビショップ
シュプリンガー・ジャパン株式会社 / ￥ 6,825 ()
在庫あり。

data <- rbinom(30,1,0.5)

plot_dbeta <- function(data){
a <- 0
b <- 0
x <- seq(0.01, 1.0, len = 500)
for (i in 1:length(data)) {
    if(data[i] == 1) { a = a + 1 } 
    else             { b = b + 1 }
    if(i>3){
    if(i==4){ plot(x, dbeta(x,a,b),type = "l",col=i,ylim=c(0,5)) }
    else    { plot(x, dbeta(x,a,b),type = "l",col=i,xlab="",ylab="",axes=F,ylim=c(0,5)) }
    par(new=T)
    }
}
}

二項分布で30個のリストを発生させて読みながらa,bのデータを更新していきつつその時のベータ分布の密度を求めつつプロットしていくという操作。

> data
 [1] 1 0 1 1 0 0 0 1 0 0 0 1 0 0 1 1 1 0 1 0 0 0 1 0 0 1 0 1 0
[30] 1

データはこんな感じで、4番目のデータからプロットしていくとだんだん二項分布に近づいていく様がみてとれる。

beta

06082009 R DMPK

ファーマコキネティクスの演習一通りやった

一通りやった。

ファーマコキネティクスを解いていく(Rで)

ファーマコキネティクス―演習による理解
杉山雄一
南山堂 / 6300円 ( 2003-08 )

2章のデコンボリューションはExcel依存だったのでスキップ。6章の薬物間相互作用と9章のトランスポーターは、今は必要なさそうなのでこれまたスキップ。

感想としてはクリアランス、吸収の章は理解が深まって良かった。PBPKモデリングとアロメトリーも同様。一方、PK-PDに関してはちょっと物足りなかった。あと、エクセルのマクロ(VBA?)依存の解答が多少目立ったので、普通はExcel必須かも。

一通りやってみて分かったことは、chemoinformaticsな自分のスタンスだと、PBPKモデリングがLeadOptimizationには親和性が高いかなというあたりと、Population PharmacokineticsよりはPK-PDをきちんと押さえないといかんなぁと。

ちゅうわけで、RでPK-PDをやる場合にはどの本読むのがいいんじゃろかね？

Pharmacokinetic-Pharmacodynamic Modeling And Simulation
Peter L., Ph.D. Bonate
Springer / 7596円 ( 2005-11 )

Pharmacometrics: The Science of Quantitative Pharmacology

Wiley-Interscience / 17858円 ( 2007-04-06 )

Dose Finding in Drug Development (Statistics for Biology and Health)

Springer / 11559円 ( 2006-03-07 )

Adaptive Design Methods in Clinical Trials (Chapman & Hall/CRC Biostatistics Series)
Shein-Chung Chow
Chapman and Hall/CRC / 8227円 ( 2006-11-16 )

06082009 R DMPK

PK8-4

Caco-2細胞のモデリング

ファーマコキネティクス―演習による理解
杉山雄一,山下伸二,加藤基浩
南山堂 / ￥ 6,300 ()
通常2～5週間以内に発送

設問1

Vmax は結合サイト数、ATPターンオーバー、細胞総数からもとまる

> 100000*1e+6*25/6.02e+23*60
[1] 2.491694e-10

つまり2.491694(nmol/min)

library(odesolve)
params = c(Vmax=2.491694*10^-3,S=4.71,Pa=5e-6,Pb=3e-6,Km=0.3,Va=1.5,Vc=0.025,Vb=3.0)
times=c(0,(1:60))

dydt <- function(t,y,p){
Vmax <- p['Vmax']
S <- p['S']
Pa <- p['Pa']
Pb <- p['Pb']
Km <- p['Km']
Va <- p['Va']
Vc <- p['Vc']
Vb <- p['Vb']
Ca <- (S*Pa*(y[2]-y[1]) + Vmax*y[2]/(Km+y[2]))/Va
Cc <- (S*(Pa*(y[1]-y[2]) + Pb*(y[3]-y[2])) - Vmax*y[2]/(Km+y[2]))/Vc
Cb <- (S*Pb*(y[2]-y[3]))/Vb
return(list(c(Ca,Cc,Cb)))
}

濃度リストをあたえるとPeffのリストを返す関数を書いておく。

calc <- function(dose_list){
l <- c()
for (dose in dose_list){
    y <- lsoda(c(Ca=dose,Cc=0.0,Cb=0.0),times,dydt,params)  
    Y <- data.frame(y)
    res <- lm(Cb ~ time,Y[c("time","Cb")][seq(41,61),])
    Peff <- res$coefficients[2]*3/(4.71*dose)
    l <- append(l,Peff)
}
return(l)
}

d <- c(0.1,0.3,0.6,1,3,6,10,30,60,100,300,600,900,1000)
peffs <- calc(d)
plot(d,peffs)

pk8-4

Pgpが飽和した後(高濃度側)ではシグモイド。

設問2

省略

05082009 R DMPK

PK8-3

膜透過性の解釈

ファーマコキネティクス―演習による理解
杉山雄一,山下伸二,加藤基浩
南山堂 / ￥ 6,300 ()
通常2～5週間以内に発送

> Mw = 350
> pKa = 9.5
> sigma = 1.0
> R = 8.62
> e_g = 1.22
> kai = 0.7
> Nav = 6.02 * 10^23
> kb = 1.38 * 10 ^ -16
> T = 273.15 + 37
> n = 0.006915
> r <- (3*Mw/(4*pi*Nav*sigma))^(1/3)
> r
[1] 5.17759e-08
> D <- kb*T/(6*pi*n*r)
> D
[1] 6.342054e-06

> r= 5.17759 # A
> F_r_R <- (1-(r/R))^2*(1-2.104*(r/R)+2.109*(r/R)^3-0.95*(r/R)^5)
> F_r_R
[1] 0.01897611
> r= 5.17759
> F_r_R <- (1-(r/R))^2*(1-2.104*(r/R)+2.109*(r/R)^3-0.95*(r/R)^5)
> F_r_R
[1] 0.0189761
> Ppara <- e_g*D*F_r_R*(kai/(1-exp(-kai)))
> Ppara
[1] 2.041593e-07
> 2.041593e-07/0.42e-6 * 100
[1] 48.60936 # %

05082009 R DMPK

PK8-2

経口投与量とAUCの関係

ファーマコキネティクス―演習による理解
杉山雄一,山下伸二,加藤基浩
南山堂 / ￥ 6,300 ()
通常2～5週間以内に発送

設問1

> data <- read.csv("/Users/kzfm/PK/pk8-2.csv")
> data$dose*1000/(data$AUCp_oral*0.8)
[1]  833.3333 1059.3220 2155.1724

設問2

CLtot(iv)は

> 10*1000/(45*0.8)
[1] 277.7778

題意よりCLh = CLtot(iv) さらにEh*Qh = CLhだから

> Eh <- 277.7778/1500
> Eh
[1] 0.1851852

> Foral <- 277.7778/833.3333
> Foral
[1] 0.3333334

> Fa <- Foral/(1-Eh)
> Fa
[1] 0.409091

設問3

Faで考える

もう5年目(wishlistありマス♡)
最近はPythonとDeepLearning
日本酒自粛中
ドラムンベースからミニマルまで
ポケモンGOゆるめ