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

(1)

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

1 パラメータ推定

1-3. データ解析の例

(2)

1-3. データ解析の例

1-3-1. ピアソン相関係数

1-3-2. ⼀致性のカッパ係数

(3)

1-3-1. ピアソン相関係数

「ベイズ統計で実践モデリング」

第5章データ解析の例

5.1 ピアソン係数

(4)

データはnペアの独⽴な観測値の対

例：特定の

薬剤の投与量

と投与からt 時間後の注⽬する

遺伝⼦の発現量

2つの変数間の線形の関係性は

ピアソンの積率相関係数r

で表現される

薬剤の投与量 t時間後の 注⽬する遺伝⼦の発現量

通常は相関係数rとその有意性（無相関の検定）

が記述される。

相関係数の事後分布を求めることで、⼀つの

数字で関係性を現すだけでなく、相関のレベル

がどの程度あるのか（信⽤区間）を求めること

ができる。

(5)

問題をグラフィカルモデルで表現

x

_i

i番⽬のデータ

μ

r

σ

μ

₁1

, σ

, μ

₂2

~ InvGamma(0.001, 0.001)

~ Gaussian(0, 0.001)

r ~ Uniform(-1, 1)

x

_i

~

MvGaussian

𝜇

_"

,𝜇

_$

,

𝜎

"$

𝑟𝜎

"

𝜎

$

𝑟𝜎

_"

𝜎

_$

𝜎

_$$ '"

(6)

尤度

i番⽬の投与量x

_i

[1]と、t時間後の発現量x

_i

[2]

x

_i

= (x

_i

[1], x

_i

[2])

t

平均ベクトル μ = (μ

₁

, μ

_２

)

t

分散共分散行列

Σ

₌

𝜎

"

$

_𝑟𝜎

"

𝜎

$

𝑟𝜎

_"

𝜎

_$

𝜎

_$

$

(

1 2𝜋

$

Σ

exp −

1

2 𝑥

2 − 𝜇

3 Σ

'"

𝑥

2 − 𝜇

4 25"

(7)

事前分布

μ

₁

, μ

₂

~ Gaussian(0, 0.001)

σ

₁

, σ

₂

~ InvGamma(0.001, 0.001)

r ~ Uniform(-1, 1)

後者は「７⼈の科学者」問題と同様、

λ=1/σ

2

_が

_{ガンマ分布}

_に

従うと仮定することに相当する。

(8)

Jagsでのモデルの記述

Correlation1.txt

# 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] <- rsigma[1]sigma[2]

T[2,1] <- rsigma[1]sigma[2]

T[2,2] <- 1/lambda[2]

TI[1:2,1:2] <- inverse(T[1:2,1:2])

}

(10)

R2jagsによるMCMCサンプリング

Correlation_1_jags.R

# clears workspace:

rm(list=ls())

# sets working directories:

setwd("C:/Users/EJ/Dropbox/EJ/temp/BayesBook/test/ParameterEstimation/DataAnalysis") library(R2jags) # 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) }

2種の観測データ

の⼀⽅を選択

