追記
単に私がアホなだけでした（ごめんなさい）
端的に言えば、Stanでは共役事前分布をCauchy分布に設定していたにもかかわらず、Turingで$\sigma$の共役事前分布が設定できていなかったせいでした（Truncateで半コーシー分布や半正規分布を設定できるのを知らなかった）。
よってコードはこう変わる。
#σの引数がなくなった。
@model function prediction(x, y, n)
Hill~Uniform(0,1.5)
Treat_EC50 ~Uniform(-10.0,-4.0)
Non_EC50 ~Uniform(-10.0,-4.0)
σ ~ truncated(Cauchy(0, 2.5), 0, Inf)
#半コーシー分布を設定するには、Cauchy分布をCauchy（mu,σ）を設定し、正の値(>=0)になるようにtruncated()で囲む。
for i in 1:n
theta = 1/(1+10^(Hill*(Treat_EC50*x[i,2]+Non_EC50*x[i,3]-x[i,1])))
y[i] ~ Normal(theta,σ)
end
end;
すると、
# Sample using HMC.
#下の関数と比較して、sigmaの引数がなくなった。
m = prediction(Independ_res, Response, n)
chain = sample(m, HMC(0.05, 10),MCMCThreads(),5000,3)
describe(chain)
plot(chain)
chain
Chains MCMC chain (5000×13×3 Array{Float64, 3}):
Iterations = 1:1:5000
Number of chains = 3
Samples per chain = 5000
Wall duration = 13.75 seconds
Compute duration = 13.38 seconds
parameters = Hill, Treat_EC50, Non_EC50, σ
internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, step_size, nom_step_size
Summary Statistics
parameters mean std naive_se mcse ess rhat ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
Hill 0.8333 0.0476 0.0004 0.0006 6410.8810 1.0001 ⋯
Treat_EC50 -6.0693 0.0595 0.0005 0.0005 8782.5684 1.0001 ⋯
Non_EC50 -6.9137 0.0482 0.0004 0.0010 2103.2754 1.0028 ⋯
σ 0.0642 0.0376 0.0003 0.0006 4240.8766 1.0004 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
Hill 0.7457 0.8021 0.8319 0.8626 0.9295
Treat_EC50 -6.1582 -6.0994 -6.0685 -6.0389 -5.9790
Non_EC50 -7.0085 -6.9450 -6.9131 -6.8819 -6.8225
σ 0.0527 0.0591 0.0631 0.0675 0.0777
またNUTSをもちいると、
# Sample using HMC.
m = prediction(Independ_res, Response, n)
chain_NUTS = sample(m, NUTS(5000,0.65),MCMCThreads(),20000,10)
describe(chain_NUTS)
plot(chain_NUTS)
chain_NUTS
Chains MCMC chain (20000×16×10 Array{Float64, 3}):
Iterations = 5001:1:25000
Number of chains = 10
Samples per chain = 20000
Wall duration = 64.83 seconds
Compute duration = 63.31 seconds
parameters = Hill, Treat_EC50, Non_EC50, σ
internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
Summary Statistics
parameters mean std naive_se mcse ess rhat ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
Hill 0.8338 0.0476 0.0001 0.0001 256181.9135 1.0000 ⋯
Treat_EC50 -6.0685 0.0461 0.0001 0.0001 267290.2231 1.0000 ⋯
Non_EC50 -6.9134 0.0476 0.0001 0.0001 280346.1057 1.0000 ⋯
σ 0.0638 0.0064 0.0000 0.0000 239582.2881 1.0000 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
Hill 0.7452 0.8013 0.8320 0.8643 0.9326
Treat_EC50 -6.1591 -6.0992 -6.0686 -6.0378 -5.9780
Non_EC50 -7.0068 -6.9450 -6.9136 -6.8819 -6.8196
σ 0.0528 0.0593 0.0633 0.0677 0.0777
Stanでの推定結果とほぼ同じになった･･････。
#Non Treat群のEC50
# 0.5% 2.5% 25% 50% 75% 97.5% 99.5%
-7.03857 -7.00727 -6.94552 -6.91365 -6.88196 -6.82053 -6.78947
#Treat群のEC50
# 0.5% 2.5% 25% 50% 75% 97.5% 99.5%
-6.18968 -6.15888 -6.09906 -6.06839 -6.03783 -5.97773 -5.94764
#Hill係数
# 0.5% 2.5% 25% 50% 75% 97.5% 99.5%
0.7192710 0.7458398 0.8014040 0.8319190 0.8643612 0.9331391 0.9701001
はじめに
前回まで、最尤法（非線形最小二乗法）を用いてパラメータ推定を行った。
今回はパラメータ推定をベイズ推定法を用いて推定する。
用いるパッケージはCmdStanとTuringである。
CmdStanでやってみる（cmdstanr)
まずはパッケージを読み込む
library(cmdstanr)
library(tidyverse)
library(loo)
library(rstan)
library(bayesplot)
library(posterior)
library(tidybayes)
library(ggmcmc)
#前処理は裏で終えた結果がdat_all
dat_all2<-na.omit(dat_all)
conc<-dat_all2$Conc
#百分率の形になっているResponseを割合に変換することで、推定精度を高めている
Obs<-dat_all2$Response/100
Group<-as.integer(dat_all2$Group)
x_new<-seq(-10,-3.5,0.1)
x_new2<-rep(x_new,2)
k<-length(x_new2)
G_new<-rep(c(levels(dat_all2$Group)),each=length(x_new))
G_new<-factor(G_new, levels=unique(G_new))
Group_new<-as.integer(G_new)
l<-nrow(dat_all2)
n<-nlevels(dat_all2$Group)
\begin{align}
\theta_i&=\frac{1}{1+10^{(Hill(EC50[G_i]*G_i-Conc_i))}}\\
\\
y_i&=\theta_i+\sigma_i
\end{align}
$y_i:Response, \theta_i:用量反応曲線から導かれる予測値 y_{pred}, \sigma_i:誤差$
//stan読み込み
data{
int<lower=0> J;
int<lower=0> N;
int<lower=1> K;
array[J] real y;
array[J] real x;
array[J] int <lower=1,upper=N> g;
array[K] real X_new;
array[K] int <lower=1,upper=N> g_new;
}
//parameterの設定
parameters{
real<lower=0,upper=10> alpha;//Hill係数の範囲指定
array[N] real<lower=-10,upper=0> beta;//EC50の範囲指定
real sigma;
}
transformed parameters{
array[J] real<lower=0,upper=1> theta;
for(j in 1:J){
theta[j] = 1/(1+10^(alpha*(beta[g[j]]-x[j])));
//beta[g[j]]:各処理群毎にパラメータを推定
//仮にalpha(Hill係数)を各処理群毎に推定したい場合は、alpha[g[j]]となる
// 100/(1+10^^(alpha*(beta[g[j]]-x[j])))が1/(1+10^(alpha*(beta[g[j]]-x[j])))となっている
}
}
model{
for(j in 1:J)
y[j] ~ normal(theta[j], sigma);
sigma ~ cauchy(0,2.5);
}
//得られたパラメータを用いて、予測区間の推定
generated quantities {
array[K] real<lower=0,upper=1> theta_new;
array[K] real y_new;
for (k in 1:K){
theta_new[k]=1/(1+10^(alpha*(beta[g_new[k]]-X_new[k])));
y_new[k] = normal_rng(theta_new[k], sigma);
}
}
Bio <- list(J=l,
N=n,
K=k,
y=Obs,
x=conc,
g=Group,
X_new=x_new2,
g_new=G_new
)
nChains <- 10
nPost <- 1000 ## Number of post-burn-in samples per chain after thinning
nBurn <- 1000 ## Number of burn-in samples per chain after thinning
nThin <- 10
nIter <- (nPost + nBurn) * nThin
nBurnin <- nBurn * nThin/2
model<-cmdstan_model("上記のStancode読み込み")
#推定
#今回はChain10で、1Chainあたり25000ほど回す
fit <- model$sample(data=Bio,iter_warmup=nBurnin,
iter_sampling=nIter,thin<-nThin,chains = nChains,
adapt_delta=0.8,
max_treedepth = 20,
seed = 123)
fitの推定が終われば、
fit$summary()
fit$cmdstan_diagnose()
fitout<-rstan::read_stan_csv(fit$output_files())
mcmc_sample2<- fitout %>% rstan::extract()
# ggplotやtidyverseで扱いやすく処理する
mcmc_samples = as_draws_df(fit$draws())
# beta1のtrace plot
mcmc_samples %>%
mutate(chain = as.factor(.chain)) %>%
ggplot(aes(x = .iteration, y = `beta[1]`, group = .chain, color = chain)) +
geom_line()
mcmc_samples %>%
mutate(chain = as.factor(.chain)) %>%
ggplot(aes(x = .iteration, y = `beta[2]`, group = .chain, color = chain)) +
geom_line()
mcmc_samples %>%
mutate(chain = as.factor(.chain)) %>%
ggplot(aes(x = .iteration, y = `alpha`, group = .chain, color = chain)) +
geom_line()
Non_TreatのEC50(beta1)のTrace plot
TreatのEC50(beta2)のTrace plot
Hill係数のTrace plot
Bio.post <- data.frame(alpha = mcmc_samples$alpha,
beta1 = mcmc_samples$`beta[1]`,
beta2 = mcmc_samples$`beta[2]`)
Bio.post %>%
mutate(ld50 = beta1) %>%
with(quantile(ld50, c(0.005,.025, .25, .5, .75, .975,0.995)))
Bio.post %>%
mutate(ld50 = beta2) %>%
with(quantile(ld50, c(0.005,.025, .25, .5, .75, .975,0.995)))
Bio.post %>%
mutate(hill = alpha) %>%
with(quantile(hill, c(0.005,.025, .25, .5, .75, .975,0.995)))
結果
#Non Treat群のEC50
# 0.5% 2.5% 25% 50% 75% 97.5% 99.5%
-7.03857 -7.00727 -6.94552 -6.91365 -6.88196 -6.82053 -6.78947
#Treat群のEC50
# 0.5% 2.5% 25% 50% 75% 97.5% 99.5%
-6.18968 -6.15888 -6.09906 -6.06839 -6.03783 -5.97773 -5.94764
#Hill係数
# 0.5% 2.5% 25% 50% 75% 97.5% 99.5%
0.7192710 0.7458398 0.8014040 0.8319190 0.8643612 0.9331391 0.9701001
mcmc_samples %>%
ggplot() +
geom_histogram(aes(x=`beta[1]`),binwidth = 0.01)
mcmc_samples %>%
ggplot() +
geom_histogram(aes(x=`beta[2]`),binwidth = 0.01)
mcmc_samples %>%
ggplot() +
geom_histogram(aes(x=`alpha`),binwidth = 0.01)
rstan::stan_plot(fitout, pars = "beta", ci_level = 0.99, outer_level = 1, show_density = T,
point_est = "mean")
用量反応曲線の信頼区間、予測区間
dat_new<-bind_cols(x_new2,G_new,Group_new)%>% mutate(Conc=x_new2,Group=`...2`)
#mcmc_sample2のthetaの値
dat_out<-apply(mcmc_sample2$theta,2,quantile,probs=c(0.025, 0.50,0.975)) %>%
#転置してデータフレームに整形する
t() %>% bind_cols(dat_all2,.)
#y_newの取り出し
y.pred.interval <-
mcmc_sample2$y_new %>%
apply(MARGIN=2, quantile, prob=c(0.025, 0.5, 0.975)) %>%
t() %>%
as_tibble()
#作図
plot1<-bind_cols(dat_new,
y.pred.interval)%>% ggplot(aes(Conc, `50%`,colour=Group,fill=Group))+
geom_ribbon(aes(ymin=`2.5%`,y=`50%`,ymax=`97.5%`), alpha=0.1)+
geom_line()+
geom_ribbon(data=dat_out,aes(x=Conc,ymin=`2.5%`,ymax=`97.5%`), alpha=0.5)+#Responseの予測区間
geom_point(data=dat_all2, aes(Conc,Response/100,colour=Group))+theme_classic()
#Y軸の範囲指定
plot2<-plot1+
scale_y_continuous(limits = c(-0.2,1.25))
plot3<-plot2 + theme(axis.title = element_text(size = 15,
face = "bold"), axis.text = element_text(colour = "black"),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 15)) +labs(y = "Response")
#信頼区間の重ね合わせ（色の濃いところ）
plot3<-ggplot(data=dat_out,aes(x=Conc,y=`50%`,colour=Group,fill=Group))+
geom_ribbon(aes(ymin=`2.5%`,ymax=`97.5%`), alpha=0.1)
Turing.jlでやってみる
まずパッケージを読み込む
using Turing
using RCall
using StatsPlots
using DataFrames
using CSV
#データの読み込み、前処理
df=CSV.read("D:R-analysis/EC50/220814_dat_NA_Omit.csv",DataFrame)
dat_r2=rcopy(R"
G1<-as.numeric(factor($df$Group,levels=unique($df$Group)))-1
G2<-(-1)*(as.numeric(factor($df$Group,levels=unique($df$Group)))-2)
dat_r2<-data.frame(cbind($df,G_1=G1,G_2=G2))
"
)
dat_r3=dat_r2[:,2:7]
>print(dat_r3)
56×6 DataFrame
Row │ Conc Response Group Group_2 G_1 G_2
│ Float64 Float64 String Int64 Float64 Float64
─────┼────────────────────────────────────────────────────────────
1 │ -10.0 0.0 No_inhibitor 1 0.0 1.0
2 │ -8.0 2.831 No_inhibitor 1 0.0 1.0
3 │ -7.523 28.571 No_inhibitor 1 0.0 1.0
4 │ -7.0 45.238 No_inhibitor 1 0.0 1.0
5 │ -6.523 66.667 No_inhibitor 1 0.0 1.0
6 │ -6.0 76.19 No_inhibitor 1 0.0 1.0
7 │ -5.523 83.333 No_inhibitor 1 0.0 1.0
8 │ -5.0 92.857 No_inhibitor 1 0.0 1.0
9 │ -4.0 92.857 No_inhibitor 1 0.0 1.0
10 │ -10.0 0.0 No_inhibitor 1 0.0 1.0
11 │ -8.0 0.0 No_inhibitor 1 0.0 1.0
12 │ -7.523 23.81 No_inhibitor 1 0.0 1.0
13 │ -7.0 42.857 No_inhibitor 1 0.0 1.0
14 │ -6.523 71.429 No_inhibitor 1 0.0 1.0
15 │ -6.0 83.333 No_inhibitor 1 0.0 1.0
16 │ -5.523 97.619 No_inhibitor 1 0.0 1.0
17 │ -5.0 102.381 No_inhibitor 1 0.0 1.0
18 │ -4.0 102.381 No_inhibitor 1 0.0 1.0
19 │ -10.0 0.0 No_inhibitor 1 0.0 1.0
20 │ -8.0 9.524 No_inhibitor 1 0.0 1.0
21 │ -7.523 35.714 No_inhibitor 1 0.0 1.0
22 │ -7.0 52.381 No_inhibitor 1 0.0 1.0
23 │ -6.523 66.667 No_inhibitor 1 0.0 1.0
24 │ -6.0 83.333 No_inhibitor 1 0.0 1.0
25 │ -5.523 90.476 No_inhibitor 1 0.0 1.0
26 │ -5.0 104.762 No_inhibitor 1 0.0 1.0
27 │ -4.0 104.762 No_inhibitor 1 0.0 1.0
28 │ -10.0 0.0 Inhibitor 0 1.0 -0.0
29 │ -7.523 2.804 Inhibitor 0 1.0 -0.0
30 │ -7.0 14.019 Inhibitor 0 1.0 -0.0
31 │ -6.523 22.43 Inhibitor 0 1.0 -0.0
32 │ -6.0 47.664 Inhibitor 0 1.0 -0.0
33 │ -5.523 64.486 Inhibitor 0 1.0 -0.0
34 │ -5.0 75.701 Inhibitor 0 1.0 -0.0
35 │ -4.523 84.112 Inhibitor 0 1.0 -0.0
36 │ -4.0 89.72 Inhibitor 0 1.0 -0.0
37 │ -3.523 89.72 Inhibitor 0 1.0 -0.0
38 │ -10.0 0.0 Inhibitor 0 1.0 -0.0
39 │ -7.523 0.0 Inhibitor 0 1.0 -0.0
40 │ -7.0 8.411 Inhibitor 0 1.0 -0.0
41 │ -6.523 28.037 Inhibitor 0 1.0 -0.0
42 │ -6.0 56.075 Inhibitor 0 1.0 -0.0
43 │ -5.523 78.505 Inhibitor 0 1.0 -0.0
44 │ -5.0 86.916 Inhibitor 0 1.0 -0.0
45 │ -4.523 95.327 Inhibitor 0 1.0 -0.0
46 │ -4.0 103.738 Inhibitor 0 1.0 -0.0
47 │ -3.523 103.738 Inhibitor 0 1.0 -0.0
48 │ -10.0 0.0 Inhibitor 0 1.0 -0.0
49 │ -7.523 5.807 Inhibitor 0 1.0 -0.0
50 │ -7.0 16.822 Inhibitor 0 1.0 -0.0
51 │ -6.523 33.645 Inhibitor 0 1.0 -0.0
52 │ -6.0 64.486 Inhibitor 0 1.0 -0.0
53 │ -5.523 81.308 Inhibitor 0 1.0 -0.0
54 │ -5.0 95.327 Inhibitor 0 1.0 -0.0
55 │ -4.523 106.542 Inhibitor 0 1.0 -0.0
56 │ -4.0 106.542 Inhibitor 0 1.0 -0.0
Independ_res=dat_r3[:,["Conc","G_1","G_2"]]
Response=dar_r3[:,"Response"]./100
モデル式の設定
@model function prediction(x, y, n, σ)
Hill~Uniform(0,1.5)
Treat_EC50 ~Uniform(-10.0,-4.0)
Non_EC50 ~Uniform(-10.0,-4.0)
for i in 1:n
theta = 1/(1+10^(Hill*(Treat_EC50*x[i,2]+Non_EC50*x[i,3]-x[i,1])))
y[i] ~ Normal(theta,σ)
end
end;
n, _ = size(Indenpend_res)
# Sample using HMC.
m = prediction(Indenpend_res, Response, n, 0.5)
#今回はsamplingを5000に設定した
chain = sample(m, HMC(0.05, 10),MCMCThreads(),5000,3)
describe(chain)
#Traceplotとdensityplotを出力
plot(chain)
chain
Chains MCMC chain (5000×12×3 Array{Float64, 3}):
Iterations = 1:1:5000
Number of chains = 3
Samples per chain = 5000
Wall duration = 17.86 seconds
Compute duration = 15.39 seconds
parameters = Hill, Treat_EC50, Non_EC50
internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, step_size, nom_step_size
Summary Statistics
parameters mean std naive_se mcse ess rhat ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
Hill 0.8732 0.3060 0.0025 0.0106 707.1129 1.0028 ⋯
Treat_EC50 -6.0830 0.3851 0.0031 0.0030 15567.6102 0.9999 ⋯
Non_EC50 -6.9281 0.4127 0.0034 0.0033 17927.2660 1.0000 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
Hill 0.3552 0.6327 0.8449 1.1058 1.4601
Treat_EC50 -6.8455 -6.3344 -6.0814 -5.8292 -5.3279
Non_EC50 -7.7564 -7.1909 -6.9229 -6.6577 -6.1354
samplingを25,000にし、chainを10に増やしたらこうなる。
Chains MCMC chain (25000×12×10 Array{Float64, 3}):
Iterations = 1:1:25000
Number of chains = 10
Samples per chain = 25000
Wall duration = 211.89 seconds
Compute duration = 372.32 seconds
parameters = Hill, Treat_EC50, Non_EC50
internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, step_size, nom_step_size
Summary Statistics
parameters mean std naive_se mcse ess rhat ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
Hill 0.8736 0.3021 0.0006 0.0029 10878.7223 1.0005 ⋯
Treat_EC50 -6.0799 0.3910 0.0008 0.0008 230923.6258 1.0000 ⋯
Non_EC50 -6.9340 0.4095 0.0008 0.0008 244891.5996 1.0000 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
Hill 0.3543 0.6385 0.8501 1.1027 1.4504
Treat_EC50 -6.8605 -6.3305 -6.0768 -5.8251 -5.3169
Non_EC50 -7.7583 -7.1913 -6.9276 -6.6679 -6.1457
CmdStanでの推定結果と比べると、推定したパラメータの幅が広い。
#Non Treat群のEC50
# 0.5% 2.5% 25% 50% 75% 97.5% 99.5%
-7.03857 -7.00727 -6.94552 -6.91365 -6.88196 -6.82053 -6.78947
#Treat群のEC50
# 0.5% 2.5% 25% 50% 75% 97.5% 99.5%
-6.18968 -6.15888 -6.09906 -6.06839 -6.03783 -5.97773 -5.94764
#Hill係数
# 0.5% 2.5% 25% 50% 75% 97.5% 99.5%
0.7192710 0.7458398 0.8014040 0.8319190 0.8643612 0.9331391 0.9701001
またNUTSで推定した場合はこうなる。
# Sample using NUTS.
m = prediction(Independ_res, Response, n, 0.5)
chain_NUTS_v = sample(m, NUTS(0.65),MCMCThreads(),25000,10)
describe(chain_NUTS_v)
#Traceplotとdensityplotを出力
plot(chain_NUTS_v)
#結果を出力
chain_NUTS_v
Chains MCMC chain (25000×15×10 Array{Float64, 3}):
Iterations = 1001:1:26000
Number of chains = 10
Samples per chain = 25000
Wall duration = 167.21 seconds
Compute duration = 268.24 seconds
parameters = Hill, Treat_EC50, Non_EC50
internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
Summary Statistics
parameters mean std naive_se mcse ess rhat ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
Hill 0.8756 0.3019 0.0006 0.0009 134989.9569 1.0000 ⋯
Treat_EC50 -6.0798 0.3931 0.0008 0.0010 150105.8533 1.0000 ⋯
Non_EC50 -6.9342 0.4106 0.0008 0.0011 136846.8861 1.0001 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
Hill 0.3579 0.6403 0.8525 1.1057 1.4476
Treat_EC50 -6.8652 -6.3297 -6.0770 -5.8241 -5.3143
Non_EC50 -7.7631 -7.1914 -6.9284 -6.6677 -6.1432
NUTSで推定しても、Stanと比べて推定したパラメータの範囲が広い。
NUTS推定においても、Stanと比べて推定したパラメータの範囲が広くなった原因として、
「Stanではiter_warmupを設定することでウォームアップの設定が出来ていた。一方、Turingでは出来ていなかった。そのため、上記のような誤差につながった。」
と考えた。
そこで、adaptationをStanと同じ5000に設定し、Samplingを20,000の条件で行った。
m = prediction(Independ_res, Response, n, 0.5)
chain_NUTS= sample(m, NUTS(5000,0.65),MCMCThreads(),20000,10)
describe(chain_NUTS)
#結果を出力
chain_NUTS
Chains MCMC chain (20000×15×10 Array{Float64, 3}):
Iterations = 5001:1:25000
Number of chains = 10
Samples per chain = 20000
Wall duration = 87.18 seconds
Compute duration = 85.67 seconds
parameters = Hill, Treat_EC50, Non_EC50
internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
Summary Statistics
parameters mean std naive_se mcse ess rhat ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
Hill 0.8750 0.3005 0.0007 0.0009 111381.1466 1.0000 ⋯
Treat_EC50 -6.0797 0.3893 0.0009 0.0011 129129.2151 1.0000 ⋯
Non_EC50 -6.9344 0.4080 0.0009 0.0012 116135.2830 1.0001 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
Hill 0.3588 0.6411 0.8515 1.1044 1.4461
Treat_EC50 -6.8583 -6.3287 -6.0765 -5.8265 -5.3201
Non_EC50 -7.7561 -7.1930 -6.9283 -6.6693 -6.1473
しかしパラメータ推定の幅を狭めることが出来なかった。Adaptationの大きさは関係がなかったようである。
階層モデルで推定してみる
次に階層モデルを用いて推定した。
まずモデル式は以下のようになる。
theta(x) = 1 / (1 + 10^(x))
@model function prediction3(x, y,σ)
Hill~Uniform(0,1.5)
β ~ MvNormal([-7,-6],[1.0 -0.5; -0.5 1.0])
tmp = map(xj->Hill*(β[xj.EC50]-xj.Conc), x)
v=theta.(tmp)
for j in eachindex(x)
y[j] ~ Normal(v[x[j].i],σ)
end
end;
この処理は、$EC_{50}$の部分$\beta$を多変量正規分布であるMvNormalを設定し、予めfunction化しておいたthetaに無名関数を用いて値を代入し、推定している。
Independ_res2=[(i=i,Conc=dat_r3[i,:Conc],EC50=ifelse(dat_r3[i,:Group]=="Inhibitor",2,1)) for i in 1:nrow(dat_r3)]
# Sample using HMC.
m3 = prediction3(Independ_res2, Response,0.5)
chain3 = sample(m3, HMC(0.05, 10),MCMCThreads(),25000,10)
describe(chain3)
#TraceplotとDensityplotを出力
plot(chian3)
#結果を出力
chain3
Chains MCMC chain (25000×12×10 Array{Float64, 3}):
Iterations = 1:1:25000
Number of chains = 10
Samples per chain = 25000
Wall duration = 134.21 seconds
Compute duration = 233.95 seconds
parameters = Hill, β[1], β[2]
internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, step_size, nom_step_size
Summary Statistics
parameters mean std naive_se mcse ess rhat ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
Hill 0.8902 0.2988 0.0006 0.0027 12074.1407 1.0007 ⋯
β[1] -6.9344 0.3623 0.0007 0.0009 152435.4121 1.0001 ⋯
β[2] -6.0672 0.3493 0.0007 0.0008 172239.0647 1.0001 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
Hill 0.3764 0.6562 0.8683 1.1186 1.4520
β[1] -7.6583 -7.1690 -6.9317 -6.6959 -6.2277
β[2] -6.7556 -6.2967 -6.0658 -5.8377 -5.3799
また、NUTSではこうなる。
# Sample using NUTS.
m3 = prediction3(Independ_res, Response, n, 0.5)
chain_NUTS_v = sample(m3, NUTS(1000,0.8),MCMCThreads(),25000,10)
describe(chain_NUTS_v)
#TraceplotとDensityplotを出力
plot(chain_NUTS)
#結果を出力
chain_NUTS
Chains MCMC chain (25000×15×10 Array{Float64, 3}):
Iterations = 1001:1:26000
Number of chains = 10
Samples per chain = 25000
Wall duration = 111.99 seconds
Compute duration = 176.65 seconds
parameters = Hill, β[1], β[2]
internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
Summary Statistics
parameters mean std naive_se mcse ess rhat ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
Hill 0.8890 0.2984 0.0006 0.0009 138314.0588 1.0000 ⋯
β[1] -6.9345 0.3614 0.0007 0.0009 185196.5259 1.0001 ⋯
β[2] -6.0694 0.3489 0.0007 0.0008 185171.6209 1.0000 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
Hill 0.3759 0.6550 0.8678 1.1171 1.4495
β[1] -7.6540 -7.1692 -6.9314 -6.6957 -6.2303
β[2] -6.7602 -6.2987 -6.0681 -5.8385 -5.3844
追記2
追記する前の記事では。Turing.jlにおいて正しくパラメータ推定できない。そのため、本記事で正しくパラメータ推定出来ると思われるStanコードとTuringのコードを最後に示す。
Stan
#パッケージの読み込み
library(cmdstanr)
library(tidyverse)
library(loo)
library(rstan)
library(bayesplot)
library(posterior)
library(tidybayes)
library(ggmcmc)
Bio <- list(J=l,
N=n,
K=k,
y=Obs,
x=conc,
g=Group,
X_new=x_new2,
g_new=G_new
)
nChains <- 10
nPost <- 1000 ## Number of post-burn-in samples per chain after thinning
nBurn <- 1000 ## Number of burn-in samples per chain after thinning
nThin <- 10
nIter <- (nPost + nBurn) * nThin
nBurnin <- nBurn * nThin/2
#モデルの指定
model<-cmdstan_model("上記のStancode読み込み")
#推定
fit <- model$sample(data=Bio,iter_warmup=nBurnin,
iter_sampling=nIter,thin<-nThin,chains = nChains,
adapt_delta=0.8,
max_treedepth = 20,
seed = 123)
//stan読み込み
data{
int<lower=0> J;
int<lower=0> N;
int<lower=1> K;
array[J] real y;
array[J] real x;
array[J] int <lower=1,upper=N> g;
array[K] real X_new;
array[K] int <lower=1,upper=N> g_new;
}
//parameterの設定
parameters{
real<lower=0,upper=10> alpha;//Hill係数の範囲指定
array[N] real<lower=-10,upper=0> beta;//EC50の範囲指定
real sigma;
}
transformed parameters{
array[J] real<lower=0,upper=1> theta;
for(j in 1:J){
theta[j] = 1/(1+10^(alpha*(beta[g[j]]-x[j])));
//beta[g[j]]:各処理群毎にパラメータを推定
//仮にalpha(Hill係数)を各処理群毎に推定したい場合は、alpha[g[j]]となる
// 100/(1+10^^(alpha*(beta[g[j]]-x[j])))が1/(1+10^(alpha*(beta[g[j]]-x[j])))となっている
}
}
model{
for(j in 1:J)
y[j] ~ normal(theta[j], sigma);
sigma ~ cauchy(0,2.5);
}
//得られたパラメータを用いて、予測区間の推定
generated quantities {
array[K] real<lower=0,upper=1> theta_new;
array[K] real y_new;
for (k in 1:K){
theta_new[k]=1/(1+10^(alpha*(beta[g_new[k]]-X_new[k])));
y_new[k] = normal_rng(theta_new[k], sigma);
}
}
Turing
@model function prediction(x, y, n)
Hill~Uniform(0,1.5)
Treat_EC50 ~Uniform(-10.0,-4.0)
Non_EC50 ~Uniform(-10.0,-4.0)
σ ~ truncated(Cauchy(0, 2.5), 0, Inf)
#半コーシー分布を設定するには、Cauchy分布をCauchy（mu,σ）を設定し、正の値(>=0)になるようにtruncated()で囲む。
for i in 1:n
theta = 1/(1+10^(Hill*(Treat_EC50*x[i,2]+Non_EC50*x[i,3]-x[i,1])))
y[i] ~ Normal(theta,σ)
end
end;
# Sample using HMC.
#下の関数と比較して、sigmaの引数がなくなった。
m = prediction(Independ_res, Response, n)
chain = sample(m, HMC(0.05, 10),MCMCThreads(),5000,3)
describe(chain)
chain
plot(chain)
# Sample using NUTS
# Sample using NUTS.
m = prediction(Independ_res, Response, n)
chain_NUTS_v = sample(m, NUTS(0.65),MCMCThreads(),5000,3)
describe(chain_NUTS_v)
#Traceplotとdensityplotを出力
plot(chain_NUTS_v)
#結果を出力
chain_NUTS_v
階層モデルでは
theta(x) = 1 / (1 + 10^(x))
@model function prediction3(x, y,σ)
Hill~Uniform(0,1.5)
β ~ MvNormal([-7,-6],[1.0 -0.5; -0.5 1.0])
tmp = map(xj->Hill*(β[xj.EC50]-xj.Conc), x)
v=theta.(tmp)
for j in eachindex(x)
y[j] ~ Normal(v[x[j].i],σ)
end
end;
Independ_res2=[(i=i,Conc=dat_r3[i,:Conc],EC50=ifelse(dat_r3[i,:Group]=="Inhibitor",2,1)) for i in 1:nrow(dat_r3)]
# Sample using HMC.
m3 = prediction3(Independ_res2, Response,0.5)
chain3 = sample(m3, HMC(0.05, 10),MCMCThreads(),25000,10)
describe(chain3)
chain3
plot(chain3)
最後に
今回はCmdStan（CmdStanrから用いた）とJuliaのTuringを用いて、Hill係数、EC50を推定した。本来、パラメータ推定をGIGOでない推定をしたいのであれば、ブートストラップ法などするべきである。それに関しては、他の文献を読んで貰えるとよいと思う。
今回StanとTuringを比較した感想として、取り敢えず雑に推定するならTuringでの推定で良いが、より正確にパラメータをベイズ法で推定するにはStanを使った方が良いのではないかと思った。
みなさん。ベイズ統計モデリングをするときは、事前分布はきちんと設定しましょう！
半コーシー分布はTruncate($\mu$,$\sigma$,lower,upper)で設定できます。
参考文献
用量反応曲線のEC50の標準偏差を出す
StanとRでベイズ統計モデリング（Wonderful R）
松浦健太郎
政治学方法論Ⅱ（2015年度前期 神戸大学法学部・法学研究科）
矢内勇生
Turingで統計モデリング
Keisuke Yanagi
Turing.jl Tutorial