バイオインフォマティクス特論4

Size: px

Start display at page:

Download "バイオインフォマティクス特論4"

きよあつしんまつ
5 years ago
Views:

1 藤博幸

2 ピアソン相関係数致性のカッパ係数時系列データにおける変化検出

3 ベイズ統計で実践モデリング 5.1 ピアソン係数第 5 章データ解析の例

4 データは n ペアの独な観測値の対例 : 特定の薬剤の投与量と投与から t 時間後の注する遺伝の発現量 2 つの変数間の線形の関係性はピアソンの積率相関係数 r で表現される t 時間後の注する遺伝の発現量薬剤の投与量通常は相関係数 r とその有意性 ( 無相関の検定 ) が記述される相関係数の事後分布を求めることでつの数字で関係性を現すだけでなく相関のレベルがどの程度あるのか ( 信区間 ) を求めることができる

5 μ r σ x i μ 1, μ 2 ~ Gaussian(0, 0.001) σ 1, σ 2 ~ InvGamma(0.001, 0.001) r ~ Uniform(-1, 1) x i ~ MvGaussian μ ",μ $, $ σ " rσ " σ $ rσ " σ $ $ σ $ '" i 番のデータ

6 i 番の投与量 x i [1] と t 時間後の発現量 x i [2] x i = (x i [1], x i [2]) t 平均ベクトル μ = (μ 1, μ 2 ) t 分散共分散行列 Σ = $ σ " rσ " σ $ rσ " σ $ $ σ $ 4 ( 1 25" 2π $ Σ exp 1 2 x 2 μ 3 Σ '" x 2 μ

7 μ 1, μ 2 ~ Gaussian(0, 0.001) σ 1, σ 2 ~ InvGamma(0.001, 0.001) r ~ Uniform(-1, 1) 後者は 7 の科学者問題と同様 λ=1/σ 2 がガンマ分布に従うと仮定することに相当する

8 # Pearson Correlation model{ # Data for (i in 1:n){ x[i,1:2] ~ dmnorm(mu[],ti[,]) } # Priors mu[1] ~ dnorm(0,.001) mu[2] ~ dnorm(0,.001) lambda[1] ~ dgamma(.001,.001) lambda[2] ~ dgamma(.001,.001) r ~ dunif(-1,1)

9 # Reparameterization sigma[1] <- 1/sqrt(lambda[1]) sigma[2] <- 1/sqrt(lambda[2]) T[1,1] <- 1/lambda[1] T[1,2] <- r*sigma[1]*sigma[2] T[2,1] <- r*sigma[1]*sigma[2] T[2,2] <- 1/lambda[2] TI[1:2,1:2] <- inverse(t[1:2,1:2]) }

10 # clears workspace: rm(list=ls()) # sets working directories: setwd("c:/users/ej/dropbox/ej/temp/bayesbook/test/parameterestimation/dataanalysis") library(r2jags) 2 種の観測データのを選択 # Choose a dataset: dataset <- 1 # The datasets: if (dataset == 1) { x <- matrix(c(.8,102, 1,98,.5,100, 0.9,105,.7,103, 0.4,110, 1.2,99, 1.4,87, 0.6,113, 1.1,89, 1.3,93), nrow=11,ncol=2,byrow=t) } if (dataset == 2) { x <- matrix(c(.8,102, 1,98,.5,100, 0.9,105,.7,103, 0.4,110, 1.2,99, 1.4,87, 0.6,113, 1.1,89, 1.3,93,.8,102, 1,98,.5,100, 0.9,105,.7,103, 0.4,110, 1.2,99, 1.4,87, 0.6,113, 1.1,89, 1.3,93), nrow=22,ncol=2,byrow=t) }

