単一事例実験データを変化点検出モデルで解析

変化点検出モデル

各フェーズ\(p\)(ベースライン期 or 介入期)における初めの観測値\(y_{p1}\)が、平均\(\hat{y}_{p1}\),分散\(\sigma^2_{\epsilon}\)の正規分布に従うとします。

\[y_{p1}\sim normal(\hat{y}_{p1},\sigma^2_\epsilon).\]

そして、以降の時点\(t\)におけるアウトカムの予測値は、

\[y_{pt}|H_{pt-1}, \Theta \sim normal(\hat{y}_{pt|(pt-1)},\sigma^2_e),\]

と正規分布に従うとします。

パラメータは下記の通りです。

• \(H_{pt-1}\): 過去の履歴
  
• \(\Theta\): パラメータベクトル
  
• \(\sigma^2_e\): ホワイトノイズ (ランダム誤差+自己相関)
• \(\sigma^2_\epsilon\): ランダム誤差
  
• \(\rho\):隣接時点間の自己相関

なお、\(\rho\),\(\sigma^2_e\)と\(\sigma^2_\epsilon\)の関係は、

\[\sigma_e=\frac{\sigma_\epsilon}{\sqrt {1-\rho^2}}, \]

で表現され、時系列データはラグ1の自己相関エラーを含む線形モデルとして下記のように表現可能となります。

変化点検出モデル+AR1(CP_AR1)のStan_code

参考論文内にJagsのコードが記載されていたので、それをStanに書き変えてみました。KMSコード(きっともっとシュッとするコード)で恐縮です。このコードでJagsの推定した結果とほぼ一致します。事前分布の設定も論文とあわせていますが、適宜変更しても良いでしょう。先に説明したモデルと、Stanで書いたモデルでパラメータ名が異なりますので念のため、再度書いておきます。わかりずらくてすみません。

• mu1: ベースラインの切片(平均)
  
• mu2: 介入期の切片 (平均)
  
• sigma: ホワイトノイズ
  
• rho: 自己相関
  
• es: 標準化効果量 (mu2-mu1/sigma)
  
• tau: 周辺化消去で求める変化点

data {
    int<lower=1>T;
    real y[T]; 
}

parameters {
    real mu1;
    real mu2;
    real<lower=0> sigma;
    real <lower=-1,upper=1> rho;
}

transformed parameters {
      vector[T] log_p;
      vector[T] x;
      vector[T] dummy;
      real mu;
      log_p = rep_vector(-log(T), T);

      for (tau in 1:T){
        for (i in 1:T){
          if (i == 1){
            x[1] = 0;
            dummy[1] = 1;
            log_p[tau] = log_p[tau] + normal_lpdf(y[1] | mu1, sigma);
            }else{
              dummy[i] = int_step(tau-i);
              x[i] = dummy[i] * mu1 + (1 - dummy[i]) * mu2;
              mu = i < tau ? x[i] : (x[i]+rho * (y[i - 1] - x[i - 1]));
              log_p[tau] = log_p[tau] + normal_lpdf(y[i] | mu, sigma);
              }
        }
    }
}



model {
    mu1 ~ normal(0,100);
    mu2 ~ normal(0,100);
    rho ~ uniform(-1,1);
    sigma ~uniform(0.1,5);
    target += log_sum_exp(log_p);

} 

//Draw the discrete parameter tau. This is highly inefficient
generated quantities {
    int<lower=1,upper=T> tau;
    real es;
    tau = categorical_rng(softmax(log_p));
    es = (mu2-mu1)/sigma;
}

このモデルで推定を行うRコードは下記の通りです。subject5,7,9の3名それぞれにモデルを適用しています。purrrのmap関数が便利です。

fit<-list()

# map関数でlist型で3つのデータの推定をする為に関数化

