追記

単に私がアホなだけでした（ごめんなさい）

端的に言えば、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