11 n <- nrow(x) # number of people/units measured data <- list("x", "n") # to be passed on to JAGS myinits <- list( list(r = 0, mu = c(0,0), lambda = c(1,1))) # parameters to be monitored: parameters <- c("r", "mu", "sigma") # The following command calls JAGS with specific options. # For a detailed description see the R2jags documentation. samples = jags(data, inits=myinits, parameters, model.file ="Correlation_1.txt", n.chains=1, n.iter=5000, n.burnin=1, n.thin=1, DIC=T) # Now the values for the monitored parameters are in the "samples" object, # ready for inspection. r <- samples$bugsoutput$sims.list$r

12 #Frequentist point-estimate of r: freq.r <- cor(x[,1],x[,2]) layout(matrix(c(1,2),1,2)) layout.show(2) #some plotting options to make things look better: par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5, font.lab = 2, cex.axis = 1.3, bty = "n", las=1) # data panel: plot(x[,1],x[,2], type="p", pch=19, cex=1) # correlation panel: plot(density(r, from=-1,to=1), main="", ylab="posterior Density", xlab="correlation", lwd=2) lines(c(freq.r, freq.r), c(0,100), lwd=2, lty=2)

13 > print(freq.r) [1] > summary(r) V1 Min. : st Qu.: Median : Mean : rd Qu.: Max. :

15 ベイズ統計で実践モデリング 5.2 致性のカッパ係数第 5 章データ解析の例

16 カッパ係数 2 つの法がどの程度同じ判断を与えるかただし判断は 2 値 (1 or 0) Poehling, Griffin and Dittus (2002) 233 名の供に対してインフルエンザの臨床検査法 1と法 2の両で感染しているとされたもの 14 名 (1, 1) 法 1では感染していることとされたが法 2では診断できなかったもの 4 名 (1, 0) 法 1では感染していないとされたが法 2では感染しているとされたもの 5 名 (0, 1) 法 1と法 2の両で感染していないとされたもの 210 名 (0, 0) 法 1 は価で法 2 は安価であるような場合に較する

17 (1, 1) の数を a (1, 0) の数を b (0, 1) の数を c (0, 0) の数を d 観測された致度 p 7 = 89: 4 偶然のみから期待される致度 p ; = 89< 89= 9 <9: =9: 4 > (a + b)/n : 法 1 が 1 である頻度 (a + c)/n : 法 2 が 1 である頻度 (b + d)/n : 法 2 が 0 である頻度 (c + d)/n : 法 1 が 0 である頻度 a + b + c + d = n カッパ係数 κ = A B 'A C "'A C -1 (a=d=0, b=c=n/2) から 1(a+d=n, b=c=0) までの値をとる経験的には 0.4 未満を偶然にべて芳しくなくは偶然よりもまあまあ良い 0.75 以上は偶然よりも優れた致度されてきた

18 ξ κ ψ β α γ α : 法 1 が 1 と判定する率 (1-α): 法 1 が 0 と判定する率 β : 法 1 が 1 と判定した時に法 2 も 1 と判定する率 γ : 法 1 が 0 と判定する時に法 2 も 0 と判定する率 π a π b π c π d y y は観測データ a, b, c, d

19 ξ κ ψ β α γ π a π b π c π d α : 法 1 が 1 と判定する率 (1-α): 法 1 が 0 と判定する率 β : 法 1 が 1 と判定した時に法 2 も 1 と判定する率 γ : 法 1 が 0 と判定する時に法 2 も 0 と判定する率 π a =αβ π b =α(1 ー β) π c =(1-α)(1-γ) π d =(1-α)γ y y は観測データ a, b, c, d

20 ξ κ ψ β α γ π a π b π c π d ξ=αβ+(1-α)γ ψ= D E 9D F D E 9D G + D F 9D H D G 9D H ξ ψ κ = 1 ψ α : 法 1 が 1 と判定する率 (1-α): 法 1 が 0 と判定する率 β : 法 1 が 1 と判定した時に法 2 も 1 と判定する率 γ : 法 1 が 0 と判定する時に法 2 も 0 と判定する率 y y は観測データ a, b, c, d

21 多項分布で尤度は表現される n! a! b! c! d! π "" 8 π < "Q π = : Q" π QQ