stan_list<-function(data){

y<-data

# 初期値に使うデータ
b1<- mean(y$outcome1[y$treatment=="baseline"],na.rm=T)
b2 <- mean(y$outcome1[y$treatment=="treatment"],na.rm=T)

time <- nrow(y)

# 初期値の指定  
initf1<- function(){
  list(mu1=rnorm(1, b1, 1),
       mu2=rnorm(1,b2,1), 
       sigma=runif(1,0.1,5),
       rho=runif(1,-1,1),
       tau=sample(3:(time-2),1))
  }
# stanによる推定の実行
rstan::sampling(CP_ar1,data=list(y=y$outcome1, T=length(y$outcome1)),init=initf1, 
         iter = 5000)

}


dat%>% filter(subject==5|subject==7|subject==9)  %>%  
  split(.$subject) %>% map(stan_list) -> fit

推定結果

収束診断の結果は割愛しますが、推定対象のパラメータはすべて\(\hat{R}\)が1.01以内で、問題なく収束していました。早速、変化点\(tau\)の推定結果を見てみましょう。

tau_plot<-function(data){
  data %>% spread_draws(tau)%>% 
    ggplot()+geom_histogram(aes(x=tau),fill="burlywood4",alpha=0.6)+theme_bw()+ 
    xlim(1,15)
}


fit %>% map(tau_plot) %>% arrangeGrob(grobs=.)->q1


q0_2<-list(q0$grobs$`5`,q0$grobs$`7`,q0$grobs$`9`, 
           q1$grobs$`5`,q1$grobs$`7`, q1$grobs$`9`)

grid.arrange(grobs=q0_2,ncol=3)

上段がローデータで、下段が変化点の推定結果です。subjetc5,9は推定された変化点が実際の介入導入時点に集中しています。一方で、subject7は、実際の介入導入時点と推定された変化点に乖離がみられます。subject7の変化点は12に集中しています。確かにローデータをみると、介入導入時点では大きな得点の変化がなく時点11までズルズルと得点が推移しています。ここでは時点11までは、ベースラインを引きずっていて、介入効果は12の時点で発現した解釈することにします。

ベースライン期と介入期の切片\(mu1,mu2\)とその差を標準偏差\(sigma\)で割った標準化効果量の推定結果をみてみましょう。標準化効果量は、generated quntitiesブロックで\(mu1,mu2,sigma\)から生成量として推定しています。

左が各フェーズの切片の事後分布、右が標準化効果量の事後分布です。上からsubject5,7,9の推定結果になっています。いずれの対象者においても、明らかに各切片の事後分布の重複が少なく、標準化効果量の推定結果も95%確信区間に0を含んでいません。この3名について、推定された変化点の前後で顕著な得点の変化があったといえます。

es_plot<-function(data){
  data%>% spread_draws(es) %>% 
    ggplot(aes(x=es,y=1,alpha=0.7))+ 
    geom_halfeyeh(fill="purple")+theme_bw()+ 
    theme(legend.position="none")
}

dif_plot<-function(data){
  data%>% gather_draws(mu1,mu2) %>% 
    ggplot(aes(x=.value,y=.variable,fill=.variable,alpha=.5))+ 
    geom_halfeyeh()+scale_fill_brewer(palette="Set1")+theme_bw()+ 
    theme(legend.position="none")
}


fit %>% map(es_plot) %>% arrangeGrob(grobs=.) ->q2
fit %>% map(dif_plot) %>% arrangeGrob(grobs=.) ->q3


grid.arrange(q3,q2,ncol=2)

最後に自己相関\(rho\)とホワイトノイズ\(sigma\)の事後分布も確認してみましょう。左が自己相関、右がホワイトノイズです。いずれも推定結果の不確実性が大きく、特に自己相関の不確実性が大きくなります。これは、測定時点数が15時点と少ないことが大きな要因の一つでしょう。

rho_plot<-function(data){
  data%>% spread_draws(rho) %>% 
    ggplot(aes(x=rho,y=1,alpha=0.6))+ 
    geom_halfeyeh(fill="darkslategrey")+theme_bw()+ 
    theme(legend.position="none")
}

sigma_plot<-function(data){
  data%>% spread_draws(sigma) %>% 
    ggplot(aes(x=sigma,y=1,alpha=0.6))+ 
    geom_halfeyeh(fill="darkslategrey")+theme_bw()+
    theme(legend.position="none")
}