(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.

(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)

(13)

> print(freq.r)

[1] -0.8109671

> summary(r)

V1

Min. :-0.9759

1st Qu.:-0.8184

Median :-0.7379

Mean :-0.7057

3rd Qu.:-0.6272

Max. : 0.1336

(14)

(15)

1-3-2. 一致性のカッパ係数

「ベイズ統計で実践モデリング」

第5章データ解析の例

5.2 ⼀致性のカッパ係数

(16)

カッパ係数

２つの⽅法が、どの程度同じ判断を与えるか

ただし判断は２値(1 or 0)

Poehling, Griffin and Dittus (2002)

233名の⼦供に対してインフルエンザの臨床検査

⽅法1と⽅法２の両⽅で感染しているとされたもの

14名 (1, 1)

⽅法１では感染していることとされたが⽅法２では診断できなかったもの

4名 (1, 0)

⽅法１では感染していないとされたが、⽅法２では感染しているとされたもの

5名 (0, 1)

⽅法１と⽅法２の両⽅で感染していないとされたもの

210名 (0, 0)

⽅法１は⾼価で、⽅法２は安価であるような場合に⽐較する

(17)

(1, 1)の数を a (1, 0)の数を b (0, 1)の数を c (0, 0)の数を d 観測された⼀致度

𝑝

₇

=

89:

4

偶然のみから期待される⼀致度

𝑝

;

=

89< 89= 9 <9: =9:

4

> ※ (a + b)/n : ⽅法１が１である頻度 (a + c)/n : ⽅法２が 1である頻度 (b + d)/n : ⽅法２が０である頻度 (c + d)/n : ⽅法１が０である頻度 a + b + c + d = n

カッパ係数 𝜅 =

A

B

'A

C

"'A

_C -1 (a=d=0, b=c=n/2)から１(a+d=n, b=c=0)までの値をとる

経験的には0.4未満を偶然に⽐べて芳しくなく、０.４〜０.７５は偶然よりもまあまあ良い、 0.75以上は偶然よりも優れた⼀致度されてきた

(18)

問題をグラフィカルモデルで表現

ξ

κ

_ψ

π

_a

π

_b

π

_c

π

_d

y

α

β

γ

yは観測データa, b, c, d

α ：⽅法１が１と判定する⽐率

(1-α）：⽅法１が 0と判定する⽐率

β ：⽅法１が１と判定した時に

⽅法２も１と判定する⽐率

γ ：⽅法１が０と判定する時に

⽅法２も０と判定する⽐率

(19)

問題をグラフィカルモデルで表現

ξ

κ

_ψ

π

_a

π

_b

π

_c

π

_d

y

α

β

γ

yは観測データa, b, c, d

α ：⽅法１が１と判定する⽐率

(1-α）：⽅法１が 0と判定する⽐率

β ：⽅法１が１と判定した時に

⽅法２も１と判定する⽐率

γ ：⽅法１が０と判定する時に

⽅法２も０と判定する⽐率

π

_a

=αβ、π

_b

=α（１ーβ）

π

_c

=(1-α)(1-γ）、π

_d

=(1-α)γ

(20)

問題をグラフィカルモデルで表現

ξ

κ

_ψ

π

_a

π

_b

π

_c

π

_d

y

α

β

γ

yは観測データa, b, c, d

α ：⽅法１が１と判定する⽐率

(1-α）：⽅法１が 0と判定する⽐率

β ：⽅法１が１と判定した時に

⽅法２も１と判定する⽐率

γ ：⽅法１が０と判定する時に

⽅法２も０と判定する⽐率

ξ＝αβ＋(1-α）γ ψ=

D

E

9D

F

D

_E

9D

_G

+

D

_F

9D

_H

D

_G

9D

_H

𝜅 =

𝜉 − 𝜓

1 − 𝜓

(21)

尤度

多項分布で尤度は表現される

𝑛!

(22)

事前分布

α, β, γ ~ Beta(1,1)

x <- seq(0, 1, 0.01)

(23)

Jagsでのモデルの記述

Kappa.txt

# 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 <- alphabeta+(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)

R2jagsによるMCMCサンプリング

Kappa_jags.R

# 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)

(27)

(28)

kappa <- samples$BUGSoutput$sims.list$kappa plot(kappa[,1],ty='l') summary(kappa) V1 Min. :0.3758 1st Qu.:0.6485 Median :0.7061 Mean :0.6998 3rd Qu.:0.7608 Max. :0.9171 quantile(kappa, c(0.025, 0.975)) 2.5% 97.5% 0.5125542 0.8481686

(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] 0.7357944

(30)

pi <- samples$BUGSoutput$sims.list$pi summary(pi[,1])

Min. 1st Qu. Median Mean 3rd Qu. Max. 0.02277 0.05033 0.05986 0.06094 0.07028 0.13221 summary(pi[,2])

Min. 1st Qu. Median Mean 3rd Qu. Max. 0.8075 0.8811 0.8947 0.8938 0.9077 0.9544

π

_a

, π

_b

, π

_c

, π

_d

の⽐較

14/(14+4+5+210) [1] 0.06008584 4/(14+4+5+210) [1] 0.01716738 5/(14+4+5+210) [1] 0.02145923 210/(14+4+5+210) [1] 0.9012876

(31)

1-3-3. 時系列データにおける変化検出

「ベイズ統計で実践モデリング」

第5章データ解析の例

5.3 時系列データにおける変化検出

(32)

時系列データにおける変化検出

近⾚外線分光データ

注意⽋陥多動性障害(ADHD)の成⼈の注意課題時の前頭葉の活動を酸化ヘモグロビンのカウ

ントによって測定

注意課題によってカウントの時系列に変化が⽣じる

この

変化点を同定

したい

1178点の時系列データ

観測データは時刻t

_i

でのカウントc

_i

である。

観測されていない変数τは変化が起こる時点

カウントは、ガウス分布に従う。分散は同じだが、

時点τで平均がμ

_１

からμ

_２

に変化する

(33)

問題をグラフィカルモデルで表現

i番⽬のサンプル

c

_i

t

i

τ

μ

₁

λ

μ

2

𝜇

_{𝜆 ~ Gamma(0.001, 0.001)}

"

, 𝜇

$

~ 𝐺𝑎𝑢𝑠𝑠𝑖𝑎𝑛 0,0.001

𝜏 ~ 𝑈𝑛𝑖𝑓𝑜𝑟𝑚(0, 𝑡

_b8c

)

𝑐

₂

~e

𝐺𝑎𝑢𝑠𝑠𝑖𝑎𝑛 𝜇

"

,𝜆 𝑖𝑓 𝑡

2

< 𝜏

𝐺𝑎𝑢𝑠𝑠𝑖𝑎𝑛 𝜇

_$

,𝜆 𝑖𝑓 𝑡

₂

≥ 𝜏

(34)

尤度

( 𝑑𝑛𝑜𝑟𝑚 𝑥

₂

, 𝜇

_"

, 𝜆

3

_h

ij

( 𝑑𝑛𝑜𝑟𝑚 𝑥

₂

, 𝜇

_$

, 𝜆

(35)

事前分布

𝜇

_"

, 𝜇

_$

~ 𝐺𝑎𝑢𝑠𝑠𝑖𝑎𝑛 0, 0.001

𝜆 ~ Gamma(0.001, 0.001)

𝜏 ~ 𝑈𝑛𝑖𝑓𝑜𝑟𝑚(0, 𝑡

_b8c

)

(36)

Jagsでのモデルの記述

ChangeDetection.txt

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以上であれば１を返す関数

配列zは0か１が⼊る

配列z1でこれに１を⾜して１か２が⼊る形にし、

これはデータが従うガウス分布の平均の添字にして、

平均の時点τでの切り替えを表現している

(38)

R2jagsによるMCMCサンプリング

ChangeDetection_jags.R

# 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(

(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)

(40)

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

(41)

samples <- jags(data, inits=myinits, parameters,

(42)

(43)

tau <- samples$BUGSoutput$sims.list$tau

summary(tau)

V1

Min. :726.2

1st Qu.:731.0

Median :731.6

Mean :732.3

3rd 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.

37.20 37.72 37.89 37.89 38.10 38.66

summary(mu[,2])

Min. 1st Qu. Median Mean 3rd Qu. Max.

29.68 30.38 30.59 30.59 30.79 31.51

(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 = 17.861, df = 1176, p-value < 2.2e-16

alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval:

6.528340 8.139576 sample estimates: mean of x mean of y

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