22 α, β, γ ~ Beta(1,1) x <- seq(0, 1, 0.01) plot(x, dbeta(x,1,1), ty='l')

23 # Kappa Coefficient of Agreement model{ # Underlying Rates # Rate Objective Method Decides "one" alpha ~ dbeta(1,1) # Rate Surrogate Method Decides "one" When Objective Method Decides "one" beta ~ dbeta(1,1) # Rate Surrogate Method Decides "zero" When Objective Method Decides "zero" gamma ~ dbeta(1,1) # Probabilities For Each Count pi[1] <- alpha*beta pi[2] <- alpha*(1-beta) pi[3] <- (1-alpha)*(1-gamma) pi[4] <- (1-alpha)*gamma

24 } # Count Data y[1:4] ~ dmulti(pi[],n) # Derived Measures # Rate Surrogate Method Agrees With the Objective Method xi <- alpha*beta+(1-alpha)*gamma # Rate of Chance Agreement psi <- (pi[1]+pi[2])*(pi[1]+pi[3])+(pi[2]+pi[4])*(pi[3]+pi[4]) # Chance-Corrected Agreement kappa <- (xi-psi)/(1-psi)

25 # clears workspace: rm(list=ls()) # sets working directories: setwd("c:/users/ej/dropbox/ej/temp/bayesbook/test/parameterestimation/dataa nalysis") library(r2jags) # choose a data set: # Influenza y <- c(14, 4, 5, 210) n <- sum(y) # number of people/units measured data <- list("y", "n") # to be passed on to JAGS

26 myinits <- list( list(alpha = 0.5, beta = 0.5, gamma = 0.5)) # parameters to be monitored: parameters <- c("kappa","xi","psi","alpha","beta","gamma","pi") # The following command calls JAGS with specific options. # For a detailed description see the R2jags documentation. samples <- jags(data, inits=myinits, parameters, model.file ="Kappa.txt", n.chains=1, n.iter=2000, n.burnin=1, n.thin=1, DIC=T) # Now the values for the monitored parameters are in the "samples" object, # ready for inspection. plot(samples)

kappa <- samples$bugsoutput$sims.list$kappa plot(kappa[,1],ty='l') summary(kappa) V1 Min. :0.3758 1st Qu.:0.6485 Median :0.

28 kappa <- samples$bugsoutput$sims.list$kappa plot(kappa[,1],ty='l') summary(kappa) V1 Min. : st Qu.: Median : Mean : rd Qu.: Max. : quantile(kappa, c(0.025, 0.975)) 2.5% 97.5%

29 Cohen の κ の点推定値 p0 <- (y[1]+y[4])/n pe <- (((y[1]+y[2]) * (y[1]+y[3])) + ((y[2]+y[4]) * (y[3]+y[4]))) / n^2 kappa.cohen <- (p0-pe) / (1-pe) kappa.cohen [1]

30 π a, π b, π c, π d の較 pi <- samples$bugsoutput$sims.list$pi summary(pi[,1]) Min. 1st Qu. Median Mean 3rd Qu. Max summary(pi[,2]) Min. 1st Qu. Median Mean 3rd Qu. Max summary(pi[,3]) Min. 1st Qu. Median Mean 3rd Qu. Max /( ) [1] /( ) [1] /( ) [1] /( ) [1] summary(pi[,4]) Min. 1st Qu. Median Mean 3rd Qu. Max

31 ベイズ統計で実践モデリング第 5 章データ解析の例 5.3 時系列データにおける変化検出

32 近外線分光データ注意陥多動性障害 (ADHD) の成の注意課題時の前頭葉の活動を酸化ヘモグロビンのカウントによって測定注意課題によってカウントの時系列に変化がじるこの変化点を同定したい 1178 点の時系列データ観測データは時刻 t i でのカウント c i である観測されていない変数 τ は変化が起こる時点カウントはガウス分布に従う分散は同じだが時点 τ で平均が μ 1 から μ 2 に変化する