fit %>% map(rho_plot) %>% arrangeGrob(grobs=.)->q4
fit %>% map(sigma_plot) %>% arrangeGrob(grobs=.)->q5

grid.arrange(q4,q5,ncol=2)

最後にsubject5,7,9それぞれの5つのパラメータの事後分布の要約統計量を示しておきます。

pi_table<-function(data){
  data %>%  gather_draws(rho,mu1,mu2,sigma,es) %>%  
    point_interval(.interval=hdi) 
}

fit %>% map(pi_table) %>% bind_rows() %>% mutate(subject=rep(c(5,7,9),each=5)) %>% kable(digits=2)

終わりに

以上のように単一事例実験デザインで得られたデータを、変化点検出モデルを使って解析してみました。2017年に提案されたとても新しいモデルで、今後適用例が増えていくと楽しいですね。個々人の変化点を推定したあとで、その変化点からベースラインと介入フェーズを指定し直して、階層モデルを使って複数名の時系列データから標準化効果量を算出する方法なども提案されています。単一事例実験の解析は、ここ数年がいくつかベイズ流のモデルが提案されてきていて興味深い限りです。それらについても別な機会に紹介できたらと思います。

Enjoy Stan !!

session info

sessionInfo()

## R version 3.5.0 (2018-04-23)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Sierra 10.12.6
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] ja_JP.UTF-8/ja_JP.UTF-8/ja_JP.UTF-8/C/ja_JP.UTF-8/ja_JP.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] bindrcpp_0.2.2      RColorBrewer_1.1-2  gridExtra_2.3      
##  [4] rstan_2.19.2        StanHeaders_2.19.0  tidybayes_1.0.3    
##  [7] forcats_0.3.0       stringr_1.4.0       dplyr_0.7.8        
## [10] purrr_0.2.5         readr_1.1.1         tidyr_0.8.2        
## [13] tibble_2.1.3        ggplot2_3.1.0       tidyverse_1.2.1    
## [16] zousan_0.0.0.170825 knitr_1.26         
## 
## loaded via a namespace (and not attached):
##  [1] httr_1.4.0                jsonlite_1.6             
##  [3] modelr_0.1.1              assertthat_0.2.1         
##  [5] highr_0.6                 stats4_3.5.0             
##  [7] ggstance_0.3.1            cellranger_1.1.0         
##  [9] yaml_2.2.0                pillar_1.4.2             
## [11] backports_1.1.5           lattice_0.20-35          
## [13] glue_1.3.1                arrayhelpers_1.0-20160527
## [15] digest_0.6.23             rvest_0.3.2              
## [17] colorspace_1.3-2          htmltools_0.4.0          
## [19] plyr_1.8.4                pkgconfig_2.0.3          
## [21] broom_0.5.0               svUnit_0.7-12            
## [23] HDInterval_0.2.0          haven_1.1.2              
## [25] bookdown_0.7              scales_1.0.0             
## [27] withr_2.1.2               lazyeval_0.2.1           
## [29] cli_1.1.0                 magrittr_1.5             
## [31] crayon_1.3.4              readxl_1.1.0             
## [33] evaluate_0.14             nlme_3.1-137             
## [35] MASS_7.3-49               xml2_1.2.0               
## [37] pkgbuild_1.0.2            blogdown_0.6             
## [39] tools_3.5.0               loo_2.0.0                
## [41] prettyunits_1.0.2         hms_0.4.2                
## [43] matrixStats_0.53.1        munsell_0.5.0            
## [45] callr_2.0.3               compiler_3.5.0           
## [47] rlang_0.4.2               debugme_1.1.0            
## [49] grid_3.5.0                RCurl_1.95-4.11          
## [51] ggridges_0.5.1            rstudioapi_0.7           
## [53] bitops_1.0-6              labeling_0.3             
## [55] rmarkdown_1.18            gtable_0.2.0             
## [57] codetools_0.2-15          inline_0.3.14            
## [59] R6_2.4.1                  lubridate_1.7.4          
## [61] bindr_0.1.1               stringi_1.4.3            
## [63] parallel_3.5.0            Rcpp_1.0.3               
## [65] tidyselect_0.2.5          xfun_0.11                
## [67] coda_0.19-1