33 μ μ 1 λ 2 μ ", μ $ ~ Gaussian 0,0.001 λ ~ Gamma(0.001, 0.001) τ ~ Uniform(0, t b8c ) c 2 ~e Gaussian μ ",λ if t 2 < τ Gaussian μ $,λ if t 2 τ τ c i t i i 番のサンプル

34 ( dnorm x 2, μ ", λ 3 h ij ( dnorm x 2, μ $, λ 3 h kj

35 μ ", μ $ ~ Gaussian 0, λ ~ Gamma(0.001, 0.001) τ ~ Uniform(0, t b8c )

36 model{ # Data Come From A Gaussian for (i in 1:n){ c[i] ~ dnorm(mu[z1[i]],lambda) } # Group Means mu[1] ~ dnorm(0,.001) mu[2] ~ dnorm(0,.001) # Common Precision lambda ~ dgamma(.001,.001) sigma <- 1/sqrt(lambda) # Which Side is Time of Change Point? for (i in 1:n){ z[i] <- step(t[i]-tau) z1[i] <- z[i]+1 } # Prior On Change Point tau ~ dunif(0,n) }

37 step 関数引数が負ならば 0 0 以上であれば 1 を返す関数配列 z は 0 か 1 がる配列 z1 でこれに 1 をして 1 か 2 がる形にしこれはデータが従うガウス分布の平均の添字にして平均の時点 τ での切り替えを表現している

38 # clears workspace: rm(list=ls()) # sets working directories: setwd("/users/toh/desktop/code/parameterestimation/dataanalysis") library(r2jags) c <- scan("changepointdata.txt") n <- length(c) t <- 1:n data <- list("c", "n", "t") # to be passed on to JAGS myinits <- list( list(mu = c(1,1), lambda = 1, tau = n/2))

39 # parameters to be monitored: parameters <- c("mu","sigma","tau") # The following command calls JAGS with specific options. # For a detailed description see the R2jags documentation. samples <- jags(data, inits=myinits, parameters, model.file ="ChangeDetection.txt", n.chains=1, n.iter=1000, n.burnin=1, n.thin=1, DIC=T) # Now the values for the monitored parameters are in the "samples" object, # ready for inspection. traceplot(samples)

40 最初の200ステップは burn inとして捨てるべき

41 samples <- jags(data, inits=myinits, parameters, model.file ="ChangeDetection.txt", n.chains=1, n.iter=1000, n.burnin=1, n.thin=1, DIC=T) samples <- jags(data, inits=myinits, parameters, model.file ="ChangeDetection.txt", n.chains=1, n.iter=1000, n.burnin=200, n.thin=1, DIC=T) n.burnin の数値を変更しこの部分だけをコピペして MCMC を実する

42 plot(samples)

43 tau <- samples$bugsoutput$sims.list$tau summary(tau) V1 Min. : st Qu.:731.0 Median :731.6 Mean : rd Qu.:733.7 Max. :743.1 時点 732 あたりで切り替わっている

44 mu <- samples$bugsoutput$sims.list$mu summary(mu[,1]) Min. 1st Qu. Median Mean 3rd Qu. Max summary(mu[,2]) Min. 1st Qu. Median Mean 3rd Qu. Max

45 τ=732 の前後のカウント数で t 検定 t.test(c[1:731],c[732:1178],var.equal=t) Two Sample t-test data: c[1:731] and c[732:1178] t = , df = 1176, p-value < 2.2e-16 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: sample estimates: mean of x mean of y

バイオインフォマティクス特論12

バイオインフォマティクス特論12 藤博幸事後予測分布パラメータの事後分布に従ってモデルがどんなデータを期待するかを予測する予測分布が観測されたデータと致するかを確認することでモデルの適切さを確認できる前回と同じ問題で事後予測をう 3-1-1. 個差を考えない場合 3-1-2. 完全な個差を考える場合 3-1-3. 構造化された個差を考える場合ベイズ統計で実践モデリング 10.1 個差を考えない場合